Molecular recognition of naphthalene diimide ligands by telomeric quadruplex-DNA: the importance of the protonation state and mediated hydrogen bonds

Abstract

In depth Monte Carlo conformational scans in combination with molecular dynamics (MD) simulations and electronic structure calculations were applied in order to study the molecular recognition process between tetrasubstituted naphthalene diimide (ND) guests and G-quadruplex (G4) DNA receptors. ND guests are a promising class of telomere stabilizers due to which they are used in novel anticancer therapeutics. Though several ND guests have been studied experimentally in the past, the protonation state under physiological conditions is still unclear. Based on chemical intuition, in the case of N-methyl-piperazine substitution, different protonation states are possible and might play a crucial role in the molecular recognition process by G4-DNA. Depending on the proton concentration, different nitrogen atoms of the N-methyl-piperazine might (or might not) be protonated. This fact was considered in our simulation in terms of a case by case analysis, since the process of molecular recognition is determined by possible donor or acceptor positions. The results of our simulations show that the electrostatic interactions between the ND ligands and the G4 receptor are maximized in the case of the protonation of the terminal nitrogen atoms, forming compact ND G4 complexes inside the grooves. The influence of different protonation states in terms of the ability to form hydrogen bonds with the sugar–phosphate backbone, as well as the importance of mediated vs. direct hydrogen bonding, was analyzed in detail by MD and relaxed force constant (compliance constant) simulations.

---

Introduction

Under physiological conditions, certain G-rich oligonucleotide sequences fold into structures called G-quadruplexes¹ (G4). These structures are formed by planar guanine quartets held together by Hoogsteen hydrogen bonds and π–π stacking interactions. Nucleotide sequences that are able to fold into G4 conformations were found in many regions of the human genome, for example in telomeres and oncogenes^2,3 and recently they were also observed in living cells.^4,5 The telomeric endings of human chromosomes for example are composed by a repetition of a noncoding hexamer (TTAGGG) protecting the chromosomes from recombination and degradation.⁶ In human somatic cells the whole telomeric DNA has a length of ca. 15 kb,⁷ with a short (150–200 nucleotides) single stranded 3′ termination. Normally, the telomeric DNA is progressively shortened during cellular replication, a process which is controlled and counter balanced by the enzyme telomerase by adding nucleotides to the end of telomeres. Nevertheless, when a living cell reaches its critically short length, it enters the senescence state and stops replication. In contrast, 80–90% of human tumor cells are characterized by telomerase overexpression. Since the mid-1990s, one anti-cancer strategy has indeed been to use small molecules able to stabilize G-quadruplex structures by means of telomerase inhibition. Consequently, a number of ligands which are able to bind selectively to the G-quadruplex motif have been studied since then.^8–11 One of the first examples, 2,6-diamidoanthraquinone, stabilizing the G-quadruplex by a “threading” intercalation mode, was published as early as 1997.⁹ In 1998, it was discovered that perylene diimides (for example, N,N′-bis[2-(1-piperidino)ethyl]-3,4,9,10-perylenetetracarboxylic diimide, also known as PIPER) are able to induce and stabilize DNA G-quadruplex structures due to which they inhibit telomerase activity too.¹⁰ Not much later, in 2001, Quarfloxin (CX-3543), a fluoroquinolone derivative, was introduced, which was in succession optimized by Cylene Pharmaceuticals. It binds to G-quadruplex DNA and selectively disrupts the interaction of rDNA G-quadruplexes with the nucleolin protein, thereby inhibiting Pol I transcription.¹¹ To date CX-3543 is the sole quadruplex-binding compound, which has entered the clinical phases I an II.

Quite recently, piperazine substituted naphthalene diimide (ND) guests, pioneered by the Neidle group, have attracted special attention. They are characterized by a common naphthalene diimide aromatic core, which is able to stack effectively to the G4 3′ quartet, while the piperazine substituents strengthen the recognition process via hydrogen bonds. Under physiological pH-conditions, the piperazine nitrogen atoms are protonated, leading to four positively charged ND arms. Experimental^12–14 and computational¹⁵ biophysical studies could demonstrate that such ligands are in indeed effective G4-binders.

