Molecular dynamics simulation analysis of structural dynamic cross correlation induced by odorant hydrogen-bonding in mouse eugenol ol- factory receptor

Abstract

Structural fluctuations and dynamic cross-correlations in the mouse eugenol olfactory receptor (Olfr73) were studied by molecular dynamics (MD) simulation to characterize the dynamic response of the protein upon ligand binding. The initial structure was generated by the artificial intelligence tool AlphaFold2 due to the current lack of experimental data. We focused on the hydrogen (H) bond of the odorant eugenol to Ser113, Asn207, and Tyr260 of the receptor protein, the importance of which has been suggested by previous experimental studies. The H-bond was not observed in docking simulations, but in subsequent MD simulations the H-bond to Ser113 was formed in 2–4 ns. The lifetime of the H-bond was in the range of 1–20 ns. On the trajectory with the most stable (20 ns) H-bond, the structural fluctuation of the α-carbon atoms of the receptor main chain was studied by calculating the root mean square fluctuations, the dynamic cross-correlation map, and the time-dependent dynamic cross-correlation. The analysis suggested a correlation transfer pathway Ser113 → Phe182 → (Leu259 or Tyr260) → Tyr291 induced by the ligand binding with a time scale of 4–6 ns.

Significance

The molecular mechanism of the functions of olfactory receptor (OR) proteins is largely unknown. One reason for this difficulty is the lack of experimental 3D structures. We used the state-of-the-art artificial intelligence tool AlphaFold2 to construct a 3D structural model of the mouse eugenol OR protein Olfr73 and performed molecular dynamics simulations of a membrane solution model. Structural fluctuations and their dynamic cross-correlations with and without the odorant binding were analyzed. The analysis proposed a correlation transfer pathway and its time scale from the ligand to a residue bound to the G protein.

Introduction

The mammalian olfactory systems are capable of detecting and discriminating between a wide range of chemically diverse molecules [1–9]. Olfactory receptors (ORs) are located on the cell membrane of olfactory nerves in the nasal cavity. Binding of odorant molecules to ORs leads to the dissociation of heterotrimeric G-proteins from the intracellular domain of the OR protein, which transduces chemical signals into neuronal electrical responses. The OR proteins belong to class A of G-protein coupled receptors (GPCRs) [10], which are characterized by seven transmembrane (TM) helices. Due to the current lack of experimental 3D structural data of OR proteins, molecular modeling studies of the OR mechanisms are limited [11–19]. Recent studies have mainly used homology modeling techniques with the available 3D structures of other class A proteins such as bovine rhodopsin and human β₂-adrenergic receptor as templates to construct 3D structural models of OR proteins.

The mouse eugenol olfactory receptor Olfr73 (mOR-EG) has been studied as a prototype of mammalian ORs. By combining computational molecular modeling with experimental techniques such as site-directed mutagenesis, eleven critical amino acid residues (Cys106, Ser113, Phe182, Phe203, Asn207, Glu208, Leu212, Phe252, Ile256, Leu259, Tyr260, see Figure 1) located in TM3, TM5 and TM6 and the second extracellular loop were identified as essential for the ligand binding [16,17]. Molecular dynamics (MD) simulations of the apo-form OR protein suggested an interesting mechanism for the transfer of structural change information to a residue (Tyr291) linked to the G-protein [18]. Another MD simulation study with a fingerprint interaction analysis showed that the binding pocket of Olfr73 is smaller but more flexible than those of non-olfactory GPCRs [19].

Figure 1 

A snapshot of the MD simulation showing the positions of the twelve key residues, Cys106^3.33, Ser113^3.40, Phe182, Phe203^5.42, Asn207^5.46, Glu208^5.47, Leu212^5.51, Phe252^6.47, Ile256^6.51, Leu259^6.54, Tyr260^6.55, and Tyr291^7.53 (the superscripts denote Ballesteros-Weinstein numbers [20,21]), relative to the odorant ligand eugenol (EG).

