Molecular Dynamics Simulation

We have used molecular dynamics simulations extensively in this project. Our entire molecular dynamics effort was centered around investigating two approaches:

1. Are the NeoFvs stable at the site of action?
2. How strong are all the FcRn binding peptide variants binding to the FcRn receptor? Are any of these better than others?

NeoFv Stability

• What are we investigating?
  We are trying to investigate whether our NeoFv maintains structural stability in its condition at the site of action, i.e., at a pH of 7.4 in blood and at a temperature of 40 degrees celsius (the body temperature of a person with severe Dengue). We ran a trajectory for 100 ns on an AlphaFold predicted structure and analyzed the RMSD of the molecule over time.
• What are we investigating?
  Molecular dynamics simulations are highly realistic simulations at the molecular scale. They can be used to identify dynamical properties of the simulated system.
• What can this tell us?
  Analyzing the RMSD of the molecule can tell us whether the molecule is stable or not in the conditions simulated.

Technicals

One of the requirements of the therapeutic drug is that it should be stable. Our NeoFvs are meant to be a therapy for Dengue fever by binding and neutralizing the Dengue virus in the bloodstream (which has a pH of 7.4). Often, our proposed drug may be administered to people undergoing severe Dengue fever, where the body temperatures may reach 40 degrees celsius.
We incorporated these facts while simulating a molecular dynamics trajectory for checking the stability of our NeoFv.
We used the GROMACS package for running all our MDs.
We converted PDB structures of our NeoFvs (generated through AlphaFold 2.0.) into gro files, We went through the following steps before running the production MD.

1. Convert the .pdb file to a .gro file.
2. Choose the OPLS-AA/M force field
3. Define a cubic box around the molecule
4. Solve the box with water.
5. Add ions to the solvated box
6. Minimize the energy of the entire box
7. Equilibrate NVT
8. Equilibrate NPT
9. Generate an MD.tpr file

Once these steps are completed, we can run the molecular dynamics simulation.

The MD trajectories were calculated using IISER Pune’s PARAM Brahma supercomputer. We received the results within 24 hours. We generated room mean square deviation (RMSD) plots for the different MDs we ran.

For the 100 ns duration, RMSDs for all the variants are below 2Å, which is the accepted value for a molecule to be called stable in molecular dynamics.

Conclusion We verified that our NeoFv variants are stable under the conditions simulated.

FcRn Binding Peptide Complex Stability and pH-dependent binding

• What are we investigating?
  We are trying to estimate the association binding energies of the FcRn-binding peptide complex. We are simulating the MD at pH of the endosome where the interaction is supposed to occur. These will be analyzed by an end-state free energy calculation program: gmx_MMPBSA^[1]. We also aim to evaluate whether introducing a histidine residue in the consensus sequence has an impact on the binding energy.
• What are we investigating?
  Molecular docking scores can tell us very little about the stability of the sampled complex. Molecular dynamics approaches, on the other hand.
• What can this tell us?
  Molecular dynamics on complexes tell us whether the complex is stable in time or not.
  Binding free energies calculated from the MDs can help us compare all the different variants of our NeoFvs.
  It can tell us whether the folded version of the linear variant stays stable in the complex, as hypothesized in the structure predictions page.

Technicals

We again used the GROMACS package to input the complex MDs.

Preparing the protein and ligand complex

Since there was already a PDB structure of an FcRn inhibitory peptide similar to ours, and we had identified the sensitive residues to be the same as well, we decided to superimpose our peptide with the existing ligand protein complex.
For the AlphaFold prediction of the linear peptide variants, not all of our predictions were accurate, as discussed in the docking page, and therefore introduced mutations to the PDB 3M17 ligand and energy minimized the FcRnBP structure through OpenBabel to be a reliable starting point. The FcRnBPs were then complexed with the FcRn receptor through superimposition with the crystal structure PDB: 3M17.
Once the ligand-protein complexes were ready, we put the structures through PropKa. PropKa is a tool to predict protonation states and charges of individual residues in the protein based on the local electrostatic field and the defined pH.
Since the Y12H mutation was to improve the pH dependence of the FcRn interaction at pH 6, our intention was to verify whether there was any abrogation of binding at pH 7.4.
We ran complexed molecular dynamics trajectories for both Y12H (CycY12H and LinY12H) mutants at pH 6.0 and 7.4.

Results and analysis with gmx_MMPBSA

Free energy calculations on the different variants of the peptide were carried out to assess the variability associated with them through gmx_MMPBSA.

Variant	pH	Mean Binding free energy (kcal/mol)	Standard Deviation
Cyc	6	-32.80	4.22
CycY12H	6	-34.06	5.43
Lin	6	-37.11	6.91
LinY12H	6	-32.20	3.44
CycY12H	7.4	-31.06	3.99
LinY12H	7.4	-24.00	7.51

Most of the results we obtained from the MD simulation were around the same range of values, and with high standard deviation. Therefore, no conclusive predictions could be generated from these. However, as expected, the MDs confirm that the different variants of the peptide bind to the FcRn complex as indicated by the negative binding free energies associated with the complex.

Plausible explanation and future work

Our results stand in opposition to similar laboratory results in literature, where the introduction of histidine, in place of tyrosine (aka the Y12H mutation) abrogates binding at pH 7.4 by the creation of a pH dependent salt bridge by histidine. We could explain this by the differential protonation states of the histidine residue at different pHs.
Protonated histidines can have three variants HIE, HID and HIP, where HIP has two protons associated with it. Most of our MD simulations were run with HID, as predicted by PropKa, which works by analyzing the local electrostatic field.
However, it is possible that the sensitive histidine on the consensus sequence may be protonated to HIP in solution when it’s free. Thus, future MDs on this should properly investigate the protonation states of histidines involved and whether there is significant influence of histidine protonation states on the binding energy of the complex.
Another way to properly estimate binding energy of a complexed structure is through metadynamics and continuous addition of gaussian energy potentials across the molecule. This can help understand the depth of the minima and hence the binding energy.^[2]

Conclusions

1. Binding energies of FcRn binding peptides to the FcRn receptor were negative.
2. There was no significant difference in binding energies of the pH dependent mutants at different pHs.
3. The linear variants in the folded conformational state remained stable throughout the duration of the MD simulation.

References

1. gmx_MMPBSA: A New Tool to Perform End-State Free Energy Calculations with GROMACS Mario S. Valdés-Tresanco, Mario E. Valdés-Tresanco, Pedro A. Valiente, and Ernesto Moreno Journal of Chemical Theory and Computation 2021 17 (10), 6281-6291 DOI
   
2. Laio, A., & Parrinello, M. (2002). Escaping free-energy minima. Proceedings of the National Academy of Sciences of the United States of America, 99(20), 12562–12566. DOI