Molecular Docking, MM-GBSA, and Molecular Dynamics Approach: 5-MeO-DMT Analogues as Potential Antidepressants

1. Introduction

Depression is now among the most severe mental disorders, as millions of people throughout the globe suffer from depression, and up to 40% of patients do not react properly to antidepressant medications ( 1 ). Depressive disorders, which are very prevalent, burdensome, and costly, are linked with significant impairments in role performance and quality of life, as well as medical comorbidity and mortality. Psychedelics are extremely potent psychoactive chemicals that change the mood as well as a variety of cognitive functions. Particularly, psychedelics are mentally safe, have little withdrawal symptoms, represent a low risk of addiction or dependency, and have no persistent physiological or psychological side effects ( 2 , 3 ). Furthermore, psychedelics have considerably less harmful side effects than traditional antidepressants, compared to frequently prescribed antidepressants, which may take several weeks to effect, deliver therapeutic results in patients more quickly, and may have beneficial, long-lasting benefits with only one pill or treatment session ( 4 , 5 ). Suicide is a potential consequence of severe depression. Approximately 700,000 individuals commit suicide each year. Suicide ranks as the fourth highest cause of mortality among persons in their prime aged 15-29. Low levels of serotonin, norepinephrine, and/or dopamine in the brain are the pathophysiological basis for depression, according to the monoamine hypothesis. Since increasing serotonin levels can assist enhance mood and alleviate depressive symptoms, the majority of research on depression has focused on this neurotransmitter ( 6 ). The decrease in central serotonin is associated with mood-tuned memory distortions, reward-related behavioural changes, and disruption of inhibitory emotional processing ( 7 ). The available data support the concept that decreased serotonin levels may produce clinical symptoms. Therefore, this may be remedied by elevating inter-synaptic serotonin concentrations. Except for the ionotropic 5-HT3R, all 5-Hydroxytryptamine receptors (5-HTRs) are members of the G protein-coupled receptor superfamily ( 8 ). The 5-HT1AR and 5-HT7R are two of the fourteen types of serotonin (5-HT) receptors that have been connected to anxiety, depression, and cognitive function through emotional learning and memory pathways ( 9 ). The 5-HT1AR is among the most well-studied 5-HTR subtypes, because of its role in anxiety-like behaviours, depression, and cognitive functions which are disrupted in numerous psychiatric diseases ( 10 - 15 ). In the majority of depressed patients, 5-HT1AR abnormalities have been reported ( 16 - 17 ). The 5-HT1AR agonists are 5-HT receptor subtype-specific and provide a potential pharmacologic strategy for the treatment of serious depression ( 8 ). The 5-HT1AR agonist exhibit strong hydrophobic contacts with residues Ile189, Trp358, Phe361, and Phe362, whereas lesser hydrophobic interactions with residues Ile113, Asp116, Val117, and Tyr390. The bound antagonist formed hydrophobic bonds with residues Val117, Ile189, Ser199, and Phe361 but not with Trp358 and Phe362 ( 18 ). The DMT (and by extension, 5-MeO-DMT) and psilocybin (and by extension, Psilocin) have structural similarities with serotonin, as do other serotonergic psychedelics. Traditional psychedelic (serotonergic) drugs influence the brain's 5-HTR. The 5-MeO-DMT was first synthesised in 1936 and isolated from the seeds of Dictyoloma incanescens DC in 1959, as well as Anadenanthera peregrina (L.) Speg ( 19 ). However, it is frequently linked to the endemic Sonoran Desert Toad (Bufo alvarius/Incilius alvarius), which is found in southwestern United States and northwest Mexico ( 20 ). The hallucinogen B. alvarius is present in the milky-white venom that the toad's parotid glands discharge ( 21 ). This toad poison includes 5-MeO-DMT, which is very beneficial as an entheogen in Central and South American spiritual, religious (shamanic), and social ceremonies, as well as bufotenine, a chemically related active metabolite ( 22 ). The enzyme responsible for o-demethylating 5-MeO-DMT to bufotenine is polymorphic cytochrome P450 2D6 (CYP2D6). When administered orally, it is rendered inactive by visceral monoamine oxidase deamination (MAO). Monoamine oxidase inhibitors (MAOIs), such as harmaline, may increase 5-MeO-psychoactive DMT's effects by increasing its interaction with 5-HT1A and 5-HT2A receptors ( 23 ). In a study conducted by Uthaug et al., 5-MeO-DMT's impact on anxiety and mood disorders was studied and found higher life satisfaction, convergent thinking, higher ratings of mindfulness, and lower ratings of depression, anxiety, and stress after administering the medication as a single inhalation of vapor from dried toad secretion ( 24 ).