Though it is not always considered explicitly, the protonation state of a ligand^16–18 and the host¹⁹ plays a key role in drug design and molecular recognition. For example, adding a proton to a neutral ligand leads to different binding modes, enhancing the hydrogen bonding aptitude with the negatively charged phosphate groups. In the case of ND G4 complexes, the strength of the non-covalent interactions, and as a result the stability of the complex, heavily depends on the protonation state, which is unknown for most solid-state PDB structures. For example, in order to design novel small-molecule ligands, it is therefore of crucial interest to (1) understand and, in a second step, (2) control the guest's protonation state by the design of novel candidates with tailored substituents. In this respect, a thorough in silico simulation might help in finding effective new binders, which need, of course, to be tested in vitro and finally in vivo.

As part of our ongoing studies on the recognition process of G4 telomeres,²⁰ here we present a computational study concerning the binding of a series of naphthalene diimide (ND) ligands 1–4 (Fig. 1) using a well-defined telomeric quadruplex model.

Fig. 1 Naphthalene diimide ligands 1–4, where n = 2–5 is the number of CH₂ groups in the chain.

Enthalpic binding strengths

In general, the methyl-piperazine moiety can be considered as a diprotic base. For 1,4-dimethylpiperazine the reported values of pK_a1 and pK_a2 at 298 K are 8.38 and 3.81, respectively.²¹ Two limiting protonation states were therefore designed in our study: (1) the internal nitrogen atoms of the four methyl-piperazine rings are protonated (termed intra), and (2) the terminal piperazine nitrogen atoms (termed ter) are protonated. In both cases the total charge of the protonated guests is +4. The principal interactions of the ND ligands with G-quadruplex DNA are shown in Fig. 2.

Fig. 2 The most important interactions for the molecular recognition process (ND2/G4 complex). Electrostatic (phosphate oxygen atoms in red) and hydrogen bond interactions (left). π–π stacking interactions between the naphthalene diimide core and the telomeric quadruplex (right).

Two structures published by Parkinson et al.¹³ were used as the starting point of our model construction. The parallel telomeric quadruplex was taken from the pdb entry 3SC8. Our ligands were modeled starting from both crystallographic structures, hydrogen atoms were added and basic amines were protonated according to our definition of the intra and ter protonation states, by hand: ligand ND1 was modeled, starting from 3SC8, and a methylene group was deleted. Ligand ND2 was taken from 3SC8 as it was. Ligand ND3 was adopted from 3T5E. Ligand 4 was modeled from 3T5E, adding one additional methylene group in this case. Crystallographic water molecules and the two internal potassium ions were retained in our liquid state simulation. The OPLS_2005 force field^22,23 and the generalized Born/solvent accessible surface (GB/SA) continuum solvation model²⁴ were used applying the Schroedinger software package.²⁵ The conformational space was scanned applying a Monte Carlo Multiple Minimum (MCMM) algorithm as implemented in the MacroModel,²⁶ while the Polak–Ribiere conjugate gradient (PRCG) algorithm was used for minimization.

After a preliminary conformational Monte Carlo search of 90 000 steps, the resulting liquid phase global minimum of the G4 receptor together with the computed OPLS_2005 energy (E_G4) represented the model for our subsequent ND binding study. Each telomeric ND(1–4) complex was generated using the following protocol: (1) the individual ND guests were posed 8–9 Å away from the receptor, keeping the diimide core parallel to the 3′ quartet of the telomere. (2) 20 000 minimization steps (steepest descent) were conducted. (3) Finally, 120 000 Monte Carlo steps were performed in order to scan the low energy region. In this final step, the internal potassium ions were constrained with a force constant of 100 kJ mol⁻¹. All other degrees of freedom were included as search parameters. The relative enthalpic interaction strengths of the ND guests (ΔE) were simply calculated as the difference between the global minimum energy of the complex (E_G4–ND), and the sum of global minima of the G4 receptor (E_G4) and the individual isolated (liquid phase) ND(1–4) guests (E_ND1; E_ND2; E_ND3; E_ND4). Entropic contributions were ignored in this part of our study. In Table 1, our computed ΔE values for the binding of ligands 1–4 are compared with the experimental melting data published by Neidle et al.¹² A higher ΔT_m value is connected with a higher stability of the complex.