In this work, we further investigate the dynamic response of Olfr73 to ligand binding. The artificial intelligence tool AlphaFold2 [22–24] is used to construct a 3D structural model of Olfr73. The dynamics of both the apo and holo forms are studied by calculating the root mean square fluctuations (RMSF), the dynamic cross-correlation map (DCCM) [25–29], and the time-dependent dynamic cross-correlation (TDDCC) from MD trajectories. The dynamic structural response induced by the ligand binding suggests pathways of correlation transfer and their corresponding time scales.

Methods

Structure Model from AlphaFold2

We used the ColabFold web server [30] to construct a 3D structural model of Olfr73. (The access date was July 8, 2022.) The default options of ColabFold (AlphaFold2 with MMseqs2) were used. The amino-acid sequence was taken from UniProt entry Q920P2.

MD Simulation of Apo Form

The structure model of Olfr73 from AlphaFold2 is embedded in a palmitoyl-oleolyl-glycelo-phosphatidyl-ethanolamine (POPE) membrane in a saline solution, using the membgene3 package of the myPresto5 software suite [31]. The number of Na⁺ and Cl^– ions was adjusted to neutralize the system. The resulting numbers of molecules and ions in the simulation box were 229 POPE, 30 Na⁺, 41 Cl^–, and 10815 H₂O. After structural optimization by minimizing the total energy with positional constraints on the main chain atoms by harmonic potentials, constant pressure and temperature (NPT) simulations were performed to heat and equilibrate the system at 310 K and 1 atm. The simulation time step was 0.5 fs. The Berendsen method [32] with a relaxation time of 100 fs was used to control the pressure and temperature. Equilibration of the box lengths and the total potential energy was achieved in a 1.1 ns simulation. The box lengths were then averaged over an additional 100 ps simulation. The resulting box size was 105.9×78.9×77.7 ų (Figure 2). Next, a 20 ns NVT (constant volume and temperature) MD simulation was performed. The Gaussian constraint method [33,34] was used to control the temperature. The simulation time step was 2 fs with the SHAKE method [35,36] to constrain the bond lengths and angles involving hydrogen atoms. The zero-dipole summation method [37,38] was used to calculate the long-range electrostatic interaction. We used the software psygene-G [38] for the MD simulation. The Amber99 [39] and Lipid14 [40] potential force-fields were used.

Figure 2 

A snapshot of the MD simulation box.

Structure Sampling and Odorant Docking

The first 1 ns trajectory of the apo form was divided into twenty windows of 50 ps each. One configuration per each window of the Olfr73 protein was randomly sampled. An eugenol molecule was docked to target points of the sampled configurations near Ser113 and Asn207, covering the spatial region surrounded by TM3, TM5 and TM6 helices. The software sievgene [41] was used for the docking simulation. The atomic charge parameters for eugenol were determined by the electrostatic potential fitting with the RHF/6-31G(d) wave function. the software GAMESS [42] was used for the electronic structure calculations.

MD Simulation of Holo Form

The sampled and docked configurations of Olfr73 and eugenol were returned to the simulation box of the membrane solution. The atomic coordinates of the membrane and the saline solution were optimized by energy minimization with the harmonic constraints to the positions of Olfr73 and eugenol. After heating the temperature to 310 K in 100 ps, the constant NVT MD simulation was performed for 10 ns. In two out of twenty trajectories, stable H-bonding between Ser113 and eugenol was observed for 2–3 ns. The 3 ns H-bond trajectory was divided into ten 300 ps windows, and one configuration per each window was randomly sampled. For each configuration, the atomic velocities were randomized in the Maxwell-Boltzmann distribution to restart the constant NVT MD simulation. Six trajectories maintained the H-bond between eugenol and Ser113 for 1–20 ns. From the 20 ns H-bond trajectory (Supplementary Figure S3), the root mean square fluctuation (RMSF), the dynamic cross-correlation map (DCCM) [25–29], and the time-dependent dynamic cross-correlation (TDDCC) were calculated.

Results and Discussion

Structure from AlphaFold2