Uthaug et al. found that inhaling synthetic 5-MeO-DMT was lowered stress biomarkers (such as cortisol), stress, and anxiety ratings ( 25 ). Studies have shown that 5-MeO-DMT has antidepressant and anxiolytic characteristics with no negative physical or psychological effects and limited potential for addiction ( 3 ).

The present study aimed to identify the binding modes of potent antidepressant analogues of 5-MeO-DMT through sequential docking protocols employing 70,000 different analogues. To discover the inherent chemical binding pattern in the 5-HT1AR active region, MM-GBSA and molecular dynamics simulations were conducted (7E2Y.pdb). This exploratory procedure is shown in the following flowchart (Figure 1).

Figure 1. Flow chart illustrating step-wise protocol for the identification of potential 5-MeO-DMT analogues as an antidepressant targeting 5-HT1AR (7E2Y.pdb)

2. Materials and Methods

2.1. Molecular docking

After being downloaded from Protein Data Bank, the protein preparation wizard of Schrodinger Suite 2020-1 was used to further construct the 3D X-ray crystal structure of the serotonin 1A receptor-Gi protein co-crystallized with serotonin (7E2Y.pdb) ( 26 ). In this 7E2Y.pdb, chain-R was used for this study which has the active site of the ligand (Figure 2), while the non-redundant chains were eliminated. By removing crystal fluids and adjusting bond ordering with hydrogen additions, the protein was prepared, and the Epik module was utilised to complete the protonation and tautomeric states of acidic and basic residues at pH 7.2 by including missing side chains and loops ( 27 - 29 ). To minimize the protein energy, we used the OPLS3e (Optimized Potentials for Liquid Simulations) molecular force-field ( 30 ), with the crystallographic heavy atoms relative standard deviation (RMSD) maintained at 0.30Å. Using a van der Waals scale of 1.0 for the receptor and a partial charge cut-off of 0.25, a grid box was created with its centre at the active site's centroid (x=103.03; y=114.79; z=108.36). In Marvin Sketch, 70.000 5MeO-DMT 3D conformers were created ( 31 ). The 5-HT1AR (7E2Y.pdb) binding site was utilized to commence a virtual workflow using the High Throughput Virtual Screening (HTVS), Standard Precision (SP), and Extra Precision (XP) modes with their default settings. To characterise the behaviour of tiny molecules at the binding site of target proteins and to comprehend fundamental biochemical processes, the molecular docking approach is used to simulate the atomic-level interaction between a small molecule and a protein. The initial steps in docking are the prediction of the ligand's conformation and pose (its projected placement and orientation inside these sites), followed by the assessment of the binding affinity ( 32 ). Using the glide gscore, glide energy data, and hydrogen bond analysis, the optimal binding pose was determined.

Figure 2. Structure of antidepressant target 5-HT1AR (7E2Y.pdb)

2.2. MM-GBSA free energy calculation studies

To assess the relative relevance of the enthalpy and entropy-associated components in the binding of the ligand-protein complex, the prime MM-GBSA technique was used to the generalized-born/surface area (GB/SA) continuum solvent model ( 33 ). The formula was used to determine the energy contributions (in kilocalories per mole) from molecular mechanics, polar solvation, and a non-polar solvation factor.

ΔGbind = Gcomplex - Gprotein - GligandΔGbind = calculated binding free energy of the complex

Gcomplex = binding free energy of the minimized complex

Gprotein = binding free energy of receptor

Gligand = binding free energy of unbound ligand

2.3. Molecular dynamics (MD) simulation

