A Hybrid Hamiltonian for the Accelerated Sampling along Experimental Restraints

In this article, we present an enhanced sampling method based on a hybrid Hamiltonian which combines experimental distance restraints with a bias dependent from multiple path-dependent variables. This simulation method determines the bias-coordinates on the fly and does not require a priori knowledge about reaction coordinates. The hybrid Hamiltonian accelerates the sampling of proteins, and, combined with experimental distance information, the technique considers the restraints adaptively and in dependency of the system’s intrinsic dynamics. We validate the methodology on the dipole relaxation of two water models and the conformational landscape of dialanine. Using experimental NMR-restraint data, we explore the folding landscape of the TrpCage mini-protein and in a second example apply distance restraints from chemical crosslinking/mass spectrometry experiments for the sampling of the conformation space of the Killer Cell Lectin-like Receptor Subfamily B Member 1A (NKR-P1A). The new methodology has the potential to adaptively introduce experimental restraints without affecting the conformational space of the system along an ergodic trajectory. Since only a limited number of input- and no-order parameters are required for the setup of the simulation, the method is broadly applicable and has the potential to be combined with coarse-graining methods.


Hybrid Hamiltonian for the enhanced sampling of protein folding
We define an ideal biasing Hamiltonian H(B) leading to the real underlying partition of the system represented by a free energy landscape (FEL) along suitable restraint vectors and state that the standard forcefield parameter set describes the Hamiltonian H(A) as a reference energy function. Therefore, any property X as a result from the sampling using the Hamiltonian H(A) leads to a corresponding probability P to match a certain value : An additional bias H(B) in the energy space as conventionally applied in enhanced sampling changes the resulting probability P to match a defined target property X to : The additional bias H(B) is applied along collective variables describing the slowest degrees of freedom of the system, as for example in umbrella sampling [1] and related methods such as Metadynamics [2], conformational flooding [3] or local elevation [4]. We recently have introduced a renormalization scheme to solve the problem that the un-biased probability P(X) can be strongly affected by the added bias H(B) [5]. In this scheme, the applied bias is renormalized to the un-biased Hamiltonian H(A) in a way that its magnitude only equals a fraction of the un-biased Hamiltonian H(A) in dependency of two coupling factors α md and α . Therefor, we introduce a renormalization of the un-biased Hamiltonian by the same factor, which results in a new Hamiltonian H(C) : where we distinguish between α md as coupling parameter renormalizing the un-biased Hamiltonian and α , which defines the coupling of the bias to the system. Alternatively, the restrained Hamiltonian H(B) R leads to : where the new hybrid Hamiltonian H(C) R consists of a reweighted experimental restraint component r. In principle, the total energy remains approximately un-affected, while other properties can be introduced through the bias H(B). In other words, we generate a new probability distribution P (X) as a result from adding a particular bias H(B) in dependency of two linear coupling factors α md and α : which also affects the probabilities P (X † ) † to find the transition state X † , leading to a behavior which can enhance the sampling : Alternatively, the restrained Hamiltonian H(B) R leads to : while the same dependency holds for the probability to find the transition state X † . Through this formalism, we introduce another enhancing fragment H(B) or H(B) R to the original Hamiltonian H(A), without affecting the total energy of the system. As we show later, we generate an accelerating Hamiltonian H(B) or H(B) R from the parallel path-increments added to the system in the form of a renormalized fluctuation. The orientiation of H(B) is explicitly determined based on a definition of a path-increment and acts on the system through a renormalized fluctuation. That procedure has the advantage that the propagation of the system remains ergodic as long as the coupling factors α md and α lie within a low range of approximately < 10 −4 [5]. Considering the fact that bonded interactions in the biomolecular forcefield can contribute with gradients > 10 4 − 10 5 kJ/mol/nm, a coupling with a factor with a magnitude of α md = α = 10 −4 corresponds to the order of magnitude of typical non-bonded interactions. We mention that we used parameters in the range α md = α = 10 −8 for DNA systems using potentials of mean force derived from PDB-data [6].