The resulting 3D structure from AlphaFold2 contained a disulfide bond between Cys98 and Cys180. The structure was compared with the more recent UniProt entry (AF-Q920P2-F1-model_v4.pdb) in the Supplementary Figure S1. The root mean square deviation of the aligned part (calculated with the PyMOL software [43]) was 0.240 Å.

Sampling of Hydrogen-Bonded Trajectories

Two representative trajectories of the distance between the O_H atom (the oxygen atom of the OH group) of eugenol and the O_γ atom of Ser113, the O_γ atom of Asp207, and the O_η atom of Tyr260 starting from the docking configurations are shown in Figure 3. Since the O-O distance shorter than 3.5 Å can be considered to form a H-bond, both panels in Figure 3 indicate that a stable H-bond is formed between eugenol and Ser113 for 2–3 ns. For Asp207, Glu208, and Tyr260, stable H-bonds with eugenol were not observed in Figure 3 or in the other eighteen sampled trajectories shown in the Supplementary Figure S2.

Figure 3 

Representative trajectories of the distance between the O_H atom of eugenol and the O atoms of Ser113, Asn207, and Tyr260, starting from the docking configurations.

Figure 4 shows a snapshot from the MD trajectory of the H-bond structure between eugenol and Ser113, indicating that eugenol is also H-bonded to Glu112. The double H-bond would be the main reason for the long lifetime (>2 ns) of the ligand binding.

Figure 4 

A snapshot of the hydrogen-bond structure between eugenol (EG) and Ser113^3.40 and Glu112^3.39.

Root-Mean-Squares Fluctuation

Next, we examine the fluctuation of the protein structure induced by the ligand binding. The root mean square fluctuation (RMSF) of the Cartesian coordinate r_i of the i-th atom is defined by

  

$${\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{F}}_i=\sqrt{⟨\mathrm{\Delta }{\mathit{r}}_i^2⟩}$$(1)

where Δr_i=r_i–⟨r_i⟩ is the deviation from the average ⟨r_i⟩. In this work, the fluctuations of the amino-acid residues are represented by those of the α-carbon atoms.

The calculated RMSFs in the apo and holo forms and their difference (holo minus apo) are shown in Figure 5. In both apo and holo forms, the fluctuation is comparatively smaller in the central regions of the TM helices. The RMSF is suppressed in the holo form for small residue numbers (1–12), in the beginning of TM3, after TM4, and in the last region of TM6. On the other hand, the RMSF is enhanced in the holo-form in the regions between TM1 and TM2, between TM3 and TM4, before TM5, between TM5 and TM6, and after TM7. Twelve key residues are taken from the left panel of Figure 5 and displayed in the right panel. Eleven of them (Cys106, Ser113, Phe182, Phe203, Asn207, Glu208, Leu212, Phe252, Ile256, Leu259, and Tyr260) are those reported in ref [17] as essential for the olfactory function, and one (Tyr291) is bonded to the G-protein [18]. A snapshot of their positions is shown in Figure 1. It can be seen in Figure 5 that the suppression of RMSF in the holo-form is remarkable for Ile256, Leu259, and Tyr260, which are located near the end or TM6.

Figure 5 

Root-mean-square fluctuation (RMSF) of the C_α atoms of Olfr73 in the apo and holo forms and their difference (holo minus apo). The left panel shows all residues. TM1–TM7 denote the transmembrane α-helices. The vertical lines below the bar plot indicate the twelve key residues shown in the right panel.

Dynamic Cross Correlation Map

The dynamic cross-correlation map (DCCM) [25–29] is a map of the cross-correlation or covariance between the fluctuations of the atomic Cartesian coordinates r_i and r_j, defined by

  

$$C_{ij}=\frac{⟨\Delta {\mathit{r}}_i\cdot \Delta {\mathit{r}}_j⟩}{\sqrt{⟨\Delta {\mathit{r}}_i^2⟩}\sqrt{⟨\Delta {\mathit{r}}_j^2⟩}}.$$(2)

