Thermodynamics and Kinetics of Nucleobase Stacking

Article pubs.acs.org/JCTC

Thermodynamics and Kinetics of Nucleobase Stacking Oligomerization Revealed by Molecular Dynamics Simulations Fabian Zeller and Martin Zacharias* Physik-Department T38, Technische Universität München, James-Franck-Str. 1, 85748 Garching, Germany ABSTRACT: The aggregation of N6N9-dimethyladenine in aqueous solution into stacked oligomers represents an important model system for oligomerization processes. Molecular dynamics simulations allowed us to extract statistically converged stacking thermodynamics as well as converged association and dissociation kinetics. The simulations confirm an oligomerization mechanism according to a random isodesmic stacking model without significant cooperative effects. Multiple oligomerization and dissociation events were used to characterize stacking and unstacking processes and intermediate states in atomic detail. Standard three-point rigid body water force field models lead to significantly accelerated dynamics and an underestimation of enthalpic and entropic contributions. The fully flexible SPCFW water model yielded close agreement with the experimental equilibrium constant, entropy contribution, and, importantly, kinetic observables. In addition, the performance of several other water models was systematically tested in stacking dimerization simulations indicating that the realistic prediction of thermodynamics and kinetics will likely benefit from the use of flexible water models.

1. INTRODUCTION Oligomerization processes are of fundamental importance in many areas of organic and physical chemistry and form the basis of many processes in biological systems, ranging from peptide amyloid formation to actin/tubulin filament formation.1−3 An ultimate goal of molecular simulations is to reliably follow such processes and study thermodynamics, kinetics, and underlying driving forces in atomic detail. Molecular dynamics (MD) simulations have been used extensively to study the free energy of drug binding to proteins4,5 or to elucidate the thermodynamics effects of mutations.6 More recently, they have also been applied to extract kinetics of drug molecule binding to receptor molecules.7 In general, sufficiently long simulations allow direct access to kinetic properties like mean ligand residence times in binding sites, residence times in intermediate states during protein folding, and many other related molecular details often inaccessible to experiments. However, little statistically relevant data on the reliability of kinetic predictions by MD simulations is available, and it is not clear if current force fields are sufficiently accurate for a realistic prediction of binding kinetics. As a model system for the study of an oligomerization process that involves multiple complex binding events, we used an aqueous solution of N6N9-dimethyladenine molecules. These compounds undergo reversible binding events in solution and form stacked oligomers. The reversible stacking of the adenine derivate was previously investigated experimentally using Nuclear Magnetic Resonance (NMR) spectroscopy8 and ultrasound absorption methods9 to extract both thermodynamic and kinetic quantities of oligomerization. On the basis of the random isodesmic binding model which assumes a single characteristic equilibrium constant for the stacking of monomers, it has been possible to extract an © 2017 American Chemical Society

equilibrium stacking constant by variation of the N6N9dimethyladenine concentration. In addition, kinetic dissociation and recombination rates could be obtained from relaxation time measurements at various N6N9-dimethyladenine concentrations9 which occur on time scales that are now accessible by MD simulations in explicit solvent. We performed multiple microsecond equilibrium simulations of an N6N9-dimethyladenine ensemble and many thousand dimerization simulations. Simulations at different temperatures allowed the accurate extraction not only of free energies of binding but also of the associated entropy and enthalpy. Since many binding and unbinding events could be observed, a statistically converged analysis of the stacking kinetics and validation of the random isodesmic model was possible. In an extensive comparison of different force field models for the aqueous solvent, accurate thermodynamic quantities and kinetic parameters for the process could only be obtained by a fully flexible water model (SPCFW). In addition to the accurate prediction of the thermodynamics and kinetics, the simulations also give insight into the molecular details of the mechanism of stacking and unstacking during the oligomerization process.

2. METHODS 2.1. Random Isodesmic Stacking Model. In the random isodesmic stacking model,8−10 which was used for the analysis of the experimental measurements, all stacking interactions are considered to be equivalent with one single equilibrium stacking constant K Ai Aj K = Ai + j (1) Received: February 13, 2017 Published: May 8, 2017 3005