Utilizing the Desmond module of Schrodinger 2020-1, LLC, New York, NY, we performed a molecular dynamics simulation to examine the binding behaviour of a highly rated molecule at the atomic level and to comprehend the molecular interaction analysis ( 34 ). The complex of 20453/7E2Y was solvated using the TIP3P water model in an orthorhombic periodic boundary condition with a 20Å buffer zone between the protein atoms and box borders. Adding 0.15M NaCl counter ions helped to neutralise the solvated system. Following that, the system’s energy was minimized using the OPLS3e force field ( 30 ), keeping the settings default. With the help of the smooth particle mesh Ewald technique, the long-range electrostatic interactions were computed with a 1e-09 tolerance. At a cut-off radius of 0.9, the van der Waals and coulomb interactions in the short-range were estimated. At 300K and 1 bar of pressure, a simulation lasting 100ns total was run under an isothermal/isobaric ensemble (NPT). At 100 and 200 ps, respectively, the Nose-Hoover chain thermostat and Martyna-Tobias-Klein barostat approaches were combined. Short-range non-bonded and long-range electrostatic forces accordingly, several time-step algorithms Reference System Propagator Algorithm (RESPA) were used at 2, 2, and 6 fs. At every 100ps intervals, the data was gathered, and the generated trajectories were examined.

3. Results and Discussion

3.1. Molecular docking and MM-GBSA free energy calculation

A total of 700,000 compounds were sketched using Marvin and evaluated using Schrodinger Suite's structure-based virtual screening workflow. Docking with HTVS was performed initially, followed by docking with SP and XP, and keeping default parameters to standard practice. After visually examining the bound postures and hydrogen bond creation of the highest-scoring hits, 22 compounds, including 5-MeO-DMT, were selected. With the aid of the MM-GBSA of the PRIME module, the amount of time spent in the post-docking phase was minimised. The glide scores and PRIME MM-GBSA energy values were tabulated in tables 1 and 2, respectively. Figure 3 depicts the 3D interaction of top hits (20453, 15538, 29999, 15945) and 5-MeO-DMT with 5-HT1AR (7E2Y.pdb).

Figure 3. 3D interaction diagram of top hits and standard in the active site of 5-HT1AR (7E2Y.pdb). (a) 20453 (b) 15538 (c) 29999 (d) 15945 (e) 5-MeO-DMT

Figures 4.1 and 4.2 depicts the 2D interactions of the top hits and 5-MeO-DMT. Hydrogen bonding, hydrophobic, and salt bridge interactions were observed in the XP docked poses of the selected hits, primarily with binding site residues Thr196, Thr121, Ser190, Ser199, Asp116, Phe362, Phe361, Pro369, Leu368, Ala365, Val117, Cys120, Tyr390, Tyr195, Ala203, Ile167, Ile189, Val364, Trp358, and Ile124 (Supplementary Table 1). The glide score (Table 1) varies from -11.41 to -6.53 kcal/mol. The hit molecule, 20453 with glide score -10.96kcal/mol (Figure 4.1. (a)) established three hydrogen bonds with the active site residues of the serotonin 1A receptor - Gi protein, nitrogen at the first position of the indole ring created a hydrogen bond with the backbone OH of Ser199 (-NH...OH-; 2.54Å). Nitrogen connected at the 7th position directly from indole produced a hydrogen bond with the side chain OH of Ser199 (-NH...OH-; 2.07Å), whereas the amide group attached to the 7th position of indole formed a hydrogen bond with the side chain OH of Ser190 (-NH...OH-; 1.97Å). Nitrogen connected to the third position of indole through the ethyl side chain produces a salt bridge interaction with Asp116. Phe361 and Phe362 interacted with benzene and the indole pyrrole ring through hydrophobic interactions. The complex was stabilized by 14 hydrophobic contacts between Pro369, Leu368, Ile189, Ala365, Val364, Tyr195, Phe362, Phe361, Trp358, Ile167, Ala203, Cys120, Val117, and Tyr390. Post docking revealed binding free energies (bind) ranging from -55.76 to -30.07 kcal/mol for the top-ranking postures of the chosen hits. As shown in Table 2, the Van der Waal energy terms (Δ Vdw) -21 to -43.33kcal/mol and hydrophobic energy terms (Δ lipo) -12.05 to -26.36 kcal/mol favour total binding energy.

Figure 4.1. 2D interaction diagram of top hits in the active site of 5-HT1AR (7E2Y.pdb)(a) 20453 (b) 20675 (c) 28334 (d) 30813

Figure 4.2. 2D interaction diagram of 5-MeO-DMT in the active site of 5-HT1AR (7E2Y.pdb)

3.2. Molecular dynamics simulation