The numerator is the average of the inner product between the coordinate shift vectors. The denominator is the normalization factor, which limits the value between –1 and +1. The DCCM is close to +1 when the atomic coordinate fluctuations are in the same direction, close to –1 when they are in the opposite direction, and close to 0 when there is no directional correlation. Analysis of the DCCM would reveal the correlated dynamics of structural fluctuations, and in particular, provide the key to identifying the correlation transfers from Ser113 H-bonded to the ligand to Tyr291 linked to the G-protein. (See the lines below Eq. (3) for the term “correlation transfer”).

The calculated DCCM for the apo and holo forms and their difference (holo minus apo) are shown in Figure 6. Their cross sections at Ser113 are shown in the left panel of Figure 7. (The corresponding plots at the other eleven key residues are displayed in Supplementary Figure S4.) The difference of the correlations is remarkable; between TM2 and TM3 (negative), in the last part of TM3 (positive), between TM4 and TM5 (negative), in the last part of TM5 (positive), between TM6 and TM7 (negative). This means that negative or positive correlations in the structural fluctuation are induced by the H-bond of the ligand eugenol to Ser113 and Glu112.

Figure 6 

Dynamic cross-correlation map (DCCM) of the C_α atoms of Olfr73 in the apo and holo forms and their difference (holo minus apo).

Figure 7 

Sections of the DCCMs in Figure 6 at Ser113 and difference of the DCCMs (holo minus apo) for the key residues.

The difference map for the 12 key residues is shown in the right panel of Figure 7. The map suggests several possible pathways of the correlation transfer. In particular, we note the correlations between Ser113 and Phe182 (negative), Phe182 and Leu259 (positive), and Leu259 and Tyr291 (negative). Leu259 and Tyr260 show similar behavior. Therefore, one of the simplest correlation transfer pathways would be Ser113 → Phe182 → (Leu259 or Tyr260) → Tyr291.

Time-Dependent Dynamic Cross Correlation

For the correlation transfer pathway identified in the previous section, we calculated the time-dependent dynamic cross-correlation (TDDCC) defined by

  

$$C_{ij}\left(t\right)=\frac{⟨\Delta {\mathit{r}}_i\left(0\right)\cdot \Delta {\mathit{r}}_j\left(t\right)⟩}{\sqrt{⟨\Delta {\mathit{r}}_i^2⟩}\sqrt{⟨\Delta {\mathit{r}}_j^2⟩}}$$(3)

This is the cross-correlation with the time delay t of the fluctuations of the atomic coordinates r_i and r_j. The DCCM of Eq. (2) is the simultaneous (t=0) part of the TDDCC. Thus, the TDDCC is a time-dependent extension of the DCCM. (To our knowledge, it has never been defined and computed in the literature). The amplitude of the DCCM measures the correlation of the structural fluctuations, and its sign indicates the parallel or antiparallel correlation. Therefore, the time evolution of the TDDCC would be informative on the time scale of the correlation transfer. (More precisely, in this paper we call the behavior of the TDDCC the “correlation transfer”). Since we compute the TDDCC from the equilibrium simulation, we assume that the essential dynamics of the non-equilibrium response induced by the ligand binding is embedded in the time correlation in the equilibrium (fluctuation-dissipation theorem).

The calculated TDDCC along the path Ser113 → Phe182 → (Leu259 or Tyr260) → Tyr291 is shown in Figure 8. The results with Leu259 and Tyr260 are almost identical. All TDDCCs change their qualitative behavior after the ligand binding. The TDDCC between Ser113 and Phe182 is overall a decaying function of time in the apo form, while it is basically an increasing function of time in the holo form. Interestingly, they have minimum or maximum at t≃6 ns. The TDDCC between Phe182 and Leu259 (or Tyr260) decays from positive to negative in 10 ns in the apo form, while in the holo form it oscillates near zero on an ns timescale. The TDDCC between Leu259 (or Tyr260) and Tyr291 increases from negative to positive in ~4 ns and stays near zero (or slowly decays to zero) in the apo form, while it decreases from positive to negative overall with oscillation in ns timescale.

Figure 8 

Time-dependent dynamic cross-correlation (TDDCC) of the apo and holo forms of Olfr73 between the key residues identified as being in the major correlation transfer pathway.