DOI: 10.1021/acs.jctc.7b00150 J. Chem. Theory Comput. 2017, 13, 3005−3011

Article

Journal of Chemical Theory and Computation

Figure 1. Trajectory snapshots of the 2 μs 0.16 M N6N9-dimethyladenine simulation with the fully flexible SPFCW water model at 298 K, indicating spontaneous aggregation and dissociation of the nucleobases.

applied to bonds involving the N6N9-dimethyladenine hydrogen atoms, a 2 fs time step was used. For the simulations in which all bonds involving hydrogen atoms were flexible, a 1 fs time step was used (a time step of 1 fs was used also in the original publication of the SPCFW17 water model.) For minimization, the steepest descent method was employed. For equilibration, the Langevin thermostat20 was used with a collision frequency of 5 ps−1. In the subsequent simulations, the temperature was controlled by the Andersen thermostat,21 using a collision frequency of 5 ps−1. For a comparative simulation, the Berendsen thermostat22 was used with a coupling time of 1 ps. The pressure of 1 bar was regulated with the Berendsen barostat in all simulations, using a coupling time of 1 ps. The nucleobases’ center-of-mass (COM) separation was saved at intervals of 10 ps. 2.4. Bulk Simulations. For the bulk simulations, 25 N6N9dimethyladenine molecules were randomly placed in a box with an edge length of 50 Å, imposing a minimum intermolecule distance of 3 Å. The systems were solvated with ∼8000 TIP3P, SPC/E, or SPCFW water molecules; energy minimized for 1000 steps; and equilibrated at 298 or 313 K for 40 ps. The six simulations were then run for 2 μs. A stacking interaction was defined to be established for center of mass distances smaller than 5 Å. If more than two intermolecular distances were within the stacking threshold for one molecule, only the two closest interactions were defined as stacking interactions. The statistical errors of the autocorrelation functions and equilibrium concentrations were calculated by the block bootstrap method23 in the form of 95% confidence intervals, using a block size corresponding to the relaxation times of the respective simulations. The statistical errors of the relaxation times and rates correspond to exponential fits to the total lower/upper CIs of the stacking concentration Cs autocorrelation functions. The statistical errors of the fits of the isodesmic model are standard deviations resulting from a nonlinear least-square fit. All snapshots were prepared using Pymol.24 2.5. Dimerization Simulations. For each investigated water model and simulation setup, 5000 simulations were run. The starting structure of each simulation was generated by placing two N6N9-dimethyladenine structures with random relative orientation at a distance of rstart = 21 Å. The systems were solvated with ∼3000 water molecules and subsequently minimized for 1000 steps and equilibrated for 40 ps with a harmonic distance restraint of k = 5 kcal/mol/Å2, keeping the two molecules at a distance of rstart during the equilibration. The

and multimer concentrations Ai =

A1i K i − 1 ,

A1K < 1

(2)

where Ai is the concentration of multimers consisting of i residues, involving i − 1 stacks. Furthermore, we note ∞

C0 =

∞

∞

∑ Aii , Cm = ∑ Ai , Cs = ∑ Ai(i − 1) i=1

i=1

i=1

(3)

where C0 is the concentration of N6N9-dimethyladenine residues, Cm is the total concentration of multimers (including monomers), and Cs is the concentration of formed stacking interactions. From equilibrium observations alone, a sequential isodesmic stacking behavior, in which stacks can only be formed or broken at the ends of the multimers, cannot be distinguished from random isodesmic stacking behavior, in which stacks can be formed between any multimers or be broken at any position. The random isodesmic stacking model however predicts one single relaxation time τ according to 1/τ = CmKR + 2KD

(4)