A molecular dynamic simulation of the docked posture of the 20453/7E2Y complex that lasted 100ns revealed that the RMSD of protein Cα atoms (Figure 5a) were stabilised upon ligand binding, with slight oscillations ranging from 2.7 to 4.3Å. Up to 100ns (Figure 5a) demonstrates that the RMSD of the ligand remained rather stable (between 2.5 and 3.5Å). The protein's root mean square fluctuations (RMSF) (Figure 5c) seems to be steady, with the exception of amino acids 140 and 190, which exhibited larger variations up to 5.4 Å, which are seen in the loops. The compound was positioned in the active pocket by forming hydrogen bonding, hydrophobic bonding, ionic and water bridged interactions with Asp116, Val117, Cys120, Ile167, Ile189, Ser190, Tyr195, Thr196, Ser199, Trp358, Phe361, Phe362, Ala365, Met377, Leu381, Ile385, and Tyr390 from the ligand interaction fraction (Figure 5b). As the analogues were showing hydrophobic interaction with Trp358 and Phe362; therefore, this predicts that the activity is not antagonistic because it is known that the bound antagonist does not form hydrophobic bonds with residues Trp358 and Phe 361²². Hence, the analogues may show probable agonistic action after binding to the target site. The nitrogen present in the first position of indole exhibited one hydrogen bond with Ser199 at 94%, the amino group present in the third position of indole attached through ethyl side chain exhibited hydrogen bonding with Asp116 at 53%, and nitrogen attached directly at the seventh position of indole participates in forming hydrogen bond with Ser199 at 47%, according to the 2D ligand interaction diagram (Figure 5d) of the 100ns simulation trajectory. The nitrogen in the amide group linked to the 7th position of indole forms two hydrogen bonds with Ser190 and Thr196 at 99% and 94% of simulation trajectory, respectively. Asp116 established substantial water bridge contacts with the amino group present in the third position, which was connected through the ethyl side chain, and it also displayed ionic interactions with the same amino acid residue.

Depression is characterised by persistent emotions of sorrow and loss of interest. In addition to emotional and behavioural changes, it may also create cognitive and physical health problems. Clinical depression is also characterised as extreme sadness. Normal everyday tasks may seem difficult, and you may feel that life is not worth fighting for. Since depression is such a prevalent mental health issue, researchers have a strong motivation to create effective antidepressants for its estimated 5% worldwide prevalence. Dopamine, norepinephrine, and serotonin are the three major monoamine neurotransmitters in the brain and imbalances in these three neurotransmitters have been related to specific symptoms of clinical depression. The present study aimed to develop 5-Meo-DMT analogues that may be used as an antidepressant since they are a strong agonist for the 5-HT1AR. A High throughput Virtual Screening approach was run for a library of 70,000 5-MeO-DMT analogues, and sequential docking was done at three accuracy levels: HTVS, SP, and XP modes. Post-docking, the top 21 hit compounds were identified for further studies. According to the XP data, the Glide gscores for the compounds varied from -11.41 to -6.53 kcal/mol. When compared to the standard 5-MeO-DMT with binding affinity -7.75 kcal/mol, most compounds showed better binding affinity. Hydrophobic energy terms (Δ lipo) in PRIME MM-GBSA experiments varied from -12.05 to -26.80kcal/mol, while van der Waals energy (ΔVdw) ranged from -21 to -43.33kcal/mol. The 20453/7E2Y complex is stabilised in the catalytic pocket by hydrogen bonding, hydrophobic, ionic, and water bridging interactions, according to a 100ns molecular dynamic simulation. As expected by molecular dynamic simulations, strong hydrogen bonds were discovered between 20453/7E2Y and Ser199 (94% and 47%, respectively), Thr196 (94%), Ser190 (99%), and Asp116 (53%). The results demonstrates that further modifications in pharmacophoric features may help to improve agonistic activity against 5HT1AR. Computerized study revealed that the compounds 20453, 20675, 28334, 30813, 15538, 29740, 2235, 48166, 15539, 29999, 14836, 15612, 15945, and 15613 may have therapeutic potential against 5HT1AR and are likely to be beneficial after further refining. The 5-MeO-DMT is a psychedelic with hallucinogenic effects, and this study is an attempt to design analogues that might have better therapeutic activity and partial/complete cessation of psychedelic side effects. The in-silico structuring method used in the present work aided in the discovery of numerous lead compounds, and this may partly help to elucidate the molecular mechanism on future assessments, such as in-vitro and in-vivo tests.