Table 1 Comparison of the computed trend for the ND guests 1–4 in terms of the calculated (OPLS_2005) enthalpic binding energy (ΔE in kJ mol⁻¹) with experimental melting data (in °C) from Table 1 in ref. 13 [G. W. Collie, R. Promontorio, S. M. Hampel, M. Micco, S. Neidle and G. N. Parkinson, J. Am. Chem. Soc., 2012, 134, 2723–2731]. Both protonation states, intra and ter, were considered in our simulation

	ND1	ND2	ND3	ND4
ΔTm ^13	19.0	28.3	24.7	23.8
ΔE intra	−274.7	−252.6	−314.9	−314.0
ΔE ter	−294.5	−420.1	−405.8	−388.2

---

Looking at Table 1, we come to the following conclusion. While protonation of the terminal protonation state is leading to a stronger enthalpic binding for all ND ligands, the experimental stability trend ND2 > ND3 > ND4 > ND1 is reproduced by our simulation of the ter state (terminal protonation) only, which seems to be the relevant protonation state under the experimental conditions. Our computations therefore seem to confirm that three to four (ND2, ND3) methylene spacer groups represent the optimal chain length. While the enthalpic penalty due to suboptimal orientation of the two methylene groups (ND1) is dramatic (>100 kJ mol⁻¹ relative to the ND2 complexation), the addition of more and more spacer methylene groups seems to moderately hamper the binding strength. Since our simulation seems to prove that (1) the terminal protonation is the relevant state under the experimental conditions and (2) the naphthalene diimide ligand ND2 is the best binder in our series, we focus on this very complex in the following.

In order to test the robustness of our results in terms of the relative proton affinity, we calculated the energy of the two protonation states for the most effective binder, ND2, using electronic structure methods (DFT; B3LYP-D3 functional and a standard 6-31g(d) basis set), starting our DFT minimization from the OPLS_2005 global minimum conformer geometry. The DFT gas phase calculation seems to reproduce our force field simulation of the liquid phase. In fact, the terminal protonation of the guest ND2 is preferred by −146.7 kJ mol⁻¹. We checked the damping effect of the solvents, running a single point calculation including implicit water (PCM) for ND2. As expected, the energy difference ter/intra drops down to −21.0 kJ mol⁻¹. Nevertheless, the trend seems to be robust: the terminal protonation is preferred over the internal.

Strengths and dynamics of the individual non-covalent interactions

In order to compare the binding mode in the solid state and in solution, we compared the crystal structure 3SC8 with the global minimum found by our conformational Monte Carlo scan in the liquid phase (ter protonation) for guest ND2 (Fig. 3).

Fig. 3 Superposition of the crystal structure (light blue) with the global minimum found for ligand ND2 ter (red).

Besides the overall analogy, there are some differences in the stacking position of the ligand, in particular the position of the ND guest is closer to the 3′ quartet and shifted horizontally in the solid state. We speculate that this difference is due to crystal packing effects. In the solid state, the ND2 guest interacts with two quadruplex receptors at the same time, shifting the minimum position of the diimide core. The same seems to be true for the second protonation state (intra; Fig. 4).

Fig. 4 Superposition of the crystal structure (light blue) with the global minimum found for ligand 2 intra (red).

In order to analyze the individual contributions of the individual non-covalent interactions to the overall recognition process, we computed the strength of the most important non-covalent interactions between the receptor and ligand 2 by means of their relaxed force constants (inverse compliance constants). For an overview of the complex network of all mediated hydrogen bonds, see Fig. S1 (ESI†). The values obtained are reported in Table 2.

Table 2 Relaxed force constants (OPLS_2005 force field in N cm⁻¹) of the hydrogen bonds between ligands and the telomeric receptor (ND2/G4) complex. In the case of terminal protonation state, there are two direct NH⋯O(phosphate) and two water mediated NH⋯O(water) hydrogen bonds, while the intra protonation state is characterized by water mediated NH⋯O(water) hydrogen bonds exclusively (black). The values of the hydrogen bonds leading from the mediating water to the phosphate oxygen are in italics