where KR and KD are, respectively, the recombination and dissociation rate constants. The sequential isodesmic stacking model, in contrast, predicts a broad spectrum of relaxation times. Sound absorption measurements of the relaxation behavior found only one dominant relaxation time, indicating random isodesmic stacking behavior9 for N6N9-dimethyladenine and also for N6N6-dimethyladenine.10 2.2. Force Fields. The N6N9-dimethyladenine structure was adopted from the ff14SB11 adenosine residue library entry. The ff14SB force field was used for the N6N9-dimethyladenine molecule. The point charges were obtained by semiempirical QM calculations on the bcc level using the antechamber module of Amber14.12 The TIP3P,13 TIP4P,14 TIP5P,15 SPC/ E,16 and SPCFW17 water models were used as implemented in Amber14. The Lennard-Jones interactions with the nucleobase atoms were calculated according to the Lorentz−Berthelot combining rules. 2.3. Simulation Parameters. All simulations were carried out with Amber14.12 The short-range interactions were cut off at 9 Å; the long-range interactions were treated with the PME method using the standard Amber14 parameters. For simulations in which the N6N9-dimethyladenine hydrogen atoms were constrained with the SHAKE algorithm18 and modified by hydrogen mass repartitioning (HMR),19 a 4 fs time step was used. For simulations with only the SHAKE algorithm 3006


Article

Journal of Chemical Theory and Computation systems were then simulated until the center of mass separation of the two molecules was larger than rmax = 26 Å. This led to the simulation of at least 750 binding events for all simulation setups. The simulated trajectories within separations of 21 Å correspond to equilibrium conditions, as every trajectory that entered the corresponding phase space was simulated until leaving it. For the calculation of the equilibrium kinetic rates, the bound state was defined by nucleobase center of mass separations smaller than 5 Å. Molecules were defined to be in the unbound state between distances of 12 and 21 Å. To obtain kinetic rates and equilibrium constants corresponding to standard states, the residence times in the unbound state were scaled by a factor [1/(c0Vunbound)] that accounts for the volume accessible at the standard state concentration c0 of 1 M.25 The 95% confidence intervals of the kinetic rates and equilibrium constants of the dimerization simulations were calculated by standard bootstrapping, as the observables are based on independent simulations which are not correlated.

Figure 2. Average concentrations Ai (95% CI) of multimers of order i at a residue concentration of C0 = 0.16 M and a temperature of 298 K. Solid lines are fits of the isodesmic stacking model.

3. RESULTS AND DISCUSSION We simulated a system of 25 N6N9-dimethyladenine molecules at a concentration of C0 = 0.16 M at 298 and 313 K, which is within the concentration and temperature range of the corresponding experiments. In these simulations, we employed the rigid body TIP3P13 and SPC/E16 as well as the fully flexible SPCFW17 water model. In the rigid body water models the bonds and bond angles are constrained to the corresponding equilibrium values.18 In the simulations, starting from a random arrangement of the 25 monomers, multimers form and dissociate spontaneously multiple times (Figure 1). For a direct comparison of the simulations with the experiments, we used the random isodesmic stacking model. A stacking interaction was defined to be established at a nucleobase center-of-mass distance of less than 5 Å, which is approximately in the sensitivity range of NMR distance measurements. While in experiments, individual concentrations of multimers Ai of different sizes can usually not be observed and several experiments at different nucleobase concentrations C0 have to be performed, in MD simulations, a direct extraction from the simulation at only one nucleobase concentration C0 is possible. The equilibrium constant K was determined by fitting the isodesmic model to the average multimer concentrations (Figure 2). Observation of the cumulative average multimer concentration Cm (Figure 3) shows excellent convergence of the equilibrium properties within the simulation time of 2 μs. As in the experiments, despite its simplicity, the isodesmic model appears to describe the equilibrium distribution of the different multimer orders reasonably well. We observe a slight underestimation of the monomer concentration by the isodesmic stacking model in the case of the SPCFW water model. However, no significant (anti)-cooperativity could be detected, which is in agreement with the experiments.9 While the TIP3P and to a lesser degree the SPC/E water models underestimate the equilibrium stacking constant K, at 298 K, the SPCFW water model predicts a value between the experimental result of an NMR study and a study using vapor pressure osometry (Table 1). To obtain the enthalpy and entropy contributions to the stacking constant K, we performed the simulations also at 313 K. As in the corresponding experiments, entropy and enthalpy were separated using the van’t Hoff equation. While the rigid

Figure 3. Cumulative average multimer concentration Cm (95% CI) at 298 K, indicating excellent convergence of the equilibrium characteristics.