Pathway dependent biasing increments : ∇H(B)
We define the bias Hamiltonian H(B) used for the accelerated sampling along multiple path increments dL ik (for pathways i and atom-indices k) as well as its modification to its principal components, which is adaptively changed into a bias H(B) R in dependency of a distance restraint r given by experimental data. We consider that the simulated system in an equilibrium simulation propagates along a pathway with the general condition that the reduced action L as function of momentum p and positions q remains constant [7,8] : Along a time-dependent MD trajectory of a system exposed to non-zero fluctuations of its momenta dp, we rewrite the equation 8 as a time-integral : where t 0 stands for the start and t 1 for the end of the simulated MD trajectory. The time-dependent integral is expressed as : We then define a differential dL(t) for a microscopic system, in which fluctuations of the momenta dp(t) occur. That local change in L(t) at the time t is defined by : We obtain the following differential at time t : The expressions in the equations 9 to 12 relate to the standard MD-case with fluctuations in the momenta dp(t) and displacements dq(t) at times t. Any arbitrary biasing technique developed to accelerate a MD simulation adds an instantaneous increment dL s (t) to dL(t) : where we obtain a modified increment dL (t) and additional changes in the momenta resulting from applied bias-energies affect the instantaneous action of a system in order to reach a faster convergence to the statistical average, i.e. the free energy landscape [2,9,10]. In our recent work, we used 2 biasing increments dL s (t) : dL ab (t) (adaptive bias MD) (in the original work, we refer to the variable dL) and dL σ (L ik (t)) (path-sampling) depending from 2 coupling time intervals τ 1 (adaptive bias MD) and τ 2 (path-sampling) in which the gradient has been evaluated [11] : We extend the formalism to the sampling within N R multiple biases (and optionally N S multiple simulations) and redefine the expression 9 to a multiple sampling in multiple bias-paths along N R multiple biases, which the system can undergo simultaneously. The biased simulation then results in the expression of a modified action expression : where the multiple biasing pathway dependent increment dL s m is defined by : 1. Methods

Renormalization
As given in the equations 3 and 4, we employ a renormalization of the vector ∇H(B) or ∇H(B) R to the underlying ensemble [5], derived from multiple biasing increments ∇H(B) ab (adaptive bias MD) and ∇H(B) σ (path-sampling). We apply the relations 3 and 4 for the potential H(B) or H(B) R acting on an individual atom , where the two parameters α md and α allow fluctuations of the coupled bias and the un-biased gradient around an average value : and where η md and η define the width of the fluctuation and ξ is a random number ∈]0, 1]. Using the definition : we add the bias ∇H(B) to the system. As mentioned before, we distinguished between α md to describe the fluctuation of the un-biased MD-gradient and α for the renormalization of the bias. Using the absolute value of H(A), we can formulate the resulting modified absolute value of the hybrid Hamiltonian |H(C)| of the biased system, which is expressed by : If α md = α and α md on the order of 0 < α md < 0.1, then 1+α md +α 2 md 1+α md ≈ 1 , and we obtain : which means that the energy remains conserved for small α md -values. The renormalization to the underlying Hamiltonian leads to a reduction of the added bias s to the norm of the applied un-biased gradient. We summarize that the magnitude of the added bias s derived from multiple path-dependent increments depends from the coupling parameters β md , β and the fluctuation parameters η md , η . Considering the fact that bonded interactions in the biomolecular forcefield can contribute with gradients > 10 4 kJ/mol/nm, a coupling with a factor with a magnitude of β md = β = 10 −4 corresponds to the order of magnitude of typical non-bonded interactions. In general, the renormalization procedure guarantees a sampling in which the occurrence of non-equilibrium configurations is prevented and has the advantage that the applied bias enhances the sampling of the system through the addition of renormalized fluctuations. That means that the underlying Hamiltonian is not overdamped through the application of a bias exceeding the un-biased fluctuations, where essential configurations are potentially not accessed. We tested the influence of the α parameters in simulations of SPC/E and TIP3P water, where we investigated the dielectric and the structural properties of water in dependency of the simulation parameters (see Figure 2 b and the supplementary material : Section -2. Water simulations. supplementary Table 1, Table 2, Figures S1 and S2).