ND2 intra complex		ND2 ter complex
NH⋯O(water)	0.50	NH⋯O(phosphate)	0.29
NH⋯O(water)	0.38	NH⋯O(water)	0.42
NH⋯O(water)	0.59	NH⋯O(phosphate)	0.45
NH⋯O(water)	0.34	NH⋯O(water)	0.53
Water⋯O(phosphate)	0.09	Water⋯O(phosphate)	0.19
Water⋯O(phosphate)	0.03	Water⋯O(phosphate)	0.14
Water⋯O(phosphate)	0.12	—	—
Water⋯O(phosphate)	0.15	—	—

---

In the case of terminal protonation state ter, the global minimum obtained from the conformational search shows (I) two direct NH⋯O(phosphate) and (II) two water mediated NH⋯O(water) hydrogen bonds, while the intra protonation state is characterized by water mediated NH⋯O(water) hydrogen bonds exclusively. Fig. 5 depicts the two kinds of hydrogen bonds observed in the global minimum of the 2 ter complex. Looking at Table 2, in terms of their stiffness (and lifetime), the direct piperazine/phosphate hydrogen bonds (N–H⋯O–phosphate) seem to be more or less of the same quality as the mediated bonds (N–H⋯O–water). The calculated relaxed force constants for these kinds of non-covalent interactions range between approximately 0.3 and 0.6 N cm⁻¹. In contrast, the calculated softness of the hydrogen bonds between the mediating water and the accepting phosphate oxygen water⋯O(phosphate) is significantly lower: the values between 0.1 and 0.2 N cm⁻¹ point to very short hydrogen bond lifetimes.

Fig. 5 Direct (left) and water mediated (right) hydrogen bonds observed for the computed (OPLS_2005) lowest energy conformer of the complex between the telomeric G4 quadruplex receptor and the terminal protonated (ter) naphthalene diimide ligand ND2.

The static picture of the relevant hydrogen bond net might be misleading in this case. In order to explain the energy difference between the two protonation states, we therefore analyzed the general radial distribution function of ND2 intra and ter guests complexed with the G4 receptor (Fig. 6).

Fig. 6 Radial distribution functions (left) and related integrals (right) of the charged nitrogen and phosphorous atoms for ligand 2. (This general trend is also visible from the liquid phase simulation of the isolated ND guests 1, see ESI,† Table S2.)

In its ter protonation state, guest ND2 is characterized by short range interactions, showing two direct hydrogen bonds. (For a comparison with the shorter ND1 guest, see the ESI,† Fig. S2.) The main difference nevertheless is visible in the long range interactions: the average positions of the terminal rings are closer to phosphate groups. The methyl-piperazine rings are locked deep inside the quadruplex grooves, leading to a stronger overall electrostatic interaction with the receptor. The different positions of the charged arms are clearly shown by the electrostatic potential surfaces for the global minimum conformer (Fig. 7), where the negatively charged phosphates are plotted in red and the positively charged nitrogen atoms in blue, respectively.

Fig. 7 The electrostatic potential for the intra (left) and ter (right) isomers of the G4 complex with ND2.

In order to further test our hypothesis concerning the differences of the hydrogen bond lifetimes in both protonation states, we performed molecular dynamics simulations (see ESI† for details of our setup). The results are summarized in Fig. 8. As can be easily seen, in strong contrast to the internal protonation, the terminal protonation state is stabilized by a larger number of direct NH⋯OP hydrogen bonds during the simulation time. The 50 ns trajectory of the terminal protonation state is characterized by at least one strong hydrogen bond, while the second and third hydrogen bonds are auxiliary by nature. During the whole 50 ns seconds of simulation, we don't see a single time step, characterized by four simultaneous hydrogen bonds.

Fig. 8 Here we have plotted the number of NH⋯OP hydrogen bonds versus the simulation time for the intra (left) and the ter (right) isomers of ligand 2.

Conclusions

