Introduction
As a brain-destroying disorder, Alzheimer disease (AD) is the most common form of dementia, which is widely seen in industrialized countries [1]. Besides its relation to the cholinergic system, AD is associated with a reduction in acetylcholine (ACh) levels in areas of the brain involved in emotional responses, behavior learning, and memory. AD can also be characterized by fibrous nerve nodes, amyloid-beta (Aβ) plaques, and degeneration of cholinergic and brainstem nerves [2]. The acetylcholine synaptic availability is reduced due to the loss of cholinergic cells in the forebrain, resulting in cognitive impairment in AD [3]. In this regard, one can symptomatically treat AD by increasing the synaptic ACh level by inhibiting the enzyme acetylcholinesterase (AChE) [4]. Thus, tacrine, donepezil, galantamine, and rivastigmine, as AChE inhibitors, can be regarded as the mainstay drugs in managing AD clinical symptoms [5, 6, 7]. Studies show that AChE can also be associated with specific non-catalytic actions, i.e., Aβ accumulation activity. AChE binds to Aβ and enhances the formation of amyloid fibril through a pool of amino acids adjacent to the Peripheral Anionic Site (PAS) of the enzyme. The molecules interacting with PAS or with both peripheral and catalytic AChE binding sites would prevent the AChE activity toward Aβ [8]. Moreover, it has also been reported that some AChE inhibitors aim at both facilitating cholinergic transmission and interfering with the deposition, synthesis, and accumulation of toxic Aβ. Therefore, the inhibition of AChE can be an important goal for effective AD management due to the increased synaptic availability of ACh and the reduction of Aβ deposition. Consequently, the dual-binding compounds (to catalytic and environmental sites) have novel therapeutic agents to effectively manage AD. Moreover, as AChE inhibitors, coumarin proposes the following hypotheses: 1) This heterocyclic compound, which is well-suited to AChE inhibitory activity, involves insaculin with cognitive functions. Besides, it can also act as a peripheral anionic site, interacting with the peripheral binding site. 2) In phenylpiperazine groups, regarding potent AChE inhibitors, the nitrogen atom acts as a positive charge center, and 3) attached to the piperazine ring, the phenyl ring acts as a choline-binding site [9].
Found in a wide variety of plant species, coumarins are naturally occurring phytochemicals with a broad range of biological and pharmacological activities, i.e., anti-tumor [10], anti-inflammatory [11], hepatoprotective [12], anti-HIV-1 [12], anti-allergic [12], antiviral [12], antimicrobial [12], antifungal [12], antiasthmatic [12], antidepressant [13], antidiabetic [13], antinociceptive [14], and antioxidant effects [15]. As indicated in the previous studies, coumarin analogs, which are chemically synthesized, show potent AChE inhibitory activities [16]. In addition, when the aromatic center of coumarin is functionalized, new analogs would be developed which can inhibit Aβ aggregation. There are also more reports on the regeneration of memory and anti-amnestic characteristic of coumarin derivatives in various models of amnesia [17]. Coumarin derivatives also protect neurons when opposed to Aβ-induced free radicals and oxidative stress [18]. 
Because of the biological and pharmacological properties of coumarin derivatives and the need to search novel candidates which would pioneer the design of AChE inhibitors, we aimed at investigating the inhibitory activity of a group of bis-coumarin against AChE inhibitors. These compounds were synthesized and characterized using the reaction between benzaldehyde and 4-hydroxycoumarin in high to excellent yields and in the presence of ZnO Nanoparticles (ZnO-NPs) as a highly effective heterogeneous catalyst in [bmim] BF4 (Figure 1).

Moreover, it is worth mentioning that synthetic chemists have paid great attention to metal oxides because they are non-toxic and stable in different conditions and have dual acid/base properties. Besides, because of the high surface-to-volume ratio, we can improve the performance of these heterogeneous catalysts using nano-sized metal oxides [14, 15, 16]. 
Recently, mineral oxides have proved helpful to chemists in the laboratory and industry due to the good activation of adsorbed compounds and reaction rate enhancement, selectivity, easier work-up, recyclability of the supports, and their eco-friendly, green, reaction conditions. Thus, the heterogeneous solid base catalysts have been recognized as potential alternatives to homogeneous organic essential catalysts [16, 17, 18, 19]. 
Thus, the heterogeneous solid base catalysts have been recognized as potential alternatives to homogeneous organic basic catalysts [19, 20, 21, 22, 23, 24]. Furthermore, to have a clear picture of the biological activity and the related molecular basis, we aimed at performing molecular docking studies and molecular dynamics (MD) simulation in a solvent system. In this sense, regarding the ChEs binding site, we assess the coumarin derivatives’ stability and their dynamic behavior. 


Materials and Methods
We purchased the needed chemicals from commercial sources and then applied them as received without further purification: Ellman’s reagent (DTNB), acetylcholinesterase, acetylthiocholine iodide, 5,50-dithiobis(2-nitrobenzoic acid), and 2-([1-benzylpiperidin-4-yl] methyl)-5,6-dimethoxy-2,3-dihydro-1H-inden-1-one (donepezil). Moreover, we used CyberScan pH 510 (Eutech Instruments, Singapore), microplate fluorescence spectrometer, single-channel electronic micropipette (Eppendorf AG, Hamburg, Germany), mQuant (BioTek Instruments, Winooski, Vermont, USA), and electronic multichannel micropipette. Besides, the melting points of the resulted compounds were determined using a Stuart scientific melting point apparatus in open capillary tubes. Furthermore, the 1H- and 13C-NMR spectra were recorded applying a 300 MHz spectrometer (Bruker DPX-300) using TMS. Regarding the column chromatography, silica gel was used. Also, the purity of the products and the progress of the reactions were examined applying a Thin-Layer Chromatography (TLC) (E. Merck Kieselgel 60 F254 layer thickness: 0.25 mm).

Synthesis of final compounds 
Regarding the significance of bis-coumarin derivatives, we aimed to find a protocol with a high yield. Given the recent advances in efficient synthetic and environmentally friendly methods, the possibility of preparing bis-coumarins was investigated using the condensation of aryl aldehydes and 4-hydroxycoumarin. Besides, we applied ZnO nanoparticles as an effective heterogeneous catalyst in [bmim]BF4. In this sense, we investigated the effect of different reaction parameters to reduce reaction period and waste materials and increase the yield. Accordingly, Table 1 presents the results of our investigation on the effect of reaction conditions, i.e., catalyst temperature, on the compound 1a synthesis using the condensation of 4-bromobenzaldehyde 1 (1 mmol) with 4-hydroxycoumarin 1 (2 mmol) as the model reaction.

Considering the obtained results, we applied different aromatic aldehydes to extend the scope of the reactions, as illustrated in Table 2.

Moreover, IR, 1H NMR, and 13C NMR were used to characterize the synthesized products, which were then compared to those mentioned in the literature.

Preparation of solutions and In vitro AChE inhibition assay
To determine the inhibitory activity, we used a substrate, i.e., acetylthiocholine iodide, to assay the tested compounds. Moreover, we applied Ellman’s method and performed the related studies in a 96-well microplate [25].

Molecular docking
As shown in Figure 2, during the docking simulations, we aimed to download the crystal structure of AChE (PDB ID:1Q83) [26], from the Protein Data Bank (PDB), as the receptor.

After removing all ligands, i.e., water molecules and ions, and adding the polar hydrogens, the partial atomic charge was calculated by adding the Kollman-united charge. The 3D structures of ligands were drawn in ChemDraw Ultra 8.0 and, then, the structures were saved as .pdb files after geometry optimization using the semi-empirical AM1 Hamiltonian, as indicated in Table 1. To prepare the ligands, we assigned Gasteiger-Marsili charges and rotatable bonds to them. The protein and ligand structures were converted to pdbqt format. Moreover, we performed the docking calculations using the docking program, Autodock Vina [27], which are based on a Genetic Algorithm (GA) [28]. Exhaustiveness was set to 100 in the docking procedure for each protein-ligand complex, and then the binding affinities of protein targets were determined (Table 3).


PLIP webserver (Technical University of Dresden) was also applied to investigate the formed complexes [29]. As shown in Figure 3, reproducing the binding pose of the co-crystallized ligand as a method for validating the docking protocols was used.

In this sense, it implied the aptness of the used setup for our docking study, which would reproduce the crystal binding model, confirmed by the law energy score=-14.7 kcal/mol between the added ligand and the docked pose. 