Figure 5a. RMSD graph for the 100ns simulation trajectory of 20453/7E2Y complex

Figure 5(b). Ligand interaction fraction for the 100ns simulation trajectory of 20453/7E2Y Complex

Figure 5(c). RMSF graph for 100ns simulation trajectory of 20453/7E2Y complex

Figure 5d. 2D interaction diagram for the 100ns simulation trajectory of 20453/7E2Y complex

Authors' Contribution

Study concept and design: R.K. and K.R, Acquisition of data: R.K. and S.Jm Analysis and interpretation of data: N.M. and P.N., Drafting of the manuscript: R.K. and K.R., Critical revision of the manuscript for important intellectual content: R.K., K.R., and M.R.I, Statistical analysis: R.K, Administrative, technical, and material support: K.R

Ethics

Not applicable.

Conflict of Interest

The authors declare that they have no conflict of interest.

References

1. Fox ME, Lobo MK. The molecular and cellular mechanisms of depression: a focus on reward circuitry. Molecular psychiatry. 2019; 24(12):1798-815.
2. Davis AK, So S, Lancelotta R, Barsuglia JP, Griffiths RR. 5-methoxy-N, N-dimethyltryptamine (5-MeO-DMT) used in a naturalistic group setting is associated with unintended improvements in depression and anxiety. The American journal of drug and alcohol abuse. 2019; 45(2):161-9.
3. Gable RS. Risk assessment of ritual use of oral dimethyltryptamine (DMT) and harmala alkaloids. Addiction. 2007; 102(1):24-34.
4. Johnson MW, Griffiths RR. Potential therapeutic effects of psilocybin. Neurotherapeutics. 2017; 14:734-40.
5. Delgado PL. Depression: the case for a monoamine deficiency. Journal of clinical Psychiatry. 2000; 61(6):7-11.
6. Hasler G. Pathophysiology of depression: do we have any solid evidence of interest to clinicians?. World psychiatry. 2010; 9(3): 155.
7. Hoyer D, Hannon JP, Martin GR. Molecular, pharmacological and functional diversity of 5-HT receptors. Pharmacology Biochemistry and Behavior. 2002; 71(4):533-54.
8. Stiedl O, Pappa E, Konradsson-Geuken Å, Ögren SO. The role of the serotonin receptor subtypes 5-HT1A and 5-HT7 and its interaction in emotional learning and memory. Frontiers in pharmacology. 2015; 6:162.
9. Heisler LK, Chu H-M, Brennan TJ, Danao JA, Bajwa P, Parsons LH, et al. Elevated anxiety and antidepressant-like responses in serotonin 5-HT1A receptor mutant mice. Proceedings of the National Academy of Sciences. 1998; 95(25):15049-54.
10. Parks CL, Robinson PS, Sibille E, Shenk T, Toth M. Increased anxiety of mice lacking the serotonin1A receptor. Proceedings of the National Academy of Sciences. 1998; 95(18):10734-9.
11. Toth M. 5-HT1A receptor knockout mouse as a genetic model of anxiety. European journal of pharmacology. 2003; 463(1-3):177-84.
12. Lucki I. Behavioral studies of serotonin receptor agonists as antidepressant drugs. The Journal of clinical psychiatry. 1991; 52:24-31.
13. Ögren SO, Eriksson TM, Elvander-Tottie E, D’Addario C, Ekström JC, Svenningsson P, et al. The role of 5-HT1A receptors in learning and memory. Behavioural brain research. 2008; 195(1):54-77.
14. Millan MJ, Agid Y, Brüne M, Bullmore ET, Carter CS, Clayton NS, et al. Cognitive dysfunction in psychiatric disorders: characteristics, causes and the quest for improved therapy. Nature reviews Drug discovery. 2012; 11(2):141-68.
15. Hervás I, Artigas F. Effect of fluoxetine on extracellular 5-hydroxytryptamine in rat brain. Role of 5-HT autoreceptors. European journal of pharmacology. 1998; 358(1):9-18.
16. Artigas F. Developments in the field of antidepressants, where do we go now?. European Neuropsychopharmacology. 2015; 25(5): 657-70.
17. Yuan S, Peng Q, Palczewski K, Vogel H, Filipek S. Mechanistic studies on the stereoselectivity of the serotonin 5‐HT1A receptor. Angewandte Chemie. 2016; 128(30):8803-7.
18. Krebs-Thomson K, Ruiz EM, Masten V, Buell M, Geyer MA. The roles of 5-HT 1A and 5-HT 2 receptors in the effects of 5-MeO-DMT on locomotor activity and prepulse inhibition in rats. Psychopharmacology. 2006; 189:319-29.
19. Malaca S, Lo Faro AF, Tamborra A, Pichini S, Busardò FP, Huestis MA. Toxicology and analysis of psychoactive tryptamines. International Journal of Molecular Sciences. 2020; 21(23):9279.
20. Hoshino T, Shimodaira K. Über die synthese des bufotenin-methyl-äthers (5-methoxy-n-dimethyl-tryptamin) und bufotenins (synthesen in der indol-gruppe. xv). Bulletin of the Chemical Society of Japan. 1936; 11(3):221-4.
21. Sherwood AM, Claveau R, Lancelotta R, Kaylo KW, Lenoch K. Synthesis and characterization of 5-MeO-DMT succinate for clinical use. ACS omega. 2020; 5(49):32067-75.
22. Uthaug M, Lancelotta R, Van Oorsouw K, Kuypers K, Mason N, Rak J, et al. A single inhalation of vapor from dried toad secretion containing 5-methoxy-N, N-dimethyltryptamine (5-MeO-DMT) in a naturalistic setting is related to sustained enhancement of satisfaction with life, mindfulness-related capacities, and a decrement of psychopathological symptoms. Psychopharmacology. 2019; 236:2653-66.
23. Weil AT, Davis W. Bufo alvarius: a potent hallucinogen of animal origin. Journal of ethnopharmacology. 1994; 41(1-2):1-8.
24. Jiang X-L, Shen H-W, Mager DE, Yu A-M. Pharmacokinetic interactions between monoamine oxidase A inhibitor harmaline and 5-methoxy-N, N-dimethyltryptamine, and the impact of CYP2D6 status. Drug Metabolism and Disposition. 2013; 41(5):975-86.
25. Xu P, Huang S, Zhang H, Mao C, Zhou XE, Cheng X, et al. Structural insights into the lipid and ligand regulation of serotonin receptors. Nature. 2021; 592(7854):469-73.
26. Greenwood JR, Calkins D, Sullivan AP, Shelley JC. Towards the comprehensive, rapid, and accurate prediction of the favorable tautomeric states of drug-like molecules in aqueous solution. Journal of computer-aided molecular design. 2010; 24(6-7):591-604.
27. Shelley JC, Cholleti A, Frye LL, Greenwood JR, Timlin MR, Uchimaya M. Epik: a software program for pK a prediction and protonation state generation for drug-like molecules. Journal of computer-aided molecular design. 2007; 21:681-91.
28. Harder E, Damm W, Maple J, Wu C, Reboul M, Xiang JY, et al. OPLS3: a force field providing broad coverage of drug-like small molecules and proteins. Journal of chemical theory and computation. 2016; 12(1):281-96.
29. Cherinka B, Andrews BH, Sánchez-Gallego J, Brownstein J, Argudo-Fernández M, Blanton M, et al. Marvin: A tool kit for streamlined access and visualization of the SDSS-IV MaNGA data set. The Astronomical Journal. 2019; 158(2):74.
30. McConkey BJ, Sobolev V, Edelman M. The performance of current methods in ligand–protein docking. Current Science. 2002; 845-56.
31. Jacobson MP, Pincus DL, Rapp CS, Day TJ, Honig B, Shaw DE, et al. A hierarchical approach to all‐atom protein loop prediction. Proteins: Structure, Function, and Bioinformatics. 2004; 55(2):351-67.
32. Scalable algorithms for molecular dynamics simulations on commodity clusters. Proceedings of the 2006 ACM/IEEE Conference on Supercomputing; 2006.
33. Sidler D, Riniker S. Fast Nosé–Hoover thermostat: molecular dynamics in quasi-thermodynamic equilibrium. Physical Chemistry Chemical Physics. 2019; 21(11):6059-70.
34. Lippert RA, Predescu C, Ierardi DJ, Mackenzie KM, Eastwood MP, Dror RO, et al. Accurate and efficient integration for molecular dynamics simulations at constant temperature and pressure. The Journal of chemical physics. 2013; 139(16)