Due to our combination of experimental data, thorough Monte Carlo conformational scans, MD simulations, DFT calculations and simulations of the mechanical binding strengths, we have analyzed (A) different protonation states of tetra-substituted naphthalene diimide (ND) guests complexed with G-quadruplex (G4) DNA receptors, and (B) the importance of water mediated vs. direct hydrogen bonding. Based on the reproduction of the relative stabilities of a series of ND ligands published by the Neidle group we come to the following conclusions: (1) terminal protonation of the naphthalene diimide seems to be the dominant state, leading to persistent, direct NH⋯OP bonds. (2) Due to our calculation of the relaxed force constants in combination with MD simulations, the direct bonding is by far stronger and longer lasting than their water mediated counterparts. In this case the static picture of the X-ray structure might be misleading. (3) In the liquid phase, besides hydrogen bonding, the longer range electrostatic interactions are maximized in the terminal protonation state, leading to a perfect fit of the methyl-piperazine rings which are locked deep inside the quadruplex grooves in this case.

Computational details

The Gaussian09 software²⁷ was used for quantum mechanical and AMBER molecular mechanics (MM) calculations. Geometry optimizations were conducted using a B3LYP functional and a 6-31g(d) basis set. The Amber 99 force field was used.²⁸ The COMPLIANCE software and its graphical interface^29,30 were used to evaluate the strength of the most important hydrogen bonds between the ligands and the quadruplex receptor.

Molecular dynamics (MD) simulations on both isomers of ligand 2 were performed using the GROMACS 4.5.3 software package^31,32 and the Amber99 force field³³ with Parmbsc0 nucleic acid torsions.^34,35 The partial atomic charges of ligand 2 were obtained by DFT calculations, while other intramolecular force-field parameters were generated using the ACPYPE software.^36–39 A triclinic box of TIP3P water molecules was added around the system to a depth of 1.5 nm on each side of the solutes. 15 K⁺ counterions were added to neutralize the negative charges of the phosphate groups, while other 22 K⁺ and Cl⁻ ions were added to set the solution ionic strength to about 0.15 M. The two internal potassium ions were kept from the Monte Carlo simulations. Explicit solvent simulations were performed within the isothermal–isobaric NPT ensemble, at a temperature of 300 K, under control of a velocity rescaling thermostat.³⁹ The particle mesh Ewald method⁴⁰ was used to describe long-range electrostatic interactions. The timestep for integration was 2 fs and all covalent bonds were constrained using the LINCS algorithm. Preliminary energy minimizations were run for 5000 steps using the steepest descend algorithm. During equilibration, systems were harmonically restrained with a force constant of 1000 kJ mol⁻¹ nm⁻², gradually relaxed into five consecutive steps of 100 ps each, to 500, 200, 100 and 50 kJ mol⁻¹ nm⁻². The MD simulations were carried out for 50 ns. The hydrogen bond networks were analyzed using the GROMACS tool g_hbond.

References