Adaptive bias MD
For the adaptive bias MD section of the algorithm, we derive a history dependent bias ∇H(B) ab of the form : using a number of N R -biases in which the bias is re-evaluated within periods of τ 1 ik for the bias with an index i and atom k. As we introduced in our previous work, we define the corresponding force F b (t) at time t, and use time-derivative of ∇H(B) ab : d dt H(B) ab =Ḣ(B) ab : +γ ik (t)(p k (t) + dp k (t))dq k (t) .
As we defined in the equation 8, the added bias has to fulfill the condition that lim t→∞ dL ab (t) t ≈ 0 in order to sample the system at equilibrium. That also implies that the averages of γ ik have to fulfill that < γ ik (t) >= 0 and dγ ik dt ≈ 0. That relation is fulfilled if : which is implemented by : where γ ik stands for the fluctuation range with the dimension of a length ([nm −1 ]) and ξ is a normally distributed random number with a weight equal 1. We used a constant value γ ik = 10 −4 in all simulations. To enhance sampling along a history-dependent pathway in adaptive bias MD, we employ a coarsening expressed by : +γ ik (t)(p k (t) + dp k (t))dq k (t)) .
By taking into account that d dτ 1 ik γ ik (t)(p k (t) + dp k (t))dq k (t)) ≈ 0 we use that formalism to define the differential over a finite time increments τ 1 ik to coarse-grain the dynamics and to increase the computational efficiency, which leads to an expression for the corresponding force in adaptive bias MD :

Path sampling
In the path-sampling, we use a definition of the reactive coordinate σ ik (t) to determine the biasing segment ∇H(B) σ as function of the increment dL σ ik (L ik )(t), which we define as [11] : dt , equal to the path integral reached by integration till time t for the bias with index i and atom k. In cartesian coordinates : where the height W i and the width δσ i are conventionally parameters chosen to provide computational efficiency and an efficent exploration of the free energy F (∇H(B) σ ). We define the bias component as : That formulation constantly drives the system to explore new configurations along the variable L ik (t) and prevents the system to revisit conformers, which have been sampled previously. In the implementation for an efficient exploration of this space, we note that our algorithm uses the definition [2] : where the times t b ik and t b ik are defined by We apply a variable height of each Gaussian added to an individual variable : where W and ∆E are constants [12] We applied constant values W = 0.1 kJ/mol (∆E = 1000 kJ/mol) in all simulations. In order to accelerate the sampling along H(B) and H(B) R , we re-evaluate the principal components of H(B) in order to sample the system along its slowest modes of the 2 segments H(B) ab and H(B) σ . Therefore, we diagonalize the matrices dL σ ik (L ik ) and dL ab ik (t). The corresponding eigenvectors with the smallest eigenvalue represent the slowest modes and are applied to the system [13].

Results : Water simulations
In this section, we discuss our results from the simulations of 2 water models : TIP3P and SPC/E, with the sampling algorithm along multiple pathways. We tested the influence of the 2 renormalization parameters β, η and η MD on the time-dependent relaxation of the total dipole moment of the system. We also varied the parameters N R and the 2 time-periods τ 1 and τ 2 and tested their influence on the dielectric properties of the 2 water models. As key quantities, we chose the radial distribution function g(r), the dielectric permittivities (static) (0) (frequency dependent) (ω) and the diffusion coefficients as a quantitative measure of the properties of water in dependency of the simulation parameters.
As stated in the main text, we consider the effect of the multiple renormalized biases and the renormalization parameters on the dynamical relaxation behavior of a time-dependent quantity X describing a system, such as the relaxation of the time-dependent total dipole moment M(t) of a water system. Any quantity X(t) in an unbiased simulation follows a time-correlation function A(X)(t), which can be described by an expansion to M monoexponential decay processes with periods τ m , which is defined by the time-behavior of a quantity X, rate constants k m and prefactors A m 0 : In contrast to the equilibrium MD case, a bias coupled to a set of N R biases to N atoms, with τ 1 and τ 2 leading to an actual acceleration in terms of a change in the time-correlation function. That results in a modified relaxation behavior, affecting all dynamical quantities (That can lead to accelerated folding times, modified diffusion constants and re-orientation kinetics of H-bonds or dipoles in the system. There is also an effect on quantities such as the static and dynamic dielectric properties), which we write as an heuristic equation (as described in our simulation results of the dielectic response of SPC/E and TIP3P water) : where N R stands for the number of renormalized biases and k i R is the rate-constant within each bias with index l, and we note that the time t changes to a relative time t rel , which scales linear with an acceleration factor ρ > 1: t rel = ρt, as long as the coupling to the underlying un-biased gradient remains sufficiently low (β ik ≤ 1 × 10 −4 , η = η MD ≈ 25 -see Section : Methods, E: 'Renormalization to the underlying Hamiltonian') in the enhanced sampling simulation. The relation 34 shows that depending on the magnitude and the parameters (β, N R and the coupling times τ 1 , τ 2 ), the Hamiltonian of the system and the time-correlation behavior in the simulation are modified, while the processes depending from the parameters A and rate constants k still are described by modified monoexponential time-dependent decays, since the renormalization and the conditions on adaptive bias MD obey the principle of action as described in the equation 5 in the main text. In other words, dynamical quantities such as dielectric quantities related to dipole fluctuations and in general fluctuation-dependent properties related to a linear response of the system can effectively be varied through the choice of the bias-parameters.