Conclusion

Structural fluctuations and dynamic cross-correlations in the mouse eugenol olfactory receptor Olfr73 were studied by molecular dynamics (MD) simulation. Among the candidate residues for H-bonding with the odorant ligand eugenol, Ser113, Asn207, and Tyr260, we find that only Ser113 forms a stable H-bond in the sampled trajectories. The lifetime of the H-bond was in the range of 1–20 ns. The structural fluctuations of the C_α atoms of the receptor main chain were investigated by calculating the root mean square fluctuations (RMSF), the dynamic cross-correlation map (DCCM), and the time-dependent dynamic cross-correlation (TDDCC). The analysis suggested a correlation transfer pathway Ser113 → Phe182 → (Leu259 or Tyr260) → Tyr291 induced by the ligand binding. The TDDCC indicated that the time scale of the correlation transfer was 4–6 ns. However, to our knowledge, there is no experimental evidence to verify these results. Time-dependent spectroscopies such as fluorescence resonance excitation transfer (FRET) [44], resonance Raman [45], and transient-grating [46] could be useful to clarify the time scale of the structural correlation transfer.

In this work, we focused on the twelve key residues identified in the previous experimental studies. The correlation transfer pathways highlighted in this work appear to be the most straightforward ones that directly involve the key residues. More extensive analysis would be required to clarify the role of other residues that may be directly or indirectly involved in the correlation transfer. In addition, the precise mechanism by which the dissociation of the G-protein is induced is still unclear. Work on these issues is ongoing.

Conflict of Interest

The authors declare no conflict of interest.

Author Contributions

CO set up and carried out the MD simulations, analyzed the data, and wrote the initial draft. KA designed the research, carried out additional calculations and analysis, and wrote the manuscript.

Data Availability

The evidence data generated in this study are available from the corresponding author on reasonable request.

