Biomolecular Simulations under Realistic Macroscopic Salt Conditions
Biomolecular simulations are typically performed in an aqueous environment where the number of ions remains fixed for the duration of the simulation, generally with either a minimally neutralizing ion environment or a number of salt pairs intended to match the macroscopic salt concentration. In contrast, real biomolecules experience local ion environments where the salt concentration is dynamic and may differ from bulk. The degree of salt concentration variability and average deviation from the macroscopic concentration remains, as yet, unknown. Here, we describe the theory and implementation of a Monte Carlo osmostat that can be added to explicit solvent molecular dynamics or Monte Carlo simulations to sample from a semigrand canonical ensemble in which the number of salt pairs fluctuates dynamically during the simulation. The osmostat reproduce the correct equilibrium statistics for a simulation volume that can exchange ions with a large reservoir at a defined macroscopic salt concentration. To achieve useful Monte Carlo acceptance rates, the method makes use of nonequilibrium candidate Monte Carlo (NCMC) moves in which monovalent ions and water molecules are alchemically transmuted using short nonequilibrium trajectories, with a modified Metropolis-Hastings criterion ensuring correct equilibrium statistics for an \xce\x94\xce\xbc, N, p, T) ensemble. We demonstrate how typical protein (DHFR and the tyrosine kinase Src) and nucleic acid (Drew-Dickerson B-DNA dodecamer) systems exhibit salt concentration distributions that significantly differ from fixed-salt bulk simulations and display fluctuations that are on the same order of magnitude as the average.
Introduction<!>Biomolecular behavior can be sensitive to salt environments<!>Current simulation practice arbitrarily fixes microscopic salt composition<!>Simulations in the semigrand canonical ensemble can mimic real salt fluctuations<!>Nonequilibrium candidate Monte Carlo (NCMC) can achieve high acceptance rates<!>An NCMC osmostat can be used alongside thermostats and barostats<!>How do salt environments vary in realistic biomolecular simulations?<!>Outline<!>An NCMC osmostat for sampling ion fluctuations in the semigrand ensemble<!>The semigrand ensemble models salt exchange with a macroscopic salt reservoir<!>Gibbs sampling provides a modular way to sample from the semigrand ensemble<!>The chemical potential \xce\x94\xce\xbc must be calibrated to model macroscopic salt concentrations<!>Free energies for salt insertion can be efficiently computed using SAMS<!>The osmostat maintains electrostatic neutrality<!>Nonequilibrium candidate Monte Carlo is used to enhance sampling efficiency<!>We use an integrator that minimizes configuration sampling bias<!>Ionic strength influences the effective salt concentration<!>Simulation packages add different amounts of salt<!>Systems considered in the study<!>General simulation details<!>Calibration of the chemical potential<!>Optimization of the NCMC protocol<!>Quantifying the scaling behavior of the osmostat<!>Estimating the efficiency of ion configuration sampling with NCMC<!>Quantifying the salt concentration around biomolecules<!>SAMS simulations and BAR estimates accurately capture salt insertion free energies.<!>The chemical potential for a macroscopic salt concentration can be reliably determined<!>The average salt concentration is highly sensitive to chemical potential<!>Different water models have distinct chemical potentials for the same salt concentration<!>NCMC greatly enhances the sampling efficiency of salt insertion and deletion moves<!>An NCMC osmostat can rapidly equilibrate the salt concentration in biomolecular systems<!>Fluctuation magnitude grows with system size and macroscopic salt concentration<!>Salt concentrations vary significantly in typical biomolecular systems<!>Simulations containing charged biomolecules can experience salt concentrations that deviate systematically from the macroscopic concentrations<!>The accuracy of heuristic salination schemes is system dependent<!>The ionic strength exceeds the salt concentration for charged macromolecules<!>The osmostat accurately represents the local salt concentration around DNA<!>The NCMC osmostat can efficiency of ion-biomolecule interactions<!>Discussion<!>Potential applications<!>Enhancing realism in molecular simulations<!>Improving osmostat efficiency<!>Conclusion
<p>Molecular dynamics simulations have proven themselves a powerful tool for studying the structure, dynamics, and function of biomolecular systems in atomic detail. Current state-of-the-art approaches simulate a small volume around the biomolecule using explicit atomistic solvent to model the local environment. To more realistically emulate electrostatic screening effects in the local solvent environment, explicit ions are generally added, both to achieve net neutrality and to mimic the macroscopic salt concentration in the in vitro or in vivo environment being studied.</p><p>Salt concentrations and ionic composition are tightly regulated in biology4. Ion composition differs between inter/intracellular environments2, tumor microenvironments5, and organelles1 (see Figure 1, left). The local ionic concentration in the environment around real biological macromolecules, however, can significantly deviate from macroscopic concentrations. Many biomolecules possess a significant net charge, and the energetic penalty for physical systems to maintain charge separation over large distances serves to recruit more or less ions from bulk to maintain charge neutrality over macroscopic lengthscales. Yet, the number of ions within the immediate vicinity may not necessarily counter the net charge of the macromolecule, as proteins can predominantly bind to ions that have the same polarity as their net charge6. Additionally, statistical fluctuations in the total number of ions in the region around the biomolecule may result in significant variance in the local salt concentration, where relative concentration fluctuations diminish slowly with increasing simulation volume (Figure 1, middle).</p><!><p>The conformations, dynamics, function, and binding of biological macromolecules can be exquisitely sensitive to the salt concentration and composition of the local environment. The Hofmeister effect, in which ions modulate the strength of the hydrophobic effect—a major driving force in protein folding and association7,8—has been known since at least the nineteenth century9-11. Biomolecular interactions involving highly charged nucleic acids—such as DNA:protein interactions critical for DNA repair12—have been observed to show sensitivity to macroscopic salt concentrations13, as have DNA:antibiotic interactions14. In the realm of pharmaceutical design, where there is great interest in engineering small molecule ligands, salt effects are known to modulate the interactions of small molecules with proteins15 or with supramolecular hosts16.</p><!><p>In contrast to real physical systems, where the local region near the biomolecule is able to exchange ions with a macroscopic reservoir at a fixed salt concentration (Figure 1, middle), simulations of biomolecules typically fix the number of salt molecules present in the simulation volume. There is a great deal of diversity in how the fixed number of added ions is typically determined: Along with the specified macroscopic ion concentration, simulation packages may make use of the total cell volume (e.g., Gromacs17), the total solvent volume excluding the biomolecular solutes (e.g., CHARMM-GUI18), or the number of water molecules (converting the ion concentration into mole or mass fraction, as in OpenMM19,20). Some simulation packages choose to use only minimal neutralizing counterions or no counterions at all, relying on uniform background neutralizing charge to allow treatment of long-range electrostatics by particle mesh Ewald (PME) methods21,22 (such as Schrödinger's FEP+ alchemical free energy calculations23). In simulation volumes large enough to mimic the inclusion of a macroscopic salt reservoir far from the biomolecular system of interest, the environment near the biomolecule may be accurately represented, but long correlation times for well-ordered ions may still hinder equilibration of the ion environment24-26.</p><!><p>Simulations in the (semi)grand canonical ensemble, however, can—at least in principle—remedy this situation by explicitly allowing one or more components (such as ions) to fluctuate over the course of the simulation via grand canonical Monte Carlo (GCMC) moves (Figure 1, right). In grand and semigrand canonical methods, simulations are placed in thermodynamic equilibrium with a theoretical reservoir of components. The simulation can exchange molecules/particles with the reservoir, and the concentration the components in the reservoir are specified by their respective chemical potentials. Before running these simulations, one first has to determine the mapping between the concentration in the reservoir and chemical potentials, a process we refer to as calibration. Sampling over ion concentrations in explicit water via straightforward GCMC is difficult: Monte Carlo insertion/deletions have to overcome long-range effects, low acceptance rates for instantaneous Monte Carlo moves, and the concentration is sensitive to small (< kBT) variations in the chemical potential. Some efforts have circumvented these issues by using implicit solvent models6,27, cavity-biased insertions in specialized solvent models 28, and explicit solvent reorganization moves29. Osmotic ensemble Monte Carlo schemes that use fractional ions and Wang-Landau approaches have also proven themselves to be useful in simulations of simple aqueous electrolytes30,31.</p><!><p>More recently, nonequilibrium candidate Monte Carlo (NCMC) has been shown to be an effective solution to the problem of low acceptance rates when inserting or deleting particles32. In contrast to an instantaneous Monte Carlo (MC) proposal in which an inserted particle is switched instantaneously on and may clash with other solvent or solute particles, in an NCMC proposal, the particle is switched on slowly as the system is allowed to relax via some form of dynamics. NCMC uses a modified acceptance criteria that incorporates the nonequilibrium work to ensure that the resulting endpoints sample from the equilibrium distribution. With well-tuned nonequilibrium protocols, NCMC acceptance rates can be astronomically higher than their instantaneous MC counterparts32. In work simulating biomolecules at constant-pH, for example, Roux and coworkers have demonstrated how NCMC is effective at achieving high acceptance rates for NCMC proposals that also transmute an ion to/from a water molecule to maintain net charge neutrality of the system33,34.</p><p>While calibration of the effective chemical potential for the water and ion forcefields and simulation parameters at hand is nontrivial, this technical challenge can be satisfyingly addressed with existing technologies: Self-adjusted mixture sampling (SAMS)35, a form of adaptive expanded ensemble sampling36, can be used to conveniently achieve uniform sampling of all relevant salt concentrations in a single simulation, while the Bennett acceptance ratio (BAR) can optimally extract estimates of the relevant free energy differences from all NCMC proposals along with good estimates of statistical error and minimal bias37-39. Independent simulations at each salt concentration could be performed separately, with nonequilibrium switching trajectories used to estimate relative free energies between different numbers of salt pairs. However, SAMS helps more rapidly decorrelate the configurations of ions and, in principle, allows a single simulation to be used for calibration.</p><!><p>Here, we present a new approach that makes use of NCMC to insert/delete salt pairs with high acceptance probability in a manner that correctly models the statistical mechanics of exchange with a macroscopic salt reservoir. The osmostat needs to be calibrated once for the specified solvent and ion models, simulation parameters, and thermodynamic conditions (temperature, pressure, pH, etc.). Following calibration, the osmostat is used in a manner similar to a Monte Carlo barostat, attempting to modify the system composition (and hence interaction potential) at regular intervals to ensure sampling from a target probability density that models a system in equilibrium with a macroscopic salt reservoir (Figure 2). Similar to a Monte Carlo barostat19,40, the osmostat moves can be integrated alongside molecular dynamics simulations and other Monte Carlo schemes to sample from equilibrium distributions with specified thermodynamic control parameters. This composability is a general feature of Markov chain Monte Carlo moves, which provide a useful framework for designing modular algorithms for biomolecular simulation41.</p><!><p>Once we have developed and validated this tool, we use it to ask biophysical questions about the nature of salt environments around biological macromolecules: What is the average salt concentration in the simulation volume, and how does it compare to bulk? Which heuristic scheme, if any, most closely approximates the local salt concentration: macroscopic concentration times total cell volume or solvent volume, or mole fraction of water molecules? How much does the local salt concentration and ionic strength vary in "typical" biomolecular simulation conditions for different classes of biomolecular systems, such as proteins and nucleic acids? And can a Monte Carlo osmostat reduce correlation times for ions over that seen in standard MD simulations, such as the slow correlation times in ion environments around nucleic acids25? We consider some test systems that represent different classes of common biomolecular simulations: TIP3P42 (and TIP4P-Ew43) water boxes, dihydrofolate reductase (DHFR), the apo kinase Src, and the Drew-Dickerson B-DNA dodecamer25 as a typical nucleic acid.</p><!><p>This paper is organized as follows: First, we review the theory behind (semi)grand canonical ensembles that model the fluctuations experienced by small subvolumes surrounding biomolecules. Second, we describe the algorithmic design of the osmostat used to allow salt concentrations to fluctuate dynamically. Finally, we apply the osmostat to address biophysical questions of interest and discuss the nature of salt distributions and their fluctuations.</p><!><p>An osmostat is like a thermostat or barostat but allows the number of salt pairs in the simulation box to change dynamically under the control of a conjugate thermodynamic parameter—here, the chemical potential of salt. Salt pairs can be thought of as being exchanged with a macroscopic reservoir, with the free energy to add or remove salt to this reservoir described by the applied chemical potential. In principle, an osmostat could be implemented by including a number of noninteracting ("ghost") molecules in the simulation volume, turning their interactions on and off to allow the number of active salt molecules to fluctuate dynamically; alternatively, new salt molecules could be introduced or removed dynamically using reversible-jump Monte Carlo (RJMC) methods44. In either case, solvent cavity formation to accommodate ions would almost certainly require nonequilibrium protocols that employ soft-core potentials and significant tuning of these insertion/deletion protocols to achieve high acceptance rates.</p><p>To simplify implementation for the ions most commonly used in biomolecular simulations (such as NaCl or KCl), we instead choose to exchange the identities of water molecules and salt ions, where our conjugate thermodynamic parameter Δμ2·H2O–NaCl (which we will abbreviate as Δμ) will represent the difference in chemical potential between withdrawing an NaCl molecule from the reservoir while returning two H2O molecules. Because solvent cavities are not being created or destroyed—only modified slightly in size—this should provide superior phase space overlap between initial and final states.</p><p>We denote the total number of water molecules and ions as N, and define the identities of the water molecules and ions with the vector θ = (θ1, θ2, …, θN) with θi ∈ {−1, 0, +1} to denote anions (θi = −1), water (θi = 0), and cations (θi = +1), respectively (with the potential to extend this to divalent ions by adding −2, +2). This choice of labeling allows us to define the total number of Na+ ions as
(1)NNa+(θ)=∑iNδ(+1,θi),
the total number of Cl− ions as
(2)NCl−(θ)=∑iNδ(−1,θi),
and the number of water molecules as
(3)NH20(θ)=∑iNδ(0,θi),
where δ(x, y) denotes the Kronecker delta, which is unity when x = y and zero otherwise, and sums run from i to N. Note that the total number of waters and ions, N ≡ NNa+(θ) + NCl−(θ) + NH20(θ), is fixed, and does not depend on θ. We define the total charge number of the biomolecules, excluding counterions, as z.</p><p>When z ≠ 0, counterions will be added to ensure that the total charge of the simulation system is zero. The system can be neutralized by any of choice of θ that satisfies n(θ) = −z, where the total charge due to ions is given by
(4)n(θ)=∑iNθi.
As neutralizing the system will lead to unequal numbers of Na+ and Cl−, we define the amount of salt as the number of neutral pairs,
(5)NNaCl(θ)≡min{NNa+(θ),NCl−(θ)}.</p><!><p>When our osmostat is combined with a scheme that samples the isothermal-isobaric (N, p, T) ensemble, we formally sample the semigrand-isothermal-isobaric ensemble (Δμ,N, p, T). The associated equilibrium probability density is given by
(6)π(x,θ;Δμ,N,p,T)=1Ξ(Δμ,N,p,T)δ(n(θ),−z)e−β[U(x,θ)+pV(x)+ΔμNNaCl(θ)],
where the Kronecker delta δ(n(θ),−z) imposes net charge neutrality, β ≡ 1/kBT is the inverse temperature, and Ξ(Δμ, N, p, T) is the normalizing constant, given by
(7)Ξ(Δμ,N,p,T)=∑θδ(n(θ),−z)∫dxe−β[U(x,θ)+pV(x)+ΔμNNaCl(θ)],
where the outer sum is over all identity vectors and the integral is over all configuration space. For brevity, the dependence of π and Ξ on z will be omitted. It is also possible to express the probability density of the system as a function of the total number of cations and anions, rather than as function of θ. This can be achieved by summing π(x, θ; Δμ, N, p, T) over all identity vectors that preserve the neutral charge of the system and NNaCl(θ) at some constant value NNaCl′:
(8)π(x,NNaCl′;Δμ,N,p,T)=∑θδ(NNaCl(θ),NNaCl′)π(x,θ;Δμ,N,p,T)∝N!NNa+′!NCl−′!NH20′!e−β[U(x;NNaCl′)+pV(x)+ΔμNNaCl′],
where U(x;NNaCl′) is the potential energy for a system with fixed particle identities that contains NNaCl′ salt pairs. The factorial prefactors account for the degeneracy number of identity vectors θ that satisfy the constraintsNNaCl(θ)=NNaCl′ and n(θ) + z = 0.</p><!><p>Gibbs sampling framework can be used to create a modular simulation scheme in which the osmostat updates molecular identities infrequently while some MCMC scheme (such as Metropolis Monte Carlo or Metropolized molecular dynamics) updates particle positions using fixed particle identities:
(9)x~π(x∣θ,N,p,T)∝e−β[U(x,θ)+pV(x)]
(10)θ~π(θ∣x,Δμ,N,p,T)∝e−β[U(x,θ)+ΔμNNaCl(θ)]
By embedding this approach in a Gibbs sampling framework, it allows the osmostat to readily be combined with other sampling schemes that make use of a Gibbs sampling framework such as replica exchange and expanded ensemble simulations45.</p><p>Instead of instantaneous MC switching to propose changes in the chemical identities θ at fixed configuration x, nonequilibrium candidate Monte Carlo (NCMC) is used to propose updates of chemical identities and positions simultaneously as sufficiently long switching trajectories can sampling efficiencies that are orders of magnitude larger than instantaneous proposals32:
(11)x~π(x∣θ,N,p,T)∝e−β[U(x,θ)+pV(x)]
(12)x,θ~π(x,θ∣N,p,T,Δμ)∝e−β[U(x,θ)+pV(x)+ΔμNNaCl(θ)]
NCMC uses a modified Metropolis-Hastings acceptance protocol in which the appropriate total work for switching is accumulated during the nonequilibrium proposal and used in the acceptance criterion.</p><!><p>Simulating a system that is in chemical equilibrium with an infinitely large saline reservoir at a specified salt concentration first requires the calibration of the chemical potential Δμ. There are multiple ways that one could compute the necessary chemical potential. For instance, one could approximate the reservoir with a sufficiently large box of water, and narrow-in on the chemical potential that produces the desired salt concentration using stochastic approximation or the density control method recommended by Speidal et al. 46. However, this requires carrying out separate calibration calculations for each desired macroscopic concentration. Instead, we aim to construct a simple calibration procedure by computing the free energies to insert salt pairs into a sufficiently large box of water. We then use these free energies to analytically compute macroscopic salt concentrations over a wide range of chemical potentials, providing a relationship that can be numerically inverted. This procedure need be done only once for a specified ion and water model, though it may need to be repeated if the method used to compute long-range electrostatic interactions is modified.</p><p>Our calibration method is similar in principle to that of Benavides et al. 47, who estimated the chemical potential of NaCl by calculating the free energy to insert NaCl to over a range of concentrations. However, unlike47—where the goal was to estimate the solubility of NaCl—our interest in estimating the chemical potential lies solely in its ability to determine the chemical potential of the osmostat saline reservoir corresponding to the desired macroscopic salt concentration in order to induce the appropriate salt distribution on microscopic simulation systems.</p><p>Our approach to calibration computes the free energies to add NNaCl ∈ {1, 2, … ,NNaCl, max} salt pairs to n initially pure box of water. We limit our free energies calculations to insert NaCl up to some maximum NNaCl, max ≪ N for practical convenience. No constraint is placed on the amount of salt that can be added in osmostat simulations—instead, the value of NNaCl, max impacts the accuracy with which the osmostat can reproduce high macroscopic salt concentrations. We define the absolute dimensionless free energy of a system with NNaCl salt pairs at pressure p and temperature T as f(NNaCl, N, p, T),
(13)f(NNaCl,N,p,T)≡−ln(Z(NNaCl,N,p,T)Z(0,N,p,T)),
where the partition function Z(NNaCl′,N,p,T) is given by
(14)Z(NNaCl′,N,p,T)=∑θδ(NNaCl(θ),NNaCl′)∫dxe−β[U(x,θ)+pV(x)]
(15)=N!NNa+′!NCl−′!NH20′!∫dxe−β[U(x;NNaCl′)+pV(x)],
where the number of water molecules NH2O′=N′−2⋅NNaCl′. For convenience, we define relative free energies as
(16)Δf(NNaCl,N,p,T)≡f(NNaCl+1,N,p,T)−f(NNaCl,N,p,T).
For simplicity, we shall use f(NNaCl) and Δf(NNaCl) as abbreviations to equations 13 and 16, respectively. The free energies f(NNaCl) can then be used to calculate the average number of salt pairs as a function of the chemical potential Δμ,
(17)〈NNaCl〉Δμ,N,p,T=Ξ(Δμ,N,p,T)−1∑NNaCl=0NNaCl,maxNNaCle−f(NNaCl)+βΔμNNaCl
where the semigrand partition function Ξ(Δμ, N, p, T) (the same one from equation 7) can be compactly written as
(18)Ξ(Δμ,N,p,T)=∑NNaCl=0NNaCl,maxe−f(NNaCl)+βΔμNNaCl
Knowledge of f(NNaCl) will also provide a convenient estimate of the macroscopic salt concentration. We define the macroscopic salt concentration as the mean salt concentration of a system in the thermodynamic limit, and derive in Appendix 2 the following expression for the macroscopic concentration that is amenable to computational analysis:
(19)〈c^〉Δμ,N,p,T=∑NNaCl=0NNaCl,maxNNaCle−f(NNaCl)+βΔμNNaCl∑NNaCl=0NNaCl,max〈V〉NNaCl,N,p,Te−f(NNaCl)+βΔμNNaCl,
where 〈V〉NNaCl,N,p,T is the average volume for a fixed NNaCl. The macroscopic concentration 〈c^〉Δμ,N,p,T is a monotonic function of the chemical potential Δμ. Therefore—provided one has estimates of f(NNaCl) and 〈V〉NNaCl,N,p,T—the value of the chemical potential Δμ(c) that yields a desired macroscopic concentration 〈c^〉Δμ,N,p,T can be obtained by numerically inverting equation 19.</p><!><p>One could estimate the free energies f(NNaCl) NNaCl ∈ {0, 1, … , NNaCl, max} using a NNaCl, max − 1 equilibrium calculations of the relative free energies Δf(NNaCl) or the recently developed grand canonical integration technique48,49. As the latter requires a priori knowledge of the approximate scaling of the chemical potential with the concentration, we instead opt to use the recently proposed self-adjusted mixture sampling (SAMS)35 method to facilitate the calculation of the free energies from a single simulation. SAMS is a development on the method of expanded ensembles36 (sometimes known as serial tempering50) and generalized Wang-Landau algorithms51,52. It is a stochastic approximation scheme that produces unbiased estimates of the free energies (unlike Wang-Landau) that—in the asymptotic limit—have the lowest variance out of all other stochastic approximation recursion schemes35. It can be used to sample over a discrete state space and simultaneously estimate the relative log-normalizing constant for each state. For our calibration simulations, the discrete states correspond to the number of salt pairs in the systems NNaCl ∈ {0, 1, … , NNaCl, max} and the log-normalizing constant are the desired free energies f(NNaCl). By dynamically altering a series of biasing potentials, one for each state, the SAMS algorithm asymptotically samples the discrete states according to user specified target weights35. When the target weights are uniform over the state space—as we choose herein to ensure the uncertainties in the estimated free energies are approximately equal—the biasing potentials are themselves estimates of the free energies f(NNaCl). Thus, SAMS can, in principle, calculate all f(NNaCl) in a single simulation more efficiently and conveniently than numerous independent equilibrium free energy calculations.</p><p>As we describe below, our osmostat employs NCMC, which allows us to calculate the salt-insertion free energies by processing all of the NCMC protocol work values in the SAMS simulations with BAR, even from the attempts that are rejected. BAR requires samples of forward and reverse work samples of salt insertion and deletion attempts to compute Δf(NNaCl) and its statistical uncertainty for NNaCl ∈ {0, 1, … , NNaCl, max}37-39. These relative free energies can then be summed to estimate f(NNaCl) and corresponding statistical uncertainties. Our calibration simulations therefore exploit the sampling efficiency of SAMS and the estimation efficiency of BAR.</p><p>In general, the chemical potential Δμ will need to be recalibrated if the practitioner changes temperature, pressure, water or ion forcefield models, nonbonded treatment, or anything that will affect f(NNaCl) or 〈V〉NNaCl,N,p,T. A sufficiently large water box must be used when calculating f(NNaCl) to reach a regime in which f(NNaCl) is insensitive to changes in simulation size; as we will show, our calibration simulations achieve this size insensitivity even for modest water boxes of a few thousand molecules.</p><!><p>To use PME21, a popular choice for accurate long-range electrostatics, charge neutrality of the entire system needs to be maintained to avoid the artifacts induced by application of a uniform background neutralizing charge22. Even if an alternative long-range electrostatics treatment is employed (e.g. reaction field electrostatics or other non-Ewald methods53), there is, in general, approximate equality between the total number of negative charges and positive charges in biological microenvironments as they approach macroscopic lengthscales (see Figure 1
left). From a purely theoretical perspective, the existence of a thermodynamic limit a system with a net charge depends on the particular details of the system54. For these reasons, we ensure that our proposals always maintain charge neutrality by inserting or deleting a neutral Na+ and Cl− pair.</p><p>We insert and delete a salt pair by converting Na+ and Cl− ions to two water molecules (see Figure 3). These moves convert the nonbonded forcefield parameters (partial changes q, Lennard-Jones radii σ, and Lennard-Jones potential well-depths ϵ) of the water and ion parameters. The Na+ and Cl− ions are given the same topology, geometry, and number of atoms as the water model used for the simulation. Irrespective of the choice of water model, the nonbonded ion parameters are placed on the water oxygen atom, and the hydrogen atoms or additional charge sites (such as in TIP4P) have their nonbonded interactions switched off. The manner in which salt and water are transmuted to one another are is described in Appendix 3. The mass of the ions is set as the same as water, which has no impact on the equilibrium configuration probability density, though it may distrupt the kinetics (which are not of interest here).</p><!><p>A benefit of exchanging ion and water nonbonded forcefield parameters is that this procedure avoids the need to create new cavities in solvent, a difficulty that significantly complicates particle creation and destruction techniques. Nevertheless, instantaneous Monte Carlo attempts to interconvert salt and water will be overwhelmingly rejected as it is highly unlikely that the dipoles of the molecules that surround a transmuted ion—usually solvent—will be orientated in a manner that favorably solvates the new charge. This effect is compounded by the long-range nature of Coulombic interactions. The acceptance probability for salt insertion and deletion would improve drastically if the dipoles and locations of the solvent could be redistributed during an MCMC attempt. Previously, Shelly and Patey developed a configuration bias Monte Carlo technique for the insertion and deletion of ions in grand canonical Monte Carlo29. Their method reorients dipoles in a shell surrounding the inserted or deleted ion, which improved the sampling efficiency by over two orders of magnitude29.</p><p>Here, we use nonequilibrium candidate Monte Carlo (NCMC) 32, a technique that is closely related to sequential Monte Carlo and annealed importance sampling55,56, to automatically relax systems around inserted or deleted ions, thereby boosting acceptance rates and sampling efficiencies to values far higher than reported elsewhere.</p><p>In NCMC, a Monte Carlo attempt is divided into a nonequilibrium protocol that drives the system through many intermediate states. Candidate configurations are generated by driving a chosen set of variables (thermodynamic or configurational) through these intermediate states whilst allowing unperturbed degrees of freedom to relax via dynamical propagation in response to the driving protocol. The total amount of work that is accumulated between interleaved steps of perturbation (of the variables of interest) and propagation (of the unperturbed degrees of freedom) is used to accept or reject the candidate configuration. Good NCMC acceptance rates can be achieved for a reasonable choice of nonequilibrium protocol; often, a parametric protocol is specified and the total protocol length (or NCMC switching time) is tuned to be long enough to ensure a system is sufficiently relaxed with respect to the completed perturbation but short enough to be efficient.</p><p>In our NCMC osmostat, the nonbonded parameters of the ions and water molecules being exchanged are linearly interpolated into a series of equally spaced alchemical states. Each perturbation step along the alchemical path was followed by a fixed number of time-steps of Langevin dynamics where the configurations of the whole system were integrated (see Figure 3). A full description of our Monte Carlo and NCMC procedure is provided in Appendix 3. Here, NCMC propagator uses the same Langevin integrator as used in equilibrium sampling to ensure there was no significant mismatch between the sampled densities. Our particular choice of Langevin integrator (described below) was used to avoid the long correlation times that results from fully Metropolized molecular dynamics integrators and to mitigate the configuration sampling bias that is incurred by unmetropolized finite time-step integrators.</p><!><p>Care must be taken to ensure that the total work is properly accumulated in NCMC, as incorrect accumulation of work or the use of alternative definitions will lead to erroneous computation of the acceptance probability and simulation results. For time reversible MCMC integrators, such as with generalized Hamiltonian Monte Carlo (GHMC), the total work is the protocol work: the sum of the instantaneous potential energy changes that result from each perturbation during the driving process57. If the system is relaxed in-between perturbations using propagators that do not leave the target distribution invariant, such as unmetropolized Langevin integrators, NCMC can drive systems to undesirable nonequilibrium steady states, whose statistics may differ from equilibrium. On top of the work that is already performed by the driving protocol, propagators that do not satisfy microscopic reversibility can also be considered to perform work on a system57. This work, known as the shadow work, must either be minimized or eliminated (i.e., via Metropolizing the dynamics) for NCMC to sample very close to, or exactly, from the target probability density.</p><p>The issue of shadow work accumulation is not limited to propagators in NCMC. Indeed, all finite time-step molecular dynamics integrators incur a discretization error that results in biased sampling when used without metropolization. While configuration sampling errors do not occur with GHMC, the correct acceptance criterion requires that the momenta of all particles are reversed upon rejection (or acceptance) of a proposal. The reversal of momenta results in a simulation 'retracing its steps', thereby significantly increasing correlation times and decreasing sampling efficiencies. Hamiltonian Monte Carlo sampling can suffer from even longer correlation times, as momenta are randomized for each trial, irrespective of whether the previous move was accepted or not. This problem can be mitigated by using GHMC reduced momentum flipping schemes that still rigorously sample from the target distribution58-60. Correlation times are minimized by GHMC schemes that do not reverse momenta at all, although this incurs sampling bias61.</p><p>Recently, Leimkuhler and Matthews have proposed an unmetropolized Langevin dynamics technique that incurs minimal configuration sampling bias62. The minimal error is achieved using a particular numerical scheme to update the positions and momenta at each time-step. Denoting half time-step velocity updates as V, half time-step position updates as R, and the addition of an Ornstein-Uhlenbeck process as O (the Brownian motion "kick"), the symmetry in the VRORV splitting scheme leads to a particularly favorable cancellation of configuration sampling error. Leimkuhler and Matthews also found that than VRORV exhibited the lowest error on configuration dependent quantities, such as the potential energy, in biomolecular simulations compared to other symmetric splittings. As Langevin dynamics with VRORV splitting samples very closely to the true configuration Hamiltonian, we expect its neglect within NCMC moves designed to sample configurational properties to induce very little error in sampled configurational densities. For this reason, we used the protocol work to accept or reject proposals from NCMC in our osmostat.</p><!><p>We are interested in quantifying the variation of the instantaneous salt concentration c in our osmostated biomolecular simulations, where
(20)c(x,θ)=1V(x)NNaCl(θ).
Although the salt concentration of the saline reservoir, i.e. the macroscopic concentration, is known precisely and controlled by the user, the presence of a biomolecule in a simulation, along with any neutralizing counterions, may lead to significant differences in the mean salt concentration in the simulation volume from the macroscopic salt concentration. In contrast, the mean salt concentration in an initially pure box of water should match the macroscopic salt concentration of the reservoir if the chemical potential used in the osmostat is accurately calibrated.</p><p>The Debye-Hückel theory of electrolytes provided an early, analytical treatment of dilute ionic solutions using continuum electrostatics. In Debye-Hückel theory, the ionic strength I of a system, which for our simulations is
(21)I(x,θ)=121V(x)(z2+∑i=1Nθi2),
is used to predict how the effective concentrations, or activities, of ions are affected by the presence of electrolytes in the solution. The key insight of Debye-Hückel theory is that—because of electrostatic screening—the ionic strength tempers the activity of ions, such that increasing the ionic strength of a solution lowers the effective concentration of electrolytes. Although Debye-Hückel theory is too simplistic to be used to accurately predict the salt concentration in biomolecular simulations, the ionic strength may still provide insight into the salt concentrations that we will observe in our osmostated simulations. Thus, we will investigate the variation of the ionic strength as well as the salt concentration. As a large charge number of the biomolecule z will dominate I for small simulation volumes, we will also consider the variation of ionic strength of the solvent only, i.e., by neglecting z2 in equation 21.</p><!><p>There is diversity in the way that current practitioners of all-atom biomolecular simulations add salt (salinate) to systems during the preparation stages of simulations. While it is common that only neutralizing counterions are added, a number of workflows elect not to add counterions at all23. Salt pairs may be added, or not added at all, and when they are added, simulation packages use differing definitions of salt concentration, such that each package can add different numbers of salt pairs to the same system even if the desired salt concentration is the same. All packages ignore the presence of neutralizing counterions when adding salt. In this study, we are concerned with quantifying the accuracy of some of the most popular salination techniques.</p><p>Given a target salt concentration of ct, a popular method to add salt—exemplified by the Gromacs package17—uses the initial volume of the system V(x0) to count the required number of pairs. We determine the number of salt pairs that would be added by this strategy as
(22)N^NaClV=⌊V(x0)ct⌋,
where ⌊y⌋ denotes the floored value of y. We are interested in assessing the accuracy of the corresponding concentration of salt c^V(x)=N^NaClV∕V(x). Preparation tools such as CHARMM-GUI18 add salt based on the initial volume of the solvent
V(x0,H2O), which we reproduce with
(23)N^NaClS=⌊V(x0,H2O)ct⌋,
to estimate the corresponding concentration c^S(x)=N^NaClS∕V(x) that would occur for all later configurations. Estimates that use strategies similar to equations 22 and 23 are sensitive to initial volume of the system; if salt is added before the volume is sufficiently equilibrated, the salt concentration during the simulation can deviate significantly from the target concentration. In contrast, packages such as OpenMM19,20, use the ratio of salt pairs to water molecules in bulk solvent to add
(24)N^NaClR=⌊NH2Oc^H2Oct⌋,
salt pairs, where c^H2O is concentration of bulk water, for which 55.4 M is used by OpenMM. The corresponding salt concentration c^R(x)=N^NaClR∕V(x), as well as c^V(x) and c^S(x) will be compared to the concentration of salt that results from the application of our osmostat to help inform future simulation strategies.</p><!><p>The primary aims of this study are to quantify and understand how the concentration of salt and ionic strength vary around typical biomolecules, to assess the accuracy of methods that insert salt in typical simulation strategies, and to ascertain whether an NCMC osmostat can decorrelate biomolecule:ion interactions faster than fixed-salt dynamics. To meet these aims, we considered four biological systems that are representative of those that are commonly simulated with molecular dynamics: pure water, dihydrofolate reductase (DHFR), the apo kinase Src, and the Drew-Dickerson B-DNA dodecamer palindromic sequence. All systems were taken from the OpenMMTools [0.11.1] set of test systems63, such that each system has a different provenance.</p><p>Dihydrofolate reductase (DHFR) is a small, globular enzyme that has frequently been used as a model system in molecular simulations. The DHFR structure used here was taken from the joint Amber-CHARMM (JAC) benchmark (obtained from the Amber 14 benchmark archive64). The protein structure was stripped of hydrogen atoms, and using tleap65, was re-protonated at pH 7 and solvated in an orthorhombic box of TIP3P waters that had a clearance of at least 10 Å. The Amber 14SB forcefield from the AmberTools 16 package was used for the protein65. As an initial relaxation of the system, the solvated system was minimized and propagated for 3 ps with Langevin dynamics.</p><p>The tyrosine kinase Src, a member of the non-receptor tyrosine kinase family, was selected for this study as an example of a prototypical drug target. The apo Src structure was taken from the OpenMMtools testsystems data set and resolvated with TIP3P in an orthorhombic box that was at least 10 Å away from the protein. As part of the preparation, the energy of system was minimized and subsequently relaxed using 3 ps of Langevin dynamics to remove any bad contacts. Further equilibration was performed as detailed below. The original system was not suitable for simulation with the osmostat as fixed neutralizing counterions were present in the system. The OpenMMtools structure was downloaded from the Protein Data Bank, identification code 1YI6, and prepared using PDBFixer66 and protonated at pH 7. The small molecule in the binding site was also removed during the preparation. The Amber 14SB forcefield from the AmberTools 16 package was used for the simulations65.</p><p>The Drew-Dickerson dodecamer (CGCGAATTGCGC) is a classic model DNA system. The B-DNA structure of the Drew-Dickerson dodecamer was downloaded from the Protein Data Bank (identification code 4C64). The structure was stripped of ions and solvated in a box of TIP3P water to ensure at least 9 Å of clearance around the DNA. To test the effect of the amount of solvent on the distribution of salt and ions, the structure was also solvated in a box of TIP3P water that had a clearance of at least 16 Å around the DNA. As with the apo kinase Src, the system was energy minimized and subsequently relaxed using 3 ps of Langevin dynamics. As described below, further equilibration was also performed. The Amber OL15 forcefield from the AmberTools 16 package was used for the DNA67.</p><!><p>Simulations were performed with OpenMM [7.1.0]20. The osmostat was implemented within the open-source package SaltSwap [0.5.2] that was written for the purpose of this publication. Simulations utilized either TIP3P42 or TIP4P-Ew43 water models, and Joung and Cheatham parameters were used for Na+ and Cl− ions68. Unless otherwise stated, the amount of salt in a simulation was initialized by salinating the system according to equation 24 with the macroscopic concentration as the target concentration ct.</p><p>For all simulations, long-range electrostatic interactions were treated with particle mesh Ewald (PME), with both direct-space PME and Lennard-Jones potentials making use of a 10 Å cutoff; the Lennard-Jones potential was switched to zero at the cutoff over a switch width of 1.5 Å to ensure continuity of potential and forces. PME used a relative error tolerance of 10−4 at the cutoff to automatically select the α smoothing parameter, and the default algorithm in OpenMM was used to select Fourier grid spacing (which selected a grid spacing of ~0.8 Å in each dimension). All bonds to hydrogen were constrained to a within a fractional error of 1 × 10−8 of the bond distances using CCMA69,70, and waters were rigidly constrained with SETTLE71. OpenMM's long-range analytical dispersion correction was used to avoid pressure artifacts from truncation of the Lennard-Jones potential. Simulations were run at 300 K with a Monte Carlo barostat with 1 atm external pressure and Monte Carlo update interval of 25 steps. Equilibrium and NCMC dynamics were propagated using high-quality Langevin integrators taken from the OpenMMTools [0.11.1] package, with a 2 fs timestep and collision rate of 1 ps−1. Integrators used deterministic forces and OpenMM's mixed single and double precision implementation. In addition to the dynamics used to prepare the systems, every simulation was briefly thermalized using 4 ps of dynamics. Where stated, additional simulation data was discarded from the start of simulations using the automatic procedure in the pymbar timeseries module as detailed in72. As described above, positions and velocities were updated using the VRORV splitting scheme (also known as BAOAB) to mitigate the configuration space error in equilibrium sampling and NCMC proposals that result from unmetropolized Langevin dynamics62</p><p>The insertion or deletion of salt was attempted every 4 ps using the procedure described in Appendix 3. All ions used the same number of atoms, topology, and geometry as the water model used in the simulation. As illustrated in Figure 3, the "insertion" of an ion was achieved by switching the nonbonded parameters of the water oxygen atom to either Na+ or Cl− and by simultaneously switching the nonbonded parameters of the water hydrogen atoms (along with any extra charge sites) to zero—the "deletion" of an ion involved the reverse procedure. With the exception of the simulations where the NCMC protocol was optimized, the NCMC protocol was 20 ps long, and consisted of 1000 perturbation steps, where each perturbation followed by 10 steps of Langevin integration with a 2 fs timestep. The pseudo-code for the entire NCMC osmostat, including how it is combined with molecular dynamics can also be found in Appendix 3. Unless otherwise stated, the NCMC protocol length is not accounted for in the reported lengths of the simulations.</p><p>The simulations were analyzed with open source scripts that used a combination of numpy 1.13.173, scipy 0.19.174, pymbar 3.0.175, MDTraj 1.8.076, VMD 1.9.477 (see Code and data availability); the saltswap conda package provided automatically installs the dependencies needed to run the simulation scripts. Plots and figures were produced using Matplotlib 2.0.278 and Inkscape 0.91.</p><!><p>The chemical potential was calibrated in cubic boxes of TIP3P water and TIP4P-Ew water. Both boxes initially had edge lengths of 30 Å with water molecules at roughly the same density as bulk water; the box of TIP3P water contained 887 molecules and the box of TIP4P-Ew water contained 886 molecules. Ten 80 ns SAMS simulations were performed on each box, and were targeted to sample uniformly over salt pairs NNaCl(θ) ∈ {0,1,…, 20}. The insertion or deletion of salt was attempted every 4 ps. Half of the simulations were initialized with 0 salt pairs, whereas the other half were initialized with 20 salt pairs. The maximum number of salt pairs NNaCl,max was chosen to be 20 in these calibration simulations because the corresponding salt concentration (roughly 1.2 M) is beyond the concentrations in biological microenvironments that are typically considered. (Note that the maximum amount of 20 salt pairs applies only to these calibration simulations—the osmostat simulations with solutes have no such maximum number of salt pair limitation.) The volumes of the boxes at each salt occupancy were recorded during the SAMS simulations in order to estimate 〈V〉NNaCl,N,p,T (henceforth abbreviated as 〈V〉NNaCl). The SAMS simulation procedure automatically provides on-line estimates of the free energies f(NNaCl), which, along 〈V〉NNaCl, are required to calibrate the chemical potential. The protocol work from all of the NCMC insertion and deletion attempts were post-processed with BAR (using the pymbar package75) to provide additional estimates of f(NNaCl) along with statistical uncertainties.</p><p>To assess whether f(NNaCl) and 〈V〉NNaCl had been accurately calculated, larger boxes of TIP3P and TIP4P-Ew water were simulated for 32 ns at a range of chemical potentials Δμ. The mean salt concentrations from the simulations were compared to concentrations predicted using equation 19 with the estimated values for f(NNaCl) and 〈V〉NNaCl. The boxes of these validation simulations were initially 50 Å in length, and contained 4085 TIP3P and 4066 TIP4P-Ew water molecules. These simulations were initialized without any salt present in the systems.</p><!><p>We consider only two parameters in optimizing the nonequilibrium protocol used in NCMC proposals: the total number of times the potential is perturbed, T, and the number of Langevin steps that occur before and after each perturbation, K. Generally, we expect the acceptance probability to increase as the overall perturbation is broken into smaller pieces—as T increases. Increasing the number of propagation steps following each perturbation, K, also improves the acceptance probability in a manner that is dependent on the computational efficiency details of the simulation code. To quantify the trade-off between acceptance probability and compute time, we define the NCMC efficiency E(T, K) as
(25)E(T,K)=〈A〉(T,K)C(T,K),
where 〈A〉(T, K) is the average acceptance probability and C(T, K) is the average computer time per insertion/deletion attempt. All simulations were performed and timed on single Nvidia GTX-1080 GPUs. The total protocol length of an NCMC attempt is equal to T × K multiplied by the Langevin integration timestep, which is 2 fs in this case.</p><p>Simulations using various NCMC protocols lengths were performed on cubic boxes of TIP3P and TIP4P-Ew that had initial edge lengths of 30 Å. The simulations sampled configurations for a total of 32 ns (excluding the NCMC sampling) and had NCMC protocol lengths up to 40 ps for different combinations of total perturbation steps T and propagation steps K. The insertion or deletion of salt was attempted every 4 ps, such that there were a total of 8000 insertion/deletion attempts for each simulation. The efficiency of each protocol E was estimated relative the efficiency of instantaneous insertion and deletion. Shelly and Patey also used the ratio of the average acceptance probability to the compute time to estimate the efficiency of their configuration bias ion insertion scheme relative to instantaneous insertions29. In this work, no effort was made to optimize the alchemical path.</p><!><p>To investigate the sampling efficiency of our osmostat under physiological conditions, DHFR was simulated with macroscopic concentrations of 100 mM, 150 mM, and 200 mM. Each simulation was 30 ns long and there were three repeats per macroscopic concentration. Equation 24 was used to add an initial amount of salt to the simulation. The timeseries module in pymbar75 was used to estimate the autocorrelation function of salt concentration as well as the integrated autocorrelation time for each macroscopic salt concentration.</p><p>It is important to establish how the distributions of salt concentration and salt numbers scale with the number of water molecules in the system and the macroscopic concentration. To this end, we simulated different sizes of water boxes with macroscopic concentrations of 100 mM, 150 mM, and 200 mM. Each simulation was repeated three times.</p><!><p>Ponomarev et al. previously used the Drew-Dickerson DNA palindromic sequence to quantify the rate of convergence of spatial ion distributions in DNA simulations25. Three osmostated simulations and three fixed-salt simulations of the Drew-Dickerson dodecamer were performed for 60 ns with a macroscopic salt concentration of 200 mM. As the insertion or deletion of salt was attempted every 4 ps, there was a total of 15,000 attempts. The fixed salt simulations used the same ion topologies and masses as those used by the osmostat, are were added to the system using the scheme summarized by equation 24. The autocorrelation of ion:phosphate interaction occupancies were estimated from the osmostated and fixed-salt simulations using the open-source analysis scripts that accompany this manuscript.</p><!><p>Three 30 ns simulations of apo Src kinase were performed, with salt insertion or deletion attempted every 4 ps, using a macroscopic concentration of 200 mM. The amount of salt that was initially added to this system was calculated using equation 24. These simulations, as well as those of TIP3P water, DHFR, and the DNA dodecamer described above, were used to analyze the distributions of salt concentration (equation 20), ionic strength (equation 21), and the concentrations of salt that would occur for the heuristic salination schemes described in equations 22, 23, and 24.</p><p>To further understand the scaling behavior of the distributions of salt concentration with system size, and to assess the extent of finite size effects on the ion spatial distributions around DNA, additional simulations were performed on DNA. The Drew-Dickerson DNA dodecamer was resolvated in a box of TIP3P water that was at least 16 Å away from the molecule. Three repeats of 45 ns long osmostated and fixed-salt simulations were performed, with the insertion or deletion of salt was attempted every 4 ps. The salt concentration distribution was estimated, as were the Na+ and Cl− spatial distributions around the DNA.</p><!><p>In order to estimate the chemical potential Δμcorresponding to a desired macroscopic salt concentration, we must have precise estimates of both free energies to insert salt into a box of water containing NNaCl salt molecules, f(NNaCl), and the average saline box volume as a function of NNaCl, 〈V〉NNaCl, for NNaCl ∈ {0, 1, … , NNaCl, max}. Figure 4 (upper left) depicts the computed relative free energy difference for inserting an additional salt pair into a box of water molecules already containing NNaCl salt molecules for both TIP3P and TIP4P-Ew for NNaCl ∈ {0, … , 19}. The relative free energies were estimated with BAR using all nonequilibrium work values for salt pair insertion/deletion NCMC proposals, irrespective of whether the proposal attempt was accepted or not, from ten SAMS simulation. Although SAMS also provides online estimates for f(NNaCl) over this same range35, these online estimates were found to have significantly higher variance than the BAR estimates (see Figure A5.1), so we make use of BAR-derived estimates of f(NNaCl) derived from SAMS simulations throughout.</p><p>The primary accuracy of the calibration simulations lies in their ability to reproduced desired salt concentrations in bulk water. Nevertheless, it is encouraging to note that calculated free energy to insert one NaCl pair in a box of TIP3P and TIP4P-Ew are broadly in agreement with previous computational estimates and experimental measurements. As implied by equation 16, the free energy to insert the first salt pair, Δf(NNaCl = 0), can be expressed as the difference in hydration free energy between NaCl and two water molecules. Assuming the hydration free energy of TIP3P and TIP4P-Ew water to be −6.3 kcal/mol79, we estimate the hydration free energy of NaCl to be −171.73 ± 0.04 kcal/mol and −170.60 ± 0.04 kcal/mol in TIP3P and TIP4P-Ew water, respectively. Using a different treatment of long-rang electrostatics but same ion parameters as this study, Joung and Cheatham calculated the individual hydration free energies of Na+ and Cl− in TIP3P and TIP4P-Ew, which can be summed to approximate the hydration free energy of NaCl68. These hydration free energies (−178.3 kcal/mol in TIP3P −177.7 kcal/mol in TIP4P-Ew) are within 5% of our estimates. For comparison, estimates of standard NaCl hydration free energies based on experimental data are −170.4 kcal/mol80, −171.8 kcal/mol81, and −177.8 kcal/mol82.</p><!><p>The salt insertion free energies and average volumes in Figure 4
upper left provide a way to relate the chemical potential Δμ to macroscopic salt concentration 〈c^〉 via equation 19. Figure 4
upper right shows the predicted macroscopic salt concentration for a range of chemical potentials Δμ computed using equation 19. The average salt concentration in a saline box 〈c〉 should equal the predicted macroscopic concentration for sufficiently large saline boxes if the chemical potential has been properly calibrated. To verify the accuracy of the calculated values for f(NNaCl) and 〈V〉NNaCl, simulations of water boxes, that initially had no salt present, were performed using an osmostat with different fixed chemical potentials and the average salt concentrations in the simulations were estimated (Figure 4; upper right). These boxes of TIP3P and TIP4P-Ew waters contained 4085 and 4066 molecules respectively, whereas the TIP3P and TIP4P-Ew boxes used to estimate f(N) and 〈V〉NNaCl contained 887 and 886 molecules respectively. As Figure 4
upper right shows, the macroscopic concentrations 〈c^〉 predicted using equation 19 fall within the statistical error of the average concentrations 〈c〉 determined from the fixed-Δμ simulations.</p><p>Although Δμ is the thermodynamic control parameter for osmostated simulations, experimental wetlab conditions instead generally specify the macroscopic salt concentration 〈c^〉 rather than Δμ. As the relationship between Δμ and 〈c^〉 is monotonic, as illustrated by Figure 4
upper right, we can numerically invert equation 19 to enable practitioners to choose the desired macroscopic salt concentration and extract the required Δμ for the osmostat to model equilibrium with the macroscopic salt concentration 〈c^〉.</p><!><p>The macroscopic salt concentration 〈c^〉Δμ for a fixed chemical potential Δμ is a highly sensitive and nonlinear function of the chemical potential (Figure 4; upper right) for both water models. Small changes to the chemical potential, on the order of 1 kT, can alter the mean concentration by hundreds of millimolar. Correspondingly, to accurately model a given macroscopic concentration c, the function Δμ(c) must be very precisely calibrated.</p><!><p>Strikingly, both the value and shape of 〈c^〉Δμ is very sensitive to choice of water model (Figure 4; upper right). For instance, a Δμ of about 316 kT results in a mean salt concentration in TIP3P water that is approximately 500 mM, compared to approximately 200 mM in TIP4P-Ew water for the same value of Δμ. These features highlight the importance of specifically calibrating the chemical potential for each water and ion model as well as estimating f(NNaCl) and 〈V〉NNaCl to a sufficient degree of precision. Figure A5.2 shows that for TIP3P and the treatment of long-rang interactions used herein, the free energies f(NNaC1) for each NNaCl ∈ {0, 1, … , 20} need to be determined to a standard error of 4 kcal/mol to consistently determine the macroscopic concentration to an inaccuracy of at least about 80 mM for 1 mM ≤ 〈c^〉 ≤ 1000 mM. The average standard error achieved in the calibration simulations for the free energies f(NNaC1) is 0.02 kcal/mol, which determines the concentration to an inaccuracy no larger than about 1 mM.</p><!><p>We estimate that instantaneous salt insertion and deletion moves have acceptance probabilities of 3.0 × 10−51 [95% CI: 5.0 × 10−66, 9.0 × 10−51] and 1.0 × 10−46 [95 % CI: 3.0 × 10−64, 4.0 × 10−46] in TIP3P and TIP4P-Ew water respectively, implying that the implementation of an osmostat is practically impossible using such naive moves. In contrast, we found that in our longest protocol, NCMC insertion/deletion attempts achieved acceptance probabilities of about 30% in TIP3P water and approximately 15% in TIP4P-Ew water (see the lower left of Figure 4). Although the acceptance probability increases monotonically with the length of the protocol, so does the computational cost and time for each attempt. The efficiency, defined in equation 25, quantifies the trade-off between the acceptance rate and computational expense. Figure 4
lower right shows that NCMC protocols in TIP3P water that are between 15 ps and 30 ps in length are the most efficient for our procedure. For this reason, all subsequent simulations used TIP3P water and a 20 ps long NCMC protocol. In addition, it was found that 10 propagation steps (at 2 fs) between each perturbation was found to be the most computationally efficient for our simulation code SaltSwap [0.5.2] and OpenMM [7.1.0] (see Figure A5.3). Further optimization of the NCMC protocol would be required for NCMC attempts in TIP4P-Ew to achieve sampling efficiencies that are competitive with those in TIP3P water.</p><!><p>Figure 5 shows example salt concentration trajectories around DHFR as well as plots of the corresponding autocorrelation functions for three biologically plausible macroscopic salt concentrations. The autocorrelation times for the three macroscopic salt concentrations are on the order of 1 ns, implying that our osmostated simulations should be at least tens of nanoseconds long to generate sufficient uncorrelated samples of salt concentrations. Importantly, the magnitude of the instantaneous salt concentration fluctuations increases with the macroscopic salt concentration, which causes an increase in the correlation time as our osmostat implementation proposes the insertion/deletion of one salt pair a at a time. As a result, more attempts are required to explore salt concentration distributions of higher variance. This suggests that inserting or deleting multiple salt pairs in each attempt could improve the sampling efficiency of our osmostat at higher macroscopic salt concentrations, though longer NCMC insertion/deletion protocols would likely be required to achieve similar acceptance probabilities.</p><!><p>Figure 6
upper left demonstrates that for a pure box of saline and fixed macroscopic salt concentration, increasing the number of molecules in the system increases both amount of salt and the spread of the salt number distribution; in contrast, Figure 6 (upper right) reveals that the distribution of the concentration remains centered around the macroscopic concentration, but the variance decreases. Both of these trends are to be expected from statistical mechanics (see Appendix 2). The salt concentration distribution for the smallest water box (with 2094 molecules) in Figure 6 (upper right) can be seen to be highly multimodal. Each peak corresponds a particular number of salt pairs in the system; there are so few water molecules in this system that changing NNaCl, by one results in a large jump in the concentration. Figure 6 (bottom left and right) highlight that for a system with a fixed number of water molecules, the number of salt pairs increases in proportion with the macroscopic concentration.</p><!><p>Figure 7 shows the distribution of salt concentration and ionic strength for 3 typical biomolecular systems: DHFR, apo Src kinase, and the Drew-Dickerson DNA dodecamer. The distributions in a box of TIP3P are also shown for reference. The fluctuations of the salt concentration around the macromolecules are substantial: 95% of all salt concentration samples fall within a range of 90.2 mM for DHFR, 87.7 mM for Src kinase, and 135.6 mM for the DNA dodecamer system. We expect these values to be indicative of the natural variation in salt concentration in the local environments of real biomolecules.</p><!><p>The DHFR, apo Src kinase, and the Drew-Dickerson DNA dodecamer structures have net charges of −11 |e|, −6 |e| and −22 |e|, respectively. The net charge of the DNA dodecamer is a result of the phosphate group on each of the nucleotides (with each of the eleven phosphate groups carrying −1 |e| charge), whereas the net charges on DHFR and Src kinase are due to an excess of glutamate and aspartate residues over arginine, histidine, and lysine residues. Neutralizing Na+ ions were added to both systems to avoid the uniform background charge that would be applied automatically with PME electrostatics. Like the other ions in our osmostat, these counterions had transmutable identities.</p><p>Figure 7 shows that in our osmostated simulations of the macromolecules, the average salt concentration is on average less than the macroscopic salt concentration. This is particularly apparent with the DNA dodecamer, which has a mean concentration of 128.0 [121.5, 134.5] mM (where the quantity in brackets denotes the 95% confidence interval of the mean concentration). The salt concentration distribution in the DHFR and Src kinase systems are centered closer to the macroscopic concentration of 200 mM, with estimated means of 174.0 [164.4, 180.4] and 176.3 [171.6, 189.5] mM, respectively. To compute these statistical estimates and confidence intervals, no data was discarded at the start of the simulation, and approximately statistically independent concentration samples were extracted using the pymbar timeseries module75.</p><p>The larger number of water molecules in the Src kinase system is partly the reason why its mean concentration is closer to the macroscopic value than the DNA dodecamer. Bulk-like conditions anchor the sampled salt concentrations about the macroscopic concentration; the more water molecules and salt pairs there are, the smaller the effect a macromolecule has on the salt concentration relative to the whole system. Figure 8
inset highlights this phenomenon with the DNA dodecamer; the mean salt concentration moves closer to the macroscopic value when more water molecules are added to the simulation.</p><!><p>On its own, the excluded volume of the macromolecule will reduce the number of salt pairs that can occupy the simulation volume compared to bulk saline. So, as we define the salt concentration as the number of salt pairs over the total volume of the system (equation 20), one would expect there to be a lower salt concentration than the macroscopic value. The preparation schemes that are typically used to add salt in fixed-salt simulations that account for this effect use either the volume of the solvent (equation 23), or the ratio of the number of salt pairs to water molecules (equation 24). As a result, these methods are closer to the mode of the concentration distributions in the osmostated simulations than the heuristic method that uses the total volume of the system (equation 22). The volume-based methods are sensitive to how equilibrated the volume is when salt is added, and, in Figure 7, the volume at the start of the production simulation was used to estimate the amount of salt that would be added with equations 22 and 23. The salt-water ratio method (equation 24) has no such volume dependence, which is partly why it is a better predictor for the salt concentration than the others.</p><!><p>In addition to the distributions of salt concentrations, Figure 7 also shows the ionic strength of the saline buffer. While the ionic strength is used in analytical models to estimate the activities of ionic species83, the only discernible common feature of the ionic strength in our simulations is that it tends to be greater than the salt concentration, which is predominantly due to the presence of neutralizing counterions. The estimated mean ionic strength of the saline buffer in the macromolecular systems are 208.2 [198.2, 213.6] mM for DHFR, 189.0 [179.5,196.4] mM for Src kinase, and 263.4 [256.6,269.8] for the DNA dodecamer. It is important to note that the calculated ionic strength can be much larger when the contribution of the macromolecule is included: the estimated ionic strengths for the whole of the DHFR, Src kinase, and DNA systems are 551.0 [541.0, 556.4] mM, 263.6 [253.8, 270.8] mM, and 3241.6 [3227.3, 3244.7] mM respectively. These high values, particularly for the DNA system, is because the ionic strength is proportional to the square of the charged number of the ionic solute. It could be more informative to consider the macromolecule and the counterions that are bound to it as a single, aggregate macro-ion, such that the contribution to the ionic strength would be lessened83; however, as there is no clear boundary between bound and unbound ions (see Figure 8), this approach is conceptually difficult.</p><!><p>The aim of our osmostat is to replicate the local ion concentrations that would occur around biomolecules when embedded in large saline reservoirs. However, the use of periodic simulation cells and the addition of neutralizing counterions constrains length scale at which charges are screened (the Debye length) to be less than or equal to the length scale of the periodic cell. An artificial constriction of the Debye length would be finite size effect that would limit the accuracy of the salt concentrations from osmostated simulations. Figure 8 shows the total charge contained within ever increasing distances from the Drew-Dickerson DNA dodecamer for two simulation box sizes. The smallest box was constructed by solvating the DNA up to a minimum distance of 9 Å away from the DNA (4296 water molecules), whereas the larger box resulted from solvating up to a distance of 16 Å from the DNA (9276 water molecules). If the Debye length was significantly affected by the periodic cell size of the smallest simulation, there would be large discrepancies between the charge distributions around the DNA of the smallest box and the larger box. Figure 8 indicates that if such discrepancies exists, they are small, and are not found to be statistically significant in our analysis.</p><p>Shown first in Figure 7 (lower right), the osmostated simulation of the Drew-Dickerson DNA dodecamer experienced significantly lower NaCl concentrations than the applied 200 mM macroscopic NaCl concentration. This difference highlights how the local ionic environment of a solute can be strikingly different from bulk saline. Increasing the amount of water in the simulation diminishes the relative effect that DNA has on perturbing the salt concentration distribution of the whole system. Figure 8 (inset), shows that increasing the number of water molecules in the system from 4296 to 9276 molecules partially masks the local salt concentration around the DNA, such that the total salt concentration over the whole system is closer to the macroscopic concentration of 200 mM.</p><!><p>To compare the computational efficiency of NCMC ion sampling to that of fixed-salt MD simulations, the autocorrelation functions of cation-phosphate interactions were estimated from the DNA dodecamer simulations. Cation-phosphate interactions were recorded as every time a cation was within 5 Å of the phosphorous atoms in adenine nucleotides. This cutoff was chosen following the DNA convergence analysis of Ponomarev et al.25. The autocorrelation function of these interactions measures the probability that a cation that is initially within the distance cutoff will also be present after a given amount of time. As our osmostat uses NCMC to add and remove ions, one would expect the osmostat interaction autocorrelation function to decay significantly faster than that from the fixed salt simulations when only considering the molecular dynamics—Figure 9 shows that this is indeed the case.</p><p>When the simulation time from NCMC is not considered, the phosphate-cation interaction autocorrelation function from the osmostat simulations decays significantly faster than the fixed salt simulations (Figure 9). The corresponding integrated autocorrelation times for osmostated simulations and fixed-salt simulations are 0.11 [0.09, 0.13] ns and 0.29 [0.23, 0.36] ns respectively. As each accepted NCMC move has propagated the configurations of the whole system, the faster decorrelation of DNA-ion interactions could be a result of these extra propagation steps, as opposed to the fact that ions are being inserted and deleted. As described in the methods, a salt insertion or deletion attempt occurs every 4 ps, and an NCMC attempt involves 20 ps of dynamics. The average acceptance probability in the DNA simulations was calculated to be 11.9 [11.7, 12.2] %. Therefore, the osmostated simulations propagate the system 1.6[≈(0.119×20 ps + 4 ps)/4 ps] times as much dynamics than fixed salt simulations. Multiplying the osmostated integrated autocorrelation time by this factor results in a value that remains significantly less than the integrated autocorrelation time from the fixed salt simulations. Figure 9
right shows the osmostated autocorrelation function when the timescale has been multiplied by the effective NCMC sampling factor (1.6). Despite the application of this factor, the fixed-salt autocorrelation function can be seen to decay significantly slower than the stretched osmostated autocorrelation function. Thus, the increased sampling efficiency observed in the osmostated simulations cannot be explained by the extra dynamics sampled in the NCMC simulations. This implies that the random insertion and deletion, not the NCMC that was used to enhance the move efficiency, is responsible for the rapid decorrelation of ion interactions observed in the DNA osmostated simulations.</p><p>The total number of NCMC timesteps (including from rejected moves) can be used to account for the additional computational burden of the NCMC osmostat in the phosphate-cation autocorrelation times. There is an additional 20 ps of dynamics for every insertion/deletion attempt, irrespective of whether the proposal was accepted or not. As each attempted is preceded by 4 ps of equilibrium dynamics, our osmostated simulations have 6 (= (20 ps + 4 ps)/4 ps) times as timestep evaluations than the fixed-salt simulations. Multiplying the mean integrated autocorrelation time from the osmostat simulations by this factor yields an effective autocorrelation of 0.65 [0.55, 0.75] ns. Although this estimate now exceeds the upper confidence interval of the fixed-salt integrated autocorrelation time (0.29 [0.23, 0.36] ns), there is only approximately 0.1 ns difference between the lower and upper confidence intervals. Figure 9 also shows the osmostat phosphate-ion autocorrelation function when the all the NCMC propagation steps (including rejected moves) are accounted for. One can see that for below ~1 ns, the 95% confidence intervals of the autocorrelation functions overlap with those of fixed-salt autocorrelation function. These results imply the dynamic NaCl sampling achieved by our osmostat has a similar cost effectiveness—with regards to ion sampling—than fixed-salt simulations, with the additional benefit of sampling realistic salt concentrations.</p><!><p>In this work, we have implemented an osmostat that dynamically samples the NaCl concentration in biomolecular simulations. The osmostat couples a simulation cell to a saline reservoir at a fixed macroscopic concentration and allows the salt concentration in the simulation to fluctuate about its equilibrium value. We have applied our osmostat to simulations of dihydrofolate reductase (DHFR), apo Src kinase, and the Drew-Dickerson B-DNA dodecamer (CGCGAATTGCGC), and found that the mean salt concentration can differ significantly from the amount salt added by common molecular dynamics methodologies. In addition, we found that the salt concentration fluctuations were large, being of the same order of magnitude as the mean. These results show that the ionic composition around biomolecules can be highly variable and system dependent.</p><p>The insertion and deletion of salt was greatly enhanced by nonequilibrium candidate Monte Carlo (NCMC), to the extent that the protocol used in our simulations was approximately 5 × 1046 times more efficient than instantaneous attempts in TIP3P water. The Drew-Dickerson B-DNA dodecamer is a palindromic sequence that facilitated a study of the convergence of ion distributions around the DNA. We found that, despite the additional computational expense of the NCMC osmostat, the sampling and computational efficiency of DNA:ion interactions remained comparable to fixed-salt simulations. However, it is important to note that made no effort to optimize the NCMC protocols beyond selecting an appropriate total switching time for NCMC moves—it is possible that further optimization of these protocols using recent techniques based on mapping geodesics in the thermodynamic metric tensor space84-88 can lead to increased efficiency.</p><!><p>While the dependence of enzyme-substrate activity on ionic strength is well documented, the impact of salt concentration on protein-ligand binding affinity is much less clear. Recently, Papaneophytou et al. performed a systematic analysis on the effect of buffer conditions on the in vitro affinity of three complexes15, finding salt concentration dependence to be system dependent and largest for complexes that formed hydrophilic interactions. Our osmostat provides the opportunity to rigorously study the impact of salt concentration on protein-ligand binding affinities in silico. We are interested to know if similar trends to what Papaneophytou et al. observed can be reproduced in all-atom binding free energies calculations, and whether binding free energy estimates differ significantly between simulations carried out with and without an osmostat. Free energy calculations on complexes whose association is sensitive to the concentration of salt are likely to be most affected by the osmostat, given the large fluctuations of concentration and the deviation from the fixed-salt values that occurred in our simulations (see Figure 7). The combination of self adjusted mixture sample (SAMS) and Bennett acceptance ratio (BAR) that we used to calibrate the chemical potential can also be used to estimate the difference between traditional and osmostated free energy calculations. If significant differences between binding free energy calculations in fixed-salt and osmostat simulations are observed, it is also possible to apply the same SAMS-BAR methodology to correct the free energy calculations that have been performed with fixed salt.</p><p>As our osmostat has been designed to reproduce realistic salt environments around biomolecules, it is well suited to study systems whose function are sensitive to the salt concentration, or biomolecules that are regulated by interactions with Na+ or Cl−. While our osmostat can efficiently sample ion binding to biomolecular surfaces, the sampling of deeply buried ion binding sites is likely to be no more than efficient than in typical molecular dynamics simulations due to the fact that our osmostat is implemented by swapping water with salt. To this end, the osmostat could be improved and generalized if position-biased insertions of fully-decoupled ghost molecules could be added to its sampling repertoire. An example of one such biasing scheme can be found in the biomolecular simulation package ProtoMS, where the grand canonical insertion and deletion of water are attempted in a pre-defined region within proteins48,49. Previously, Song and Gunner studied the interplay between protein conformation, residue pKas, and ion binding affinity using a grand canonical ion insertion scheme within the MCCE framework6. Their work provided structural insight into the often tight-coupling between ion and proton affinity as well as the pH sensitivity of ion binding, and highlights the power of specialized ion sampling schemes to rationalize and understand experimental measurements. The insertion of decoupled ghost molecules—while it would likely require more highly optimized alchemical protocols for insertion—would also permit generalizing the method to more complex salt or buffer molecules or other excipients.</p><!><p>Because the pKa of protein residues are dependent on the ionic strength of the medium, a natural extension of the osmostat is to combine it with constant-pH simulations in explicit water. Previously, Chen and Roux coupled protonation state changes with the insertion and deletion of ions to maintain electrostatic neutrality33,34. The application of an osmostat to such transformations would allow for the macroscopic ion concentration—as well pH—to be rigorously maintained, and could be implemented in modular MCMC scheme that updates protonation states and ion identities in tandem.</p><p>This work only considers the concentration of NaCl, but both the formalism we introduce in the Theory section and the flexibility SaltSwap code-base can be readily extended to sample over biologically relevant salt mixtures by including additional monovalent species such as K+ and divalent species like Ca2+. More complex ions or buffer molecules, such as HCO3− would require a more significant extension to code (such as the insertion of ghost particles described earlier), and could be implemented by using a softcore alchemical NCMC pathway that converts the molecule between fully interacting and noninteracting states.</p><p>The combination of a multicomponent osmostat with a constant-pH methodology would allow for realistic physiological conditions to be better approximated in molecular simulations. While it is well appreciated that pathological tissue can be found with altered pH—tumor microenvironments can have low pH, while cancer cells can have elevated pH, for example89—pathologies can also disrupt healthy ion compositions5. The ability to reproduce specific ionic concentrations as well as pH would open the possibility of using molecular simulations to target compounds to specific microenvironments or achieve selectivity via salt-dependent environmental differences. Indeed, Spahn et al. recently used molecular simulations to develop an analgesic that selectively targets the μ-opioid receptors in damaged, low pH, tissues90.</p><!><p>We have demonstrated that our implementation of the NCMC osmostat was sufficient to sample equilibrium distributions of ions around biomolecules in practical simulation times. We have not yet extensively optimized the osmostat for computational or algorithmic efficiency beyond exploring NCMC protocol lengths (Figure 4 and Figure A5.3), such that there a number of ways that the computational efficiency could be further improved.</p><p>In our current implementation, which only proposes insertion/deletion of a single salt pair in each proposal, the correlation time for the instantaneous salt concentration increases with increasing system size as the size of the equilibrium fluctuations also grow in terms of total numbers of ions (Figure 5). Inserting or deleting multiple ion pairs—likely using longer specialized NCMC protocols tuned to the number of ions being inserted or deleted—could help maintain efficiency. Adaptive MCMC proposals, currently in widespread use in the Bayesian inference community (e.g., PyMC91), could be used to automatically tune the number of ions proposed to be deleted or inserted based on the current concentration and the history of the sampler, provided care was taken to ensure the adaptation method maintained ergodicity and ensured the target density was properly sampled92. One of the earliest adaptive scheme was originally validated on unimodal distributions93, such that a discretized variant could be well suited to sampling the number of pairs.</p><p>Acceptance rates can also be increased by using proposals that do not simply select ions at random, but instead select ions that are more easily inserted/deleted based on some rapidly-evaluated surrogate (such as their instantaneous Monte Carlo acceptance probabilities or the electrostatic potential on water and ion sites), provided this biased selection probability is accounted for in a modified Metropolis-Hastings acceptance criteria.</p><p>There is a great deal of potential to improve the efficiency of the NCMC protocol used for the insertion and deletion proposals. The current work uses a linear interpolation of the salt and water nonbonded parameters as the alchemical path and perturbations steps that are equally spaced with respect to the parameters, primarily because this is the simplest scheme to implement. The only optimization carried out here was tuning the total protocol length to be sufficiently long to achieve high acceptance rates but not so long that the overall efficiency would be diminished by further extending the protocol length (Figure 4). Optimized NCMC protocols can reduce protocol switching times required to achieve high acceptance rates, thereby increasing overall efficiency. The ability to quantify the thermodynamic length of the nonequilibrium protocol allows the problem of protocol optimization to be tackled rigorously. The thermodynamic length (an application of the Fisher-Rao metric to statistical mechanics94) is a natural, albeit abstract, measure of the distance traversed by a system during a thermodynamic driving process84.</p><p>Within this framework, optimal NCMC protocols are given by geodesics in a Riemannian metric tensor space86. The thermodynamic length of the NCMC protocol can be estimated in separate equilibrium simulations spaced along the alchemical path, or estimated directly from the protocol work values of the NCMC switching trajectories, including those from rejected proposals85. For optimizing a preselected alchemical path, spacing the perturbation steps to be equidistant with respect to the thermodynamic length can improve acceptance rates by reducing the total variance of the protocol work. As optimal paths are geodesics in thermodynamic space, the most efficient alchemical path for the insertion or deletion will likely be a nonlinear, rather than linear, interpolation of the water and ion nonbonded parameters. Previous efforts to optimize nonequilibrium paths have included directly solving for the geodesic87, sampling the protocol from an ensemble88, and by restricting the optimization to a family of functional forms95. The close relationship between thermodynamic length and the dissipation along the path also suggests that restricting the propagated dynamics to only the first few layers of the solvation shell around the transmuted molecules could also improve the NCMC protocol.</p><!><p>The philosophy of this work is that increasing the realism of biomolecular simulations will aid structural inference and improve the quantitative accuracy of predictions. We believe that the NCMC osmostat we have presented here will be a useful tool for probing the interactions of ions and biomolecules under more physiological conditions than considered in traditional molecular dynamics simulations. It is our hope that the application of the osmostat to protein-ligand binding free energy calculations and extending the method to more comprehensive ion compositions will improve its utility even further.</p>