1. S. A. Ohnmacht and S. Neidle, Bioorg. Med. Chem. Lett., 2014, 24, 2602–2612 , and references herein.
2. A. K. Todd, M. Johnston and S. Neidle, Nucleic Acids Res., 2005, 33, 2901–2907.
3. J. L. Huppert and S. Balasubramanian, Nucleic Acids Res., 2007, 35, 406–413.
4. G. Biffi, M. Di Antonio, D. Tannahill and S. Balasubramanian, Nat. Chem., 2014, 6, 75–80.
5. G. F. Salgado, C. Cazenave, A. Kerkour and J.-L. Mergny, Chem. Sci., 2015, 6, 3314–3320.
6. J. A. Hackett, D. M. Feldser and C. W. Greider, Cell, 2001, 106, 275–286.
7. T. R. Cech, Angew. Chem., Int. Ed., 2000, 39, 34–43.
8. For a review, see: Y. Xu, Chem. Soc. Rev., 2011, 40, 2719–2740.
9. D. Sun, B. Thompson, B. E. Cathers, M. Salazar, S. M. Kerwin, J. O. Trent, T. C. Jenkins, S. Neidle and L. H. Hurley, J. Med. Chem., 1997, 40, 2113–2116.
10. C. Sissi, L. Lucatello, A. P. Krapcho, D. J. Maloney, M. B. Boxer, M. V. Camarasa, G. Pezzoni, E. Mentac and M. Palumbo, Bioorg. Med. Chem., 2007, 15, 555–562.
11. D. Drygin, et al. , Cancer Res., 2009, 69, 7653–7661.
12. M. Micco, G. W. Collie, A. G. Dale, S. A. Ohnmacht, I. Pazitna, M. Gunaratnam, A. P. Reszka and S. Neidle, J. Med. Chem., 2013, 56, 2959–2974.
13. G. W. Collie, R. Promontorio, S. M. Hampel, M. Micco, S. Neidle and G. N. Parkinson, J. Am. Chem. Soc., 2012, 134, 2723–2731.
14. S. M. Hampel, A. Sidibe, M. Gunaratnam, J.-F. Riou and S. Neidle, Bioorg. Med. Chem. Lett., 2010, 20, 6459–6463.
15. G. Prato, S. Silvent, S. Saka, M. Lamberto and D. Kosenkov, J. Phys. Chem. B, 2015, 119, 3335–3347.
16. S. Kalyaanamoorthy and Y.-P. P. Chen, J. Mol. Graphics Modell., 2013, 44, 44–53.
17. S. Zhu, M. F. Brown and S. E. Feller, J. Am. Chem. Soc., 2013, 135, 9391–9398.
18. D. P. Oehme, R. T. C. Brownlee and D. J. D. Wilson, J. Comput. Chem., 2012, 33, 2566–2580.
19. S. K. Sadiq, D. W. Wright, O. A. Kenway and P. V. Coveney, J. Chem. Inf. Model., 2010, 50, 890.
20. J. Grunenberg, G. Barone and A. Spinello, J. Chem. Theory Comput., 2014, 10, 2901–2905.
21. F. Khalili, A. Henni and A. L. L. East, J. Chem. Eng. Data, 2009, 54, 2914–2917.
22. W. L. Jorgensen, D. S. Maxwell and J. Tirado-Rives, J. Am. Chem. Soc., 1996, 118, 11225–11236.
23. J. L. Banks, H. S. Beard, Y. Cao, A. E. Cho, W. Damm, R. Farid, A. K. Felts, T. A. Halgren, D. T. Mainz, J. R. Maple, R. Murphy, D. M. Philipp, M. P. Repasky, L. Y. Zhang, B. J. Berne, R. A. Friesner, E. Gallicchio and R. M. Levy, J. Comput. Chem., 2005, 26, 1752–1780.
24. W. Hasel, T. F. Hendrickson and W. C. Still, Tetrahedron Comput. Methodol., 1988, 1, 103–116.
25. Maestro, version 9.2, Schrödinger, LLC, New York, NY, 2011.
26. MacroModel, version 9.9, Schrödinger, LLC, New York, NY, 2011.
27. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, Ö. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford, CT, 2009.
28. W. D. Cornell, P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell and P. A. Kollman, J. Am. Chem. Soc., 1995, 117, 5179–5197.
29. K. Brandhorst and J. Grunenberg, Chem. Soc. Rev., 2008, 37, 1558–1567.
30. K. Brandhorst and J. Grunenberg, J. Chem. Phys., 2010, 132, 184101.
31. D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C. Berendsen, J. Comput. Chem., 2005, 26, 1701–1718.
32. B. Hess, C. Kutzner, D. van der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
33. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
34. A. Pérez, I. Marchán, D. Svozil, J. Sponer, T. E. Cheatham 3rd, C. A. Laughton and M. Orozco, Biophys. J., 2007, 92, 3817–3829.
35. A. T. Guy, T. J. Piggot and S. Khalid, Biophys. J., 2012, 103, 1028–1036.
36. A. W. S. da Silva and W. F. Vranken, BMC Res. Notes, 2012, 5, 367.
37. J. Wang, W. Wang, P. A. Kollman and D. A. Case, J. Mol. Graphics Modell., 2006, 25, 247–260.
38. J. Wang, P. Cieplak and P. A. Kollman, J. Comput. Chem., 2000, 21, 1049–1074.
39. G. Bussi, D. Donadio and M. Parinello, J. Phys. Chem., 2007, 126, 14101–14107.
40. T. Darden, D. York and L. Pedersen, J. Phys. Chem., 1993, 98, 10089.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5cp05576h

---