MD simulations protocol
To efficiently design the new inhibitors, specific information regarding the interaction properties of inhibitors is required. Two complexes with the best IC50 (L14-ACHE and L15-ACHE) values were subjected to molecular dynamics simulation to understand the interacting behavior and dynamics of these compounds. The selected complexes from the docking study were subjected to 100 ns MD simulations using GROMACS 2020 [27]. The ligand topologies were obtained from the CGenFF server. As an input file for MD simulations, the best ligand-protein complexes with the lowest docking score from the docking study and favorable interactions were selected. The CHARMM36 (C36) force field model was used for all MD simulations. After solvation in the box, to make the whole system neutral, 7 Na+ ions were introduced to the system. Each system was energy-minimized and equilibrated to achieve the appropriate volume. Overall, 10000 frames were obtained from each production simulation. The results of root mean square deviation RMSD, RMSF, and Rg were evaluated to confirm the protein stability for the two systems.

MM/PB SA binding-free-energy calculations
A stable short region from the simulated trajectory of the docked complexes was extracted to estimate the binding energy according to the Molecular Mechanics-Poisson Boltzmann Surface Area (MM-PBSA) method applying g_mmpbsa. This method is commonly applied in biomolecular interaction studies, and scoring function in computational drug design is expressed as Equation 1:

1. ΔGbinding=Gcomplex−(Gprotein+Gligand)

Where Gprotein and Gligand can be regarded as the free energies of the protein and ligand in solvent and Gcomplex depicts the total free energy of the protein-ligand complex, respectively. This approach consists of three energetic terms Equation 2:

2. Gx=⟨EMM⟩+⟨Gsolvation⟩−TS

The first onesignifies the changing potential energy in the vacuum, including bonded terms, i.e., torsion, bond and angle energies, and nonbonded terms, i.e., electrostatic and van der Waals interactions Equation 3:

3. EMM=Ebonded+Enonbonded=Ebonded+(EvdW+Eelec) 

The other term, (Gsolvation), is quantified by the polar and nonpolar solvation energies to calculate the desolvation of the different species in the gas phase (Equation 4).

4. Gsolvation=Gpolar+Gnonpolar

The third term (TS) would signify the entropic contribution to the free energy where S and T are the entropy and the temperature, respectively. This analysis estimates the energies combined with the protein-ligand binding during the MD simulation. However, because of the relatively high computational demand to calculate the entropy, the conformational entropy was not considered.


Results 
Chemistry and biological activity
We aimed at synthesizing, purifying, and characterizing the target compounds. TLC was used, and the melting point was used to assess the purity. Besides, we evaluated the compounds as positive controls in vitro for AChE inhibitory activities using the modified element method [25]. Furthermore, considering the use of donepezil, we could calculate the inhibition ratio against AChE, based on the product’s response in the mentioned reaction. Table 3 presents the results of the inhibitory activity of the mentioned coumarin derivatives. The inhibitory activity of the following compounds, L2, L5, L6, L7, L11, and L3, was moderate with IC50 values of 0.100 and 0.900 μM. Nevertheless, the inhibition was affected by adding F, Cl, NO2 groups on position 4, compared to the compounds L1, L2, L5, L6, L7, L11, and L3. The synthesized compounds, L1, L4 L8, L9, and L10, exhibited good inhibitory activity with the following IC50 values 0.800-0.770 μM. Moreover, L14 and L15 compounds showed considerable inhibitory activities with IC50 values of 0.20 and 0.50 μM, which can be due to the presence of the amines groups. Thus, compound L14 is the highest with an IC50 of 0.20 μM. 

Molecular docking
We aimed at performing the molecular docking simulation over the synthesized compounds on the AChE protein. In this sense, Table 3 indicates the sorted results from the lowest to the highest docking score. The docking score of these drugs (between -10 kcal/mol and 12.6 kcal/mol) would indicate their effective interaction with the receptor by the hydrogen-bonding interactions with the binding pocket residues. The 3D diagram for the docking pose of ligands in the AChE binding site can be seen in Figure 4.

Previous studies propose the Trp286 as the key residue in the binding pocket, which can be related to the adhesion function of AChE. Two complexes of L14-AChE and L15-AChE were important for us mainly because they have the best-measured IC50 . The docking conformation of L14 shows three hydrogen bonds with TYR124, and seven hydrophobic interactions with TYR72, LEU76, TRP286, PHE297, TYR337, and PHE338 and π-stacking interactions with TYR72, TYR124, TRP286, and TYR341 (Figure 2b and 2c). This ligand interacts with AChE with a low docking score of −11.5, indicating its essential affinity to interact with the receptor rather than other studied compounds. The binding mode for compound L15 can be summarized by 5-π-stacking interactions with TYR72, TYR124, TRP286, TYR341, and 8-hydrophobic interactions with TYR72, LEU76, GLU285, TRP286, PHE297, TYR337, PHE338, and TYR341. Figure 2 illustrates the binding mode of L14 and L15 in the active site of the AChE.

Discussion
Molecular dynamics
We used 100 ns molecular dynamics simulations to evaluate the stability of (L14, L15)-AChE complex. In this sense, binding interactions of the docking complexes were provided with a system embedded with water molecules, temperature, and pressure. As seen in Figure 3, the individual root mean square deviations (RMSDs) of the complexes were plotted on the simulation time, which indicates the converged L14-AChE and L15-AChE complexes to an equilibration state in the last 40 ns and 60 ns of simulation, respectively. Besides, stability during the simulation was also experienced. Regarding the binding free energies in the active site, quantitative estimates were obtained using the MM-PBSA free energy calculations. The energy values related to the 14 and 15 ligands at various time levels are presented in Table 4.

Compared to L15, the difference in the van der Waals energy levels reveals a more specific molecular shape for L14. This result shows that the L14 complex is more stable in the enzyme’s active site. Moreover, the results from MM-PBSA and RMSD show that the L14-AChE and L15-AChE complexes are highly stable. The amino acid fluctuations in the RMSF plot (Figure 4) reveal that, during the simulation, the following domains of 70-90, 255-265, 280-295, and 330-355 possessed the highest fluctuation. In sum, under NPT (National Standard Taper Pipe Threads Size NPT Chart) conditions, even though protein-ligand complexes were stabilized, the greatest effect on the enzyme came from L14. To further evaluate the conformational changes in AChE, we performed a Gyration radius (Rg) analysis to clarify the role of protein compactness (Figure 4b). As seen in both complexes, Rg of AChE at 0 ns was 2.34 nm depicting that, throughout the simulation time, we witnessed a conserved AChE conformation. 

Conclusion
In conclusion, the inhibitory activity of a group of bis-coumarin derivatives was tested to study novel candidates which would act as AChE inhibitors. The results showed that some of the tested compounds indicated acceptable activity against AChE. Besides, the L14 compound demonstrated significant activity, having an IC50 value of 2.0 μM. Moreover, the binding state of coumarin derivatives was rationalized using intra-silica studies regarding stability analysis, IC50 values, binding interactions, and binding score. The principal stabilizing forces were hydrogen bonds and halogen and hydrophobic interactions in AChE. Herein, it is suggested that these compounds can be regarded as strong candidates against AChE, used to treat Alzheimer disease. Moreover, it is worth mentioning that further studies must be performed to establish the efficacy and safety of these products.

Data availability
The datasets generated during and or analyzed during the current study are available from the corresponding author on reasonable request.

Ethical Considerations
Compliance with ethical guidelines
All ethical principles were considered in this study.

Funding
This research did not receive any grant from funding agencies in the public, commercial, or non-profit sectors.

Authors' contributions
Conceptualization and supervision: Zhila Zare Akbari; Methodology: Zhila Zare Akbari, Ladan Edjali, and Moosa Eshaghi; Writing the original draft, review, and editing: Zhila Zare Akbari, Moosa Eshaghi, and Ladan Edjali; Data curation, investigation, formal analysis, investigation: All authors. 

Conflict of interest
The authors declared no competing interests.

Acknowledgments
The authors wish to express their special thanks to Malihe Akhavan from Mazandaran University of Medical Sciences for providing laboratory facilities to carry out this research.


References

1. Jung M, Park M. Acetylcholinesterase inhibition by flavonoids from Agrimonia pilosa. Molecules. 2007; 12(9):2130-9. [DOI:10.3390/12092130] [PMID] [PMCID]
2. Al-Aboudi A, Odeh H, Khalid A, Naz Q, Choudhary MI, Atta-Ur-Rahman. Butyrylcholinesterase inhibitory activity of testosterone and some of its metabolites. J Enzyme Inhib Med Chem. 2009; 24(2):553-8. [DOI:10.1080/14756360802236393] [PMID]
3. Kozurkova M, Hamulakova S, Gazova Z, Paulikova H, Kristian P. Neuroactive multifunctional tacrine congeners with cholinesterase, anti-amyloid aggregation and neuroprotective properties. Pharmaceuticals. 2011; 4(2):382-418. [DOI:10.3390/ph4020382] [PMCID]
4. Zhao T, Ding KM, Zhang L, Cheng XM, Wang CH, Wang ZT. Acetylcholinesterase and butyrylcholinesterase inhibitory activities of β-carboline and quinoline alkaloids derivatives from the plants of genus Peganum. J Chem. 2013; 2013:717232. [DOI:10.1155/2013/717232]
5. Gordh T, Wahlin Å. Potentiation of the neuromuscular effect of succinylcholine by tetrahydro-amino-acridine. Acta Anaesthesiol Scand. 1961; 5(2):55-61. [DOI:10.1111/j.1399-6576.1961.tb00083.x] [PMID]
6. Shen ZX. Rationale for diagnosing deficiency of ChEs and for applying exogenous HuChEs to the treatment of diseases. Med Hypotheses. 2008; 70(1):43-51. [DOI:10.1016/j.mehy.2007.04.035] [PMID]
7. Loizzo MR, Tundis R, Menichini F, Menichini F. Natural products and their derivatives as cholinesterase inhibitors in the treatment of neurodegenerative disorders: An update. Curr Med Chem. 2008; 15(12):1209-28. [DOI:10.2174/092986708784310422] [PMID]
8. Hoerr R, Noeldner M. Ensaculin (KA-672. HCl): A multitransmitter approach to dementia treatment. CNS Drug Rev. 2002; 8(2):143-58. [DOI:10.1111/j.1527-3458.2002.tb00220.x] [PMID] [PMCID]
9. Fallarero A, Oinonen P, Gupta Sh, Blom P, Galkin A, Mohan CG, et al. Inhibition of acetylcholinesterase by coumarins: The case of coumarin 106. Pharmacol Res. 2008; 58(3-4):215-21. [DOI:10.1016/j.phrs.2008.08.001] [PMID]
10. Huang XY, Shan ZJ, Zhai HL, Su L, Zhang XY. Study on the anticancer activity of coumarin derivatives by molecular modeling. Chem Biol Drug Des. 2011; 78(4):651-8. [DOI:10.1111/j.1747-0285.2011.01195.x] [PMID]
11. Hwang JK, Noh EM, Moon SJ, Kim JM, Kwon KB, Park BH, et al. Emodin suppresses inflammatory responses and joint destruction in collagen-induced arthritic mice. Rheumatology. 2013; 52(9):1583-91. [DOI:10.1093/rheumatology/ket178] [PMID]
12. Vasconcelos JF, Teixeira MM, Barbosa-Filho JM, Agra MF, Nunes XP, Giulietti AM, et al. Effects of umbelliferone in a murine model of allergic airway inflammation. Eur J Pharmacol. 2009; 609(1-3):126-31. [DOI:10.1016/j.ejphar.2009.03.027] [PMID]
13. De Almeida Barros TA, De Freitas LAR, Filho JMB, Nunes XP, Giulietti AM, De Souza GE, et al. Antinociceptive and anti-inflammatory properties of 7-hydroxycoumarin in experimental animal models: Potential therapeutic for the control of inflammatory chronic pain. J Pharm Pharmacol. 2010; 62(2):205-13. [DOI:10.1211/jpp.62.02.0008] [PMID]
14. Mohammadi M, Khodamorady M, Tahmasbi B, Bahrami K, Ghorbani-Choghamarani A. Boehmite nanoparticles as versatile support for organic-inorganic hybrid materials: Synthesis, functionalization, and applications in eco-friendly catalysis. J Ind Eng Chem. 2021; 97:1-78. [DOI:10.1016/j.jiec.2021.02.001]

15. Kostova I, Bhatia S, Grigorov P, Balkansky S, Parmar VS, Prasad AK, et al. Coumarins as antioxidants. Curr Med Chem. 2011; 18(25):3929-51. [DOI:10.2174/092986711803414395] [PMID]
16. Rollinger JM, Hornick A, Langer T, Stuppner H, Prast H. Acetylcholinesterase inhibitory activity of scopolin and scopoletin discovered by virtual screening of natural products. J Med Chem. 2004; 47(25):6248-54. [DOI:10.1021/jm049655r] [PMID]
17. Iwaki S, Nakamura T, Koyama J. Inhibitory effects of various synthetic substrates for aminopeptidases on phagocytosis of immune complexes by macrophages. J Biochem. 1986; 99(5):1317-26. [DOI:10.1093/oxfordjournals.jbchem.a135599] [PMID]
18. Dinparast L, Hemmati S, Alizadeh AA, Zengin G, Samadi Kafil H, Bahadori MB, et al. An efficient, catalyst-free, one-pot synthesis of 4H-chromene derivatives and investigating their biological activities and mode of interactions using molecular docking studies. J Mol Struct. 2020; 1203:127426. [DOI:10.1016/j.molstruc.2019.127426]
19. Akhavan M, Foroughifar N, Pasdar H, Khajeh-Amiri AR, Bekhradnia AR. Copper(II)-complex functionalized magnetite nanoparticles: a highly efficient heterogeneous nanocatalyst for the synthesis of 5-arylidenthiazolidine-2,4-diones and 5-arylidene-2-thioxothiazolidin-4-one. Transit Met Chem. 2017; 42(6):543-52. [DOI:10.1007/s11243-017-0159-3]
20. Esam Z, Akhavan M, Bekhradnia AR. One-pot multicomponent synthesis of novel 2-(piperazin-1-yl) quinoxaline and benzimidazole derivatives, using a novel sulfamic acid functionalized Fe3O4 MNPs as highly effective nanocatalyst. Appl Organomet Chem. 2021; 35(1):e6005. [DOI:10.1002/aoc.6005]
21. Esam Z, Akhavan M, Bekhradnia AR, Mohammadi M, Tourani S. A novel magnetic immobilized para-aminobenzoic acid-Cu(II) Complex: A green, efficient and reusable catalyst for aldol condensation reactions in green media. Catal Letters. 2020; 150(11):3112-31. [DOI:10.1007/s10562-020-03216-w]
22. Akhavan M, Foroughifar N, Pasdar H, Bekhradnia AR. Green synthesis, biological activity evaluation, and molecular docking studies of aryl alkylidene 2, 4-thiazolidinedione and rhodanine derivatives as antimicrobial agents. Comb Chem High Throughput Screen. 2019; 22(10):716-27. [DOI:10.2174/1386207322666191127103122] [PMID]
23. Zare-Akbari Zh, Dastmalchi S, Edjlali L, Dinparast L, Es’haghi M. A novel nanomagnetic solid acid catalyst for the synthesis of new functionalized bis-coumarin derivatives under microwave irradiations in green conditions. Appl Organomet Chem. 2020; 34(7):e5649. [DOI:10.1002/aoc.5649]
24. Ellman GL, Courtney KD, Andres Jr V, Featherstone RM. A new and rapid colorimetric determination of acetylcholinesterase activity. Biochem Pharmacol. 1961; 7(2):88-90-IN1,91-5. [DOI:10.1016/0006-2952(61)90145-9]
25. Bourne Y, Kolb HC, Radić Z, Sharpless KB, Taylor P, Marchot P. Freeze-frame inhibitor captures acetylcholinesterase in a unique conformation. Proc Natl Acad Sci U S A. 2004; 101(6):1449-54. [DOI:10.1073/pnas.0308206100] [PMID] [PMCID]
26. Trott O, Olson AJ. AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem. 2010; 31(2):455-61. [DOI:10.1002/jcc.21334] [PMID] [PMCID]
27. Lindahl E, Hess B, van der Spoel D. GROMACS 3.0: A package for molecular simulation and trajectory analysis. Mol Model Annu. 2001; 7(8):306-17. [DOI:10.1007/s008940100045]
28. Berendsen HJC, van der Spoel D, van Drunen R. GROMACS: A message-passing parallel molecular dynamics implementation. Comput Phys Commun. 1995; 91(1-3):43-56. [DOI:10.1016/0010-4655(95)00042-E]
29. Van Der Spoel D, Lindahl E, Hess B, Groenhof G, Mark AE, Berendsen HJC. GROMACS: Fast, flexible, and free. J Comput Chem. 2005; 26(16):1701-18. [DOI:10.1002/jcc.20291] [PMID]