A preliminary version of this work was deposited in arXiv (2311.02654) on 5th Nov. 2023, and bioRxiv (https://doi.org/10.1101/2023.11.05.565671) on 6th Nov. 2023.

Acknowledgments

We thank Dr. S. Katada and Prof. K. Touhara for useful discussions. This work was supported by JSPS KAKENHI Grant Number JP19K22173.

References

• [1]   Buck,  L.,  Axel,  R. A novel multigene family may encode odorant receptors: A molecular basis for odor recognition. Cell 65, 175–187 (1991). https://doi.org/10.1016/0092-8674(91)90418-x
• [2]   Kraft,  P.,  Bajgrowicz,  J. A.,  Denis,  C.,  Fráter,  G. Odds and trends: Recent developments in the chemistry of odorants. Angew. Chem. Int. Ed. Engl. 39, 2980–3010 (2000). https://doi.org/10.1002/1521-3773(20000901)39:17<2980::AID-ANIE2980>3.0.CO;2-%23
• [3]   Firestein,  S. How the olfactory system makes sense of scents. Nature 413, 211–218 (2001). https://doi.org/10.1038/35093026
• [4]   Zozulya,  S.,  Echeverri,  F.,  Nguyen,  T. The human olfactory receptor repertoire. Genome Biol. 2, (2001). https://doi.org/10.1186/gb-2001-2-6-research0018
• [5]   Zhang,  X.,  Firestein,  S. The olfactory receptor gene superfamily of the mouse. Nat. Neurosci. 5, 124–133 (2002). https://doi.org/10.1038/nn800
• [6]   Serizawa,  S.,  Miyamichi,  K.,  Nakatani,  H.,  Suzuki,  M.,  Saito,  M.,  Yoshihara,  Y., et al. Negative feedback fegulation ensures the one receptor-one olfactory neuron rule in mouse. Science 302, 2088–2094 (2003). https://doi.org/10.1126/science.1089122
• [7]   Malnic,  B.,  Godfrey,  P. A.,  Buck,  L. B. The human olfactory receptor gene family. Proc. Natl. Acad. Sci. U.S.A. 101, 2584–2589 (2004). https://doi.org/10.1073/pnas.0307882100
• [8]   Reed,  R. R. After the holy grail. Cell 116, 329–336 (2004). https://doi.org/10.1016/s0092-8674(04)00047-9
• [9]  Sell, C. S. Chemistry and the Sense of Smell (Wiley, New Jersey, 2014).
• [10]   Yang,  D.,  Zhou,  Q.,  Labroska,  V.,  Qin,  S.,  Darbalaei,  S.,  Wu,  Y., et al. G protein-coupled receptors: structure- and function-based drug discovery. Signal Transduct. Target. Ther. 6, 7 (2021). https://doi.org/10.1038/s41392-020-00435-w
• [11]   DiPizio,  A.,  Niv,  M. Y. Computational studies of smell and taste receptors. Isr. J. Chem. 54, 1205–1218 (2014). https://doi.org/10.1002/ijch.201400027
• [12]   Doszczak,  L.,  Kraft,  P.,  Weber,  H.-P.,  Bertermann,  R.,  Triller,  A.,  Hatt,  H., et al. Prediction of perception: Probing the hOR17-4 olfactory receptor model with silicon analogues of bourgeonal and lilial. Angew. Chem. Int. Ed. Engl. 46, 3367–3371 (2007). https://doi.org/10.1002/anie.200605002
• [13]   Gelis,  L.,  Wolf,  S.,  Hatt,  H.,  Neuhaus,  E. M.,  Gerwert,  K. Prediction of a ligand-binding niche within a human olfactory receptor by combining site-directed mutagenesis with dynamic homology modeling. Angew. Chem. Int. Ed. Engl. 51, 1274–1278 (2011). https://doi.org/10.1002/anie.201103980
• [14]   Floriano,  W. B.,  Vaidehi,  N.,  Goddard,  W. A.,  Singer,  M. S.,  Shepherd,  G. M. Molecular mechanisms underlying differential odor responses of a mouse olfactory receptor. Proc. Natl. Acad. Sci. U.S.A. 97, 10712–10716 (2000). https://doi.org/10.1073/pnas.97.20.10712
• [15]   Floriano,  W. B. Making sense of olfaction through predictions of the 3-D structure and function of olfactory receptors. Chem. Senses 29, 269–290 (2004). https://doi.org/10.1093/chemse/bjh030
• [16]   Katada,  S.,  Hirokawa,  T.,  Oka,  Y.,  Suwa,  M.,  Touhara,  K. Structural basis for a broad but selective ligand spectrum of a mouse olfactory receptor: Mapping the odorant-binding site. J. Neurosci. 25, 1806–1815 (2005). https://doi.org/10.1523/jneurosci.4723-04.2005
• [17]   Baud,  O.,  Etter,  S.,  Spreafico,  M.,  Bordoli,  L.,  Schwede,  T.,  Vogel,  H., et al. The mouse eugenol odorant receptor: Structural and functional plasticity of a broadly tuned odorant binding pocket. Biochemistry 50, 843–853 (2010). https://doi.org/10.1021/bi1017396
• [18]   Baud,  O.,  Yuan,  S.,  Veya,  L.,  Filipek,  S.,  Vogel,  H.,  Pick,  H. Exchanging ligand-binding specificity between a pair of mouse olfactory receptor paralogs reveals odorant recognition principles. Sci. Rep. 5, 14948 (2015). https://doi.org/10.1038/srep14948
• [19]   Yuan,  S.,  Dahoun,  T.,  Brugarolas,  M.,  Pick,  H.,  Filipek,  S.,  Vogel,  H. Computational modeling of the olfactory receptor Olfr73 suggests a molecular basis for low potency of olfactory receptor-activating compounds. Commun. Biol. 2, 141 (2019). https://doi.org/10.1038/s42003-019-0384-8
• [20]   Pándy-Szekeres,  G.,  Caroli,  J.,  Mamyrbekov,  A.,  Kermani,  A. A.,  Keserü,  G. M.,  Kooistra,  A. J., et al. GPCRdb in 2023: state-specific structure models using AlphaFold2 and new ligand resources. Nucleic Acids Res. 51, D395–D402 (2022). https://doi.org/10.1093/nar/gkac1013
• [21]   Ballesteros,  J. A.,  Weinstein,  H. Integrated methods for the construction of three-dimensional models and computational probing of structure-function relations in G protein-coupled receptors. Methods Neurosci. 25, 366–428 (1995). https://doi.org/10.1016/S1043-9471(05)80049-7
• [22]   Jumper,  J.,  Evans,  R.,  Pritzel,  A.,  Green,  T.,  Figurnov,  M.,  Ronneberger,  O., et al. Highly accurate protein structure prediction with AlphaFold. Nature 596, 583–589 (2021). https://doi.org/10.1038/s41586-021-03819-2
• [23]   Akdel,  M.,  Pires,  D. E. V.,  Pardo,  E. P.,  Jänes,  J.,  Zalevsky,  A. O.,  Mészáros,  B., et al. A structural biology community assessment of AlphaFold2 applications. Nat. Struct. Mol. Biol. 29, 1056–1067 (2022). https://doi.org/10.1038/s41594-022-00849-w
• [24]   Yang,  Z.,  Zeng,  X.,  Zhao,  Y.,  Chen,  R. AlphaFold2 and its applications in the fields of biology and medicine. Signal Transduct. Target. Ther. 8, 115 (2023). https://doi.org/10.1038/s41392-023-01381-z
• [25]   McCammon,  J. A. Protein dynamics. Rep. Prog. Phys. 47, 1 (1984). https://doi.org/10.1088/0034-4885/47/1/001
• [26]   da Penha Valente,  R. P.,  de Souza,  R. C.,  de Medeiros Muniz,  G.,  Ferreira,  J. E. V.,  de Miranda,  R. M.,  e Lima,  A. H. L., et al. Using accelerated molecular dynamics simulation to elucidate the effects of the T198F mutation on the molecular flexibility of the West Nile virus envelope protein. Sci. Rep. 10, 9625 (2020). https://doi.org/10.1038/s41598-020-66344-8
• [27]   Parida,  P. K.,  Paul,  D.,  Chakravorty,  D. The natural way forward: Molecular dynamics simulation analysis of phytochemicals from Indian medicinal plants as potential inhibitors of SARS-CoV-2 targets. Phytother. Res. 34, 3420–3433 (2020). https://doi.org/10.1002/ptr.6868
• [28]   Wang,  Z.,  Hu,  B.,  An,  Y.,  Wang,  J. Exploring the resistance mechanisms of distal D835V mutation in FLT3 to inhibitors. Oxid. Med. Cell. Longev. 2022, 3720026 (2022). https://doi.org/10.1155/2022/3720026
• [29]   Dash,  R.,  Munni,  Y. A.,  Mitra,  S.,  Choi,  H. J.,  Jahan,  S. I.,  Chowdhury,  A., et al. Dynamic insights into the effects of nonsynonymous polymorphisms (nsSNPs) on loss of TREM2 function. Sci. Rep. 12, 9378 (2022). https://doi.org/10.1038/s41598-022-13120-5
• [30]   Mirdita,  M.,  Schütze,  K.,  Moriwaki,  Y.,  Heo,  L.,  Ovchinnikov,  S.,  Steinegger,  M. ColabFold: making protein folding accessible to all. Nat. Methods 19, 679–682 (2022). https://doi.org/10.1038/s41592-022-01488-1
• [31]   Fukunishi,  Y.,  Mikami,  Y.,  Nakamura,  H. The filling potential method: A method for estimating the free energy surface for protein-ligand docking. J. Phys. Chem. B 107, 13201–13210 (2003). https://doi.org/10.1021/jp035478e
• [32]   Berendsen,  H. J. C.,  Postma,  J. P. M.,  van Gunsteren,  W. F.,  DiNola,  A.,  Haak,  J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81, 3684–3690 (1984). https://doi.org/10.1063/1.448118
• [33]   Evans,  D. J.,  Hoover,  W. G.,  Failor,  B. H.,  Moran,  B.,  Ladd,  A. J. C. Nonequilibrium molecular dynamics via Gauss’s principle of least constraint. Phys. Rev. A 28, 1016–1021 (1983). https://doi.org/10.1103/physreva.28.1016
• [34]   Nosé,  S. Constant temperature molecular dynamics methods. Prog. Theor. Phys. Suppl. 103, 1–46 (1991). https://doi.org/10.1143/ptps.103.1
• [35]   Ryckaert,  J.-P.,  Ciccotti,  G.,  Berendsen,  H. J. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23, 327–341 (1977). https://doi.org/10.1016/0021-9991(77)90098-5
• [36]   Kasahara,  K.,  Ma,  B.,  Goto,  K.,  Dasgupta,  B.,  Higo,  J.,  Fukuda,  I., et al. myPresto/omegagene: a GPU-accelerated molecular dynamics simulator tailored for enhanced conformational sampling methods with a non-Ewald electrostatic scheme. Biophys. Physicobiol. 13, 209–216 (2016). https://doi.org/10.2142/biophysico.13.0_209
• [37]   Fukuda,  I.,  Yonezawa,  Y.,  Nakamura,  H. Molecular dynamics scheme for precise estimation of electrostatic interaction via zero-dipole summation principle. J. Chem. Phys. 134, 164107 (2011). https://doi.org/10.1063/1.3582791
• [38]   Mashimo,  T.,  Fukunishi,  Y.,  Kamiya,  N.,  Takano,  Y.,  Fukuda,  I.,  Nakamura,  H. Molecular dynamics simulations accelerated by GPU for biological macromolecules with a non-Ewald scheme for electrostatic interactions. J. Chem. Theory Comput. 9, 5599–5609 (2013). https://doi.org/10.1021/ct400342e
• [39]   Wang,  J.,  Cieplak,  P.,  Kollman,  P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 21, 1049–1074 (2000). https://doi.org/10.1002/1096-987x(200009)21:12<1049::AID-JCC3>3.0.CO;2-F
• [40]   Dickson,  C. J.,  Madej,  B. D.,  Skjevik,  Å. A.,  Betz,  R. M.,  Teigen,  K.,  Gould,  I. R., et al. Lipid14: The Amber lipid force field. J. Chem. Theory Comput. 10, 865–879 (2014). https://doi.org/10.1021/ct4010307
• [41]   Fukunishi,  Y.,  Mikami,  Y.,  Nakamura,  H. Similarities among receptor pockets and among compounds: Analysis and application to in silico ligand screening. J. Mol. Graph. Model 24, 34–45 (2005). https://doi.org/10.1016/j.jmgm.2005.04.004
• [42]   Schmidt,  M. W.,  Baldridge,  K. K.,  Boatz,  J. A.,  Elbert,  S. T.,  Gordon,  M. S.,  Jensen,  J. H., et al. General atomic and molecular electronic structure system. J. Comput. Chem. 14, 1347–1363 (1993). https://doi.org/10.1002/jcc.540141112
• [43]  Schrödinger, L, DeLano W. (2020). PyMOL. http://www.pymol.org/pymol
• [44]   Nguyen,  T. A.,  Puhl,  H. L.,  Hines,  K.,  Liput,  D. J.,  Vogel,  S. S. Binary-FRET reveals transient excited-state structure associated with activity-dependent CaMKII - NR2B binding and adaptation. Nat. Commun. 13, 6335 (2022). https://doi.org/10.1038/s41467-022-33795-8
• [45]   Yamashita,  S.,  Mizuno,  M.,  Mizutani,  Y. High suitability of tryptophan residues as a spectroscopic thermometer for local temperature in proteins under nonequilibrium conditions. J. Chem. Phys. 156, 075101 (2022). https://doi.org/10.1063/5.0079797
• [46]   Terazima,  M. Spectrally silent protein reaction dynamics revealed by time-resolved thermodynamics and diffusion techniques. Acc. Chem. Res. 54, 2238–2248 (2021). https://doi.org/10.1021/acs.accounts.1c00113