body models underestimate the enthalpy as well as the entropy changes, the contributions are well predicted by the flexible SPCFW model. In particular, the probability distribution of bond angles in the SPCFW model might change upon binding (for released waters and waters surrounding the complex) which can affect entropic and enthalpic contributions. The more realistic description of the degrees of freedom of the water molecules might be essential for an accurate representation of the entropy/enthalpy contribution to the stack formation. To extract kinetic rates experimentally, the relaxation behavior of the stack concentration was measured by ultrasound absorption. From the MD simulations, the relaxation time was extracted using the autocorrelation function of the stack concentration. As observed in the sound-absorption experiment,9 the relaxation behavior can be reasonably well described by a single relaxation time (Figure 4). We observe a quasi instantaneous decay of the stack concentration autocorrelation within the time scale of the trajectory time resolution of 10 ps to about 0.6−0.8. We attribute this decay to fast fluctuations of the stacking concentrations due to the numerical stacking threshold of 5 Å and therefore used a scaling factor for the exponential fit. The relaxation time and the rates predicted by the TIP3P model are up to 1 order of magnitude faster than the experimental results (Table 1). The SPC/E model performs slightly better, but also significantly overestimates the stacking rates. In contrast, the relaxation time resulting from the SPCFW model leads to a remarkably close agreement of calculated kinetic rates compared to experimental values. An overestimation of the rates is expected for the TIP3P model, 3007


Article


Table 1. Thermodynamic Quantities Obtained from 25 Nucleobase Ensemble MD Simulations with the TIP3P, SPC/E, and the SPCFW Water Modelsa TIP3P

SPC/E

SPCFW

osom./sound9

NMR8

−7.3 ± 0.9 −17.0 ± 3.0 298 K 41.6 ± 1.1 4.9 [4.2, 6.2] 16 [12, 20] 3.8 [2.8, 4.6] 313 K 23.1 ± 0.5 2.7 [2.3, 3.8] 21 [14, 26] 9.4 [6.3, 11]

−8.7 ± 1.5 −21.6 ± 3.0

−7.2 ± 0.6 −17.9 ± 1.8

45.6

23.6

ΔH ΔS

[kcal/mol] [cal/mol/K]

−5.3 ± 0.1 −12.8 ± 0.5

−5.3 ± 0.3 −13.0 ± 1.1

K τ kR kD

[1/M] [ns] [108/Ms] [107/s]

12.3 ± 0.1 0.9 [0.8, 1.2] 46 [33, 54] 38 [28, 44]

15.1 ± 0.2 1.2 [1.0, 1.6] 38 [28, 45] 25 [19, 30]

K τ kR kD

[1/M] [ns] [108/Ms] [107/s]

8.0 ± 0.1 0.8 [0.7, 1.1] 41 [29, 49] 51 [37, 61]

9.8 ± 0.1 0.9 [0.8, 1.2] 41 [30, 48] 42 [31, 49]

9.3 5.0 18.0

15.0

20 6.6

a

Comparison to equilibrium quantities obtained from vapor pressure osometry,9 NMR spectroscopy,8 and kinetic data obtained by sound absorption measurements.9 All measured quantities are based on the random isodesmic stacking model. Simulations were performed at a N6N9dimethyladenine concentration of C0 = 0.16 M; experiments were performed within a residue concentration range of C0 = 0.05−0.5 M. Values in square brackets are statistical 95% confidence intervals, the errors for the equilibrium values are standard deviations as resulting from fits of the isodesmic stacking model.

Figure 4. Autocorrelation function AF of the stack concentration Cs at a nucleobase concentration of C0 = 0.16 M and 298 K. Dashed lines are fits to an exponential function cet/τ with scaling factor c.

because it also overestimates the self-diffusion coefficient by more than a factor of 2 (Table 2). Surprisingly, the SPC/E Table 2. Self-Diffusion Coefficient Ds for the Different Water Models at Approximately Standard Conditions, Reported by Previous MD Studies and Experiments Ds [10−5 cm2 s−1]

TIP3P

TIP4P

TIP5P

SPC/E

SPCFW

exp.

5.3017

3.3126

2.6226

2.4117

2.3217

2.327

model, which yields a self-diffusion coefficient close to the experimental measurement, shows only modest improvement with respect to the stacking kinetics compared to TIP3P. The overestimation of the recombination rates by the rigid body water models is therefore unlikely to be caused exclusively by diffusive effects. Note that the stability of solvent cages around solutes or barriers to break solvent structures in order to initiate association can be modulated by the water flexibility and influence dissociation and association kinetics. The simulations allow for a detailed analysis of the stacking configurations. Figure 5 shows that relative fluctuations of the nucleobases mainly occur along translations or rotations parallel to the base plane. Principal stacking and dissociation pathways are likely to follow these low free energy regions. The relative orientation of two stacked bases is dominated by the positions

Figure 5. Configuration probability histograms from the simulations with the SPCFW water model at 298 K projected on relative (A) vertical (z) and horizontal (r) translation and (B) torsion (ϕ) and tilt (ϵ) rotation of two stacking bases. Trajectory snapshots illustrate prevalently populated configurations.

of the methyl groups. Therefore, no significant population of a stacking configuration resembling the arrangement in B-DNA (about 34° relative rotation) was observed. The three water models used for the bulk simulations showed significant differences with respect to equilibrium and kinetic properties. For a more detailed comparison of a larger set of water models and the possible influence of different simulation techniques like HMR or SHAKE, we chose to limit 3008


Article


Figure 6. (A) Calculated free energy change vs distance between two N6N9-dimethyladenine molecules using different water models. The bound state was here defined by base distances smaller than 5 Å and the unbound state by nucleobase distances larger than 12 Å. For each water model, 5000 trajectories were started at a base−base distance of 21 Å and terminated when the distance became larger than 26 Å. This led to the simulation of at least 750 stacking events, yielding statistically converged binding and dissociation rates. The convergence of the calculated equilibrium binding constant, dissociation rate, and association rate is illustrated for the different water models in B−D, respectively. The indicated MD time step corresponds to different simulation conditions: HMR and SHAKE applied to nonwater hydrogen atoms (4 fs), SHAKE applied to nonwater hydrogen atoms (2 fs), fully flexible nonwater hydrogen atoms (1 fs). All water models except SPCFW are rigid body water models.

where N is the number of observed transitions, Tbulk is the total time spent in the standard state, and Tbound is the total time spent in the bound/stacked state. Different MD-simulation time steps were employed corresponding to different simulation conditions. Simulations with Hydrogen Mass Repartitioning (HMR)19 and constraints on covalent bonds involving hydrogens in the nucleobases (SHAKE18) allowed for a time step of 4 fs, whereas a time step of 2 fs was used in the case of using only SHAKE and omitting HMR. For simulations with unmodified nucleobase hydrogen atoms, a time step of 1 fs was used. Figure 6A shows the potential of mean force (PMF), equilibrium binding constants and on- and off-rates of N6N9dimethyladenine dimerization for the whole set of tested water models and simulation parameters. The cumulative average values indicate statistical convergence of the thermodynamic properties within 5000 simulations (observing ca. 750 binding events). The dimerization simulations TIP3P 4 fs, SPC/E 4 fs, and SPCFW 1 fs indicate the same relative behavior as has been found in the oligomerization simulations described above (Figure 6B,C,D). The use of the Berendsen thermostat instead of the Andersen thermostat did not lead to different rates within the statistical uncertainties for the TIP3P model. The application of HMR to the nucleobase hydrogen atoms slightly accelerates the kinetics in combination with the TIP3P model. This effect was not observed for the SPC/E water model. The influence of constraining the bond lengths of bonds involving

the simulations to the dimerization process. Instead of a single bulk ensemble simulation, independent trajectories with only two nucleobases were simulated. These simulations were started at base distances of 21 Å and stopped at base distances of 26 Å. Using the same distance threshold for stack formation, the equilibrium constants for N6N9-dimethyladenine dimerization are significantly higher for TIP3P, SPC/E, and SPCFW than the equilibrium constants resulting from the isodesmic model bulk simulations. It is important to realize that this is not a contradiction to the observed equivalence of different stacking interactions according to the isodesmic stacking model, as it can be attributed to the fact that in the bulk, beyond the threshold of 5 Å, favorable interactions with other nucleobases exist that lower the free energy of the unstacked state (Figure 6, attraction up to ca. 10 Å). It indicates, however, that the equilibrium constant resulting from the conditions in the bulk experiments and simulations are not concentration independent. Thus, as the dimerization simulations are not directly comparable to the bulk system, for a relative comparison of the water models we analyze the dimerization standard binding free energy25 and corresponding kinetic rates (stacked state closer than 5 Å, unstacked state larger than 12 Å, where attractive interactions are negligible, see Figure 6). The kinetic dimerization rates are now defined according to28 kon = N /Tbulk , koff = N /Tbound

(5) 3009


Article


(2) Straub, J. E.; Thirumalai, D. Toward a molecular theory of early and late events in monomer to amyloid fibril formation. Annu. Rev. Phys. Chem. 2011, 62, 437−463. (3) Cooper, J. A. The Role of Actin Polymerization in Cell Motility. Annu. Rev. Physiol. 1991, 53, 585−605. (4) Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 2234. (5) Cavasotto, C. N.; Bomblies, R.; Luitz, M.; Zacharias, M. In Silico Drug Discovery and Design: Theory, Methods, Challenges, and Applications; CRC Press: Boca Raton, FL, 2015; pp 313−335. (6) Kollman, P. Free energy calculations: Applications to chemical and biochemical phenomena. Chem. Rev. 1993, 93, 2395−2417. (7) Dror, R. O.; Pan, A. C.; Arlow, D. H.; Borhani, D. W.; Maragakis, P.; Shan, Y.; Xu, H.; Shaw, D. E. Pathway and mechanism of drug binding to G-protein-coupled receptors. Proc. Natl. Acad. Sci. U. S. A. 2011, 108, 13118−13123. (8) Schimmack, W.; Sapper, H.; Lohmann, W. Stacking Interactions of nucleobases: NMR-Investigations I. Selfassociation of N6,N9Dimethyladenine and N6-Dimethyl-N9-Ethyladenine. Biophys. Struct. Mech. 1975, 1, 113−120. (9) Pörschke, D.; Eggers, F. Thermodynamics and Kinetics of BaseStacking Interactions. Eur. J. Biochem. 1972, 26, 490−498. (10) Heyn, M. P.; Nicola, C. U.; Schwarz, G. Kinetics of the basestacking reaction of N6-dimethyladenosine. An ultrasonic absorption and dispersion study. J. Phys. Chem. 1977, 81, 1611−1617. (11) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins: Struct., Funct., Genet. 2006, 65, 712−725. (12) Case, D.; Babin, V.; Berryman, J.; Betz, R.; Cai, Q.; Cerutti, D.; T.E. Cheatham, I.; Darden, T.; Duke, R.; Gohlke, H.; Goetz, A.; Gusarov, S.; Homeyer, N.; Janowski, P.; Kaus, J.; Kolossváry, I.; Kovalenko, A.; Lee, T.; LeGrand, S.; Luchko, T.; Luo, R.; Madej, B.; Merz, K.; Paesani, F.; Roe, D.; Roitberg, A.; Sagui, C.; Salomon-Ferrer, R.; Seabra, G.; Simmerling, C.; Smith, W.; Swails, J.; Walker, R.; Wang, J.; Wolf, R.; Wu, X.; Kollman, P. AMBER 14; University of California: San Francisco, 2014. (13) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (14) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926−935. (15) Mahoney, M. W.; Jorgensen, W. L. A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions. J. Chem. Phys. 2000, 112, 8910− 8922. (16) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. The missing term in effective pair potentials. J. Phys. Chem. 1987, 91, 6269−6271. (17) Wu, Y.; Tepper, H. L.; Voth, G. A. Flexible simple point-charge water model with improved liquid-state properties. J. Chem. Phys. 2006, 124, 024503. (18) 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. 1977, 23, 327−341. (19) Hopkins, C. W.; Le Grand, S.; Walker, R. C.; Roitberg, A. E. Long-Time-Step Molecular Dynamics through Hydrogen Mass Repartitioning. J. Chem. Theory Comput. 2015, 11, 1864−1874. (20) Pastor, R. W.; Brooks, B. R.; Szabo, A. An Analysis of the Accuracy of Langevin and Molecular Dynamics Algorithms. Mol. Phys. 1988, 65, 1409−1419. (21) Andersen, H. C. Molecular dynamics simulations at constant pressure and/or temperature. J. Chem. Phys. 1980, 72, 2384−2393. (22) 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. 1984, 81, 3684−3690.

nucleobase hydrogen atoms by SHAKE was found to be negligible (see TIP3P 1 fs). The water models TIP4P,14 TIP5P,15 and SPC/E 4 fs yield rates between the rates found for TIP3P 4 fs and SPCFW 1 fs, with the TIP4P model performing closest to SPCFW. We conclude that also more complex four-point or five-point rigid body water models do not yield kinetic rates as close to experimental measurements as the flexible water model.

4. CONCLUSION The association of multiple N6N9-dimethyladenine molecules forms an excellent model system for a systematic investigation of oligomerization processes. In the present study, we demonstrate that current MD simulations based on a molecular mechanics force field are capable of predicting thermodynamics and kinetics of stacking oligomerization in near quantitative agreement with experiments. It has also been possible to identify the mechanistic basis of the stacking oligomerization and to confirm a random isodesmic stacking association behavior. The strategy of analyzing bulk oligomer simulations supplemented by multiple dimerization simulations can serve as a basis for tackling more complex oligomerization processes in future studies. The simulations indicate that for N6N9dimethyladenine, neither the thermodynamics nor the kinetics of the stacking of an additional base are affected by the stacking state of the partner molecules. On the molecular level, the unstacking and stacking follow distinct pathways that involve mostly motions parallel to the base planes. The good agreement with experiments, however, critically depends on the solvent model and was only obtained with the fully flexible SPCFW water model. With this model not only the kinetics but also the entropic and enthalpic contributions to binding were correctly predicted. The rigid body TIP3P and SPCE/E water models as well as more advanced rigid water models yield significantly accelerated kinetics. In the case of the rigid TIP3P and SPC/E models, the magnitudes of both calculated entropic and enthalpic contributions are underestimated. It is likely that this result is also important for the simulation of other binding and folding processes. Often, reasonable agreement in terms of free energies but not for entropic and enthalpic contributions is achieved which can especially affect the understanding of cooperativity.

■

AUTHOR INFORMATION

Corresponding Author

*Phone: +49 89 289 12335. Fax: +49 89 289 12444. E-mail: [email protected]. ORCID

Martin Zacharias: 0000-0001-5163-2663 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank the Deutsche Forschungsgemeinschaft (DFG/SFB863 project A10) for financial support. This work was also supported by the LRZ supercomputer center (Leibniz Rechenzentrum) through grant pr48po.

■

REFERENCES

(1) Sipe, J. D.; Cohen, A. S. Review: History of the Amyloid Fibril. J. Struct. Biol. 2000, 130, 88−98. 3010


Article

Journal of Chemical Theory and Computation (23) Kunsch, H. R. The Jackknife and the Bootstrap for General Stationary Observations. Ann. Statist. 1989, 17, 1217−1241. (24) The PyMOL Molecular Graphics System, version 1.7.4; Schrödinger LLC: New York, 2010. (25) Gilson, M. K.; Given, J. A.; Bush, B. L.; McCammon, J. A. The Statistical-Thermodynamic Basis for Computation of Binding Affinities: A Critical Review. Biophys. J. 1997, 72, 1047−1069. (26) Mahoney, M. W.; Jorgensen, W. L. Diffusion constant of the TIP5P model of liquid water. J. Chem. Phys. 2001, 114, 363−366. (27) Krynicki, K.; Green, C. D.; Sawyer, D. W. Pressure and temperature dependence of self-diffusion in water. Faraday Discuss. Chem. Soc. 1978, 66, 199−208. (28) Pastor, R. W.; Karplus, M. Inertial effects in butane stochastic dynamics. J. Chem. Phys. 1989, 91, 211−218.

3011


Thermodynamics and Kinetics of Nucleobase Stacking

Recommend Documents