Water simulations : System preparation and analysis
We used SPC/E and TIP3P water simulations in order to test the influence of the new simulation method on the correlation behavior of the system and to validate the property of the algorithm to change dynamic properties while structural properties of water remain approximately unaffected.
For the sampling of dielectric properties of SPC/E and TIP3P water, we filled a box with dimensions 1 × 1 × 1 nm 3 with 33 water molecules (SPC/E or TIP3P). We applied PME electrostatics with a cut-off of 0.4 nm, while van der Waals interactions have been calculated using a shift function with the same cutoff. We checked the size-dependency of our results in simulations of a system with size of 5 × 5 × 5 nm 3 filled with 4124 water molecules (PME electrostatics with a real-space cut-off of 0.8 nm, and a cut-off for the Lennard-Jones 12-6 interactions of 0.6 nm). We used the same system and applied identical parameters to test replica exchange MD (REMD) using 4 replicas at 300 K with an exchange interval every 1,000 time-steps. For the REMD simulations, we omitted the determination of frequency dependent dielectric properties. All simulations have been carried out in the NVT-ensemble using the Nosé-Hoover thermostat (i.e. the Nosé-Hoover equations of motion) with a coupling interval of τ c = 1 ps and a reference temperature of 300 K. The small water-systems have been propagated over 20 ns, while the larger reference systems have been sampled over approximately 5 ns.
For the biased water-simulations, we used a β coupling parameter of 5 × 10 −4 , 20 replicas and τ 1 =5 ps, τ 2 =10 ps. The variational parameters have been chosen with a width that no variation of β and τ 1 , τ 2 has been performed. The gradient optimization have been switched off. We systematically varied η MD and η from 0 to 500.0. For the calculation of frequency dependent properties of the dielectric constant (ω), we calculated the autocorrelation function using the g_dipoles module. We normalized the faster decay of the dipole-autocorrelation function from the biased simulations to the decay in the MD-simulations of SPC/E and TIP3P water using the relative positions of the first zero-value of the autocorrelation functions. In other words, we determined the acceleration factor η of the biased-simulation relative to the un-biased simulation in terms of the decay of the dipole autocorrelation function from 1 to 0 (see Figure 1S), which means that the effective acceleration in the reorientation of the water dipoles introduced by the bias is considered as general acceleration η in the sampling of the underying system. That prodedure leads to the correct frequency dependent dielectric permittivity in comparison with the experiment (see section : Results and Discussion, Water simulations, Dielectric properties).
The static dielectric constant (0) have been calculated using the Kirkwood factor g k as implemented in the gromacs-4.5.5 analysis module g_dipoles, using [14] : and The frequency dependent dielectric constant (ω) has been calculated using the g_dielectric module from GROMACS [15]. That formalism uses a normalized autocorrelation function Φ(t) : and uses the Fourier-Laplace transform [16] : For the case of the sampling along N R multiple biases, we refer to equation 34 and state that the static permittivity can be modified in dependency of the coupling parameters η , β and the number of multiple biases N R . For the 2 processes as described in the equations 35 and 38, the correlation functions then effectively yield different permittivity values. As starting parameters for the fit to the decay of the autocorrelation function, we used a double exponential function with the fit parameters A = 0.5, τ 1 =10 ps, τ 2 =1 ps, (0) resulting from the first analysis using the Kirkwood factor, RF =78.5 and a smoothing over 10 data points.

Dielectric properties
We determined the static permittivities (0) ( (ω) at frequency ω = 0) as function of simulation time in biased and unbiased simulations for SPC/E and TIP3P water. Through the investigation of the dependency from the η MD and the η -parameters on the static permittivity, we determined the property of our algorithm to influence the correlation behavior of the underlying system. We used a static amount of bias-replicas and a constant coupling constant β in order to measure the dependency from the fluctuation dependent parameters η MD and η . Our results are summarized in Table 2.
In the biased simulations of SPC/E water, a parameter of η MD =500, η =1.0 yielded the lowest value of (0)=53. We reach a maximum value in the static permittivity of (0)=119 using parameters η MD =100.0, η =100.0, and obtain a permittivity value of (0)=78.4 using the parameters η MD =30.0, Fig. S1. Normalized dipole autocorrelation functions from MD-and biased simulations of SPC/E and TIP3P water. We scaled the decay time of the autocorrelation result of biased SPC/E and TIP3P water to determine the linear acceleration factor η by which the reorientation of the dipole of water is accelerated (see section Methods in the main text). We find that η =10 yields an optimal fit to the decay of the unbiased simulation and also leads to a good agreement with the experimental frequency dependent permittivity (ω) (see main text). η =30.0 independent from the dimension of the system (1 nm 3 and 125 nm 3 ). In the MD-simulations of the systems with dimensions 1 × 1 × 1 nm 3 and 5 × 5 × 5 nm 3 , we determine values for the permittivity of (0)=68.7 and (0)=66.5 (see Figure 2S c). We mention that the results contain an associated error of approximately ∆ (0) ≈ 1-2.
Using biased simulations of the TIP3P-water model in systems with the same dimensions, we obtain the lowest value for the static permittivity of (0)=55 using η MD =500, η =1.0 . In contrast, we measure a maximum value in the static permittivity of (0)=92 with the parameters η MD =1.0, η =1.0. We obtain a permittivity value of (0)=78.0 using η MD =110.0, η =1.0. In the MD-simulations of the 1 × 1 × 1 nm 3 and the 5 × 5 × 5 nm 3 box we determine values for the permittivity of (0)=93 and (0)=102 (see Figure 2S d). We yielded approximately an identical value in a 5 × 5 × 5 nm 3 box with the same parameters.
We analyzed the frequency dependent permittivities (ω), and observed a decay in (ω) spec from a value of 78.4 to 2.69 within the frequency range from 8E-3 to 1 THz in the simulation of SPC/E water using η MD =30.0, η =30.0. The same function decays 4 times slower in the simulation of TIP3P-water using η MD =110.0, η =1.0 in the frequency range from 8E-3 to 4 THz. In the unbiased MD-simulation of SPC/E water, we observe a decay of (ω) spec from 72 to 3.2 in the frequency range from 2.5E-3 to 7.5E-2 THz, while we measure a decay from 92.6 to 3.1 within the frequency range from 8.7E-3 to 4 THz in the unbiased MD-simulation of TIP3P-water. In the case of (ω) (i.e. the imaginary part of the frequency dependent permittivity), we observe a strong effect of the bias on the location of the maximum along the frequency axis. In the simulation of SPC/E water using η MD =30.0, η =30.0, we find a maximum value of (ω) spec =39.7 at 0.11 THz. For TIP3P water simulated with the parameters η MD =110.0, η =1.0, we find a maximum value of (ω) spec =31.47 at 0.107 THz. The experimental maximum is located at approximately (ω) spec =37.24 and =0.12 THz [17]. In contrast, in the unbiased MD-simulations of SPC/E and TIP3P-water, we obtain maxima in (ω) spec =39.69 at ω=0.031 THz (SPC/E) and (ω) spec =44.5 at ω=0.19 THz (TIP3P). In the analysis of the REMD simulations, we omitted this analysis for the REMD trajectories of SPC/E and TIP3P water. Experimental values for the dielectric spectrum ( (ω)) taken from ref. [17]. For the REMD simulations, we used 4 replicas at an identical target temperature (see section Methods) in order to sample the properties of water at 300 K.

Structural and dynamic properties
We analyzed the site-site specific radial distribution functions (RDF) of SPC/E and TIP3P-water obtained from biased and un-biased simulations of the system with a dimension of 125 nm 3 . We measured the RDFs between water-oxygen, water-oxygen (OW-OW), water-oxygen, water-hydrogen (OW-HW) and water-hydrogen, water-hydrogen (HW-HW). Our results on the structural and dynamic properties are summarized in Table 2.
In unbiased MD-simulations of SPC/E water the RDFs contain maxima ((water oxygen-oxygen)(OW-OW) 0.276 nm (1), 0.45 nm (2) and 0.69 nm (3) [17][18][19][20][21]. s . In the MD-simulation of TIP3P water, we measure a self-diffusion coefficient D with a value of 2.3045 (± 0.1547) 10 −5 cm 2 s , while the enhanced simulation of the same system using η MD =120.0 and η'=1.0 results in a value for D of 17.0078 (± 5.2154) 10 −5 cm 2 s (see Figure 2S f). The self-diffusion coefficient observed in the REMD simulations exceeds the values from the sequential runs, where we observe a value of 10.3925 10 −5 cm 2 s for SPC/E water and a value of 50.5028 10 −5 cm 2 s for TIP3P water.
SPC/E-water are in good agreement with their results. Investigations on the dielectric properties of water are of specific interest and especially the SPC/E water model has been found to describe frequency dependent dielectric permittivities shifted by -0.1 THz to the experimental maximum of (ω) [17], while a value of 70.7 [23] for the static permittivity value (0) has been reported. That value is shifted in comparison to the experiment by a value of approximately 8, while the experimental value is at 78.4 at room temperature [21]. In contrast to the SPC/E water model, static permittivities in the range from 93 to 104 have been observed for the TIP3P water model, while its peak for in the dielectric spectrum is shifted by approximately +0.1 THz in relation to the experimental value [17,24]. In terms of the static and frequency dependent permittivity value, our result from the biased simulation of SPC/E water is in good agreement with the experiment. That also holds for the REMD simulation in the case of the static permittivity, while we omitted the frequency dependent analysis due to the dependency of the dynamics of the system from the exchange frequency. The dielectric behavior of water plays an important role in the behavior of biomolecules, while static dielectric constants for biomolecular systems ranging from 10 to 20 have been reported [27,28]. Especially the non-polar region of proteins, which is shielded from the solvent in its folded state, can contain values for (0) of approximately 5 [29]. However, here it remains questionable to speak of a dielectric constant, which is a macroscopic quantity and might not be reducable to a range below the nm length-scale. Recent studies investigated the effect of ions and the folding of alanine-pentapeptide on the dielectric relaxation behavior of the underlying systems [30,31].
We conclude that the presented algorithm can change the time-dependent properties of SPC/E and TIP3P water through a scaling of η MD and η , as we stated in the main body of the text in the Section : Methods. Especially the parameter η MD can modify the time-dependent decay of the total dipole moment in comparison with the conventional MD propagation. We found that the magnitude η affects the fluctuations of the applied bias averaged and also can change with the number of parallel replicas, as we find in the REMD simulations. For the TIP3P model, the charge density has a high value that a large value for η MD =120.0 has to be used to scale down the static permittivity to a value of 78. That scaling has a significant effect on the water structure in comparison with the experimental result and shifts the location of maxima in the RDF. In the case of the SPC/E water, a scaling of η to larger values yields the correct static permittivity value, an approximately correct water-structure in comparison with the experiment as well as a comparatively good agreement in its dielectric relaxation behavior with experimental results. As we find in the REMD case, we observe that the combination of REMD with the multiple path-sampling can sample the structural properties of water in terms of the radial distribution function and the static permittivities in the same way as the sequential technique.

Restraint parameters
• Experimental NMR-NOESY restraint data used for the restraint simulations of the TrpCage miniprotein from ref. [33] (restraint-index, residue number 1, atom-name, residue number 2, atom-name, distance (Å)): • Restraint data set used for the restraint simulations of the Killer Cell Lectin-like Receptor Subfamily B Member 1A from ref. [32] (restraint-index, residue number 1, atom-name, residue number 2, atom-name, distance (Å)). We note that we considered each distance-value to be significantly lower than the chemical cross-linking distances as described before (7.5 Å). The distance value only affects the sign of the overlappling segment of the path-increment dL in contrast to harmonic potentials: