Recovering position-dependent diffusion from biased molecular dynamics simulations.

Recovering position-dependent diffusion from biased molecular dynamics simulations Ajasja Ljubetič, Iztok Urbančič, and Janez Štrancar Citation: The Journal of Chemical Physics 140, 084109 (2014); doi: 10.1063/1.4866448 View online: http://dx.doi.org/10.1063/1.4866448 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/8?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Combining molecular dynamics and an electrodiffusion model to calculate ion channel conductance J. Chem. Phys. 141, 22D519 (2014); 10.1063/1.4900879 Peptide dynamics by molecular dynamics simulation and diffusion theory method with improved basis sets J. Chem. Phys. 140, 104910 (2014); 10.1063/1.4867788 Investigation of polarization effects in the gramicidin A channel from ab initio molecular dynamics simulations J. Chem. Phys. 137, 205106 (2012); 10.1063/1.4768247 BDflex: A method for efficient treatment of molecular flexibility in calculating protein-ligand binding rate constants from Brownian dynamics simulations J. Chem. Phys. 137, 135105 (2012); 10.1063/1.4756913 Diffusive dynamics on multidimensional rough free energy surfaces J. Chem. Phys. 127, 135101 (2007); 10.1063/1.2775444

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 128.42.202.150 On: Sun, 23 Nov 2014 23:54:44

THE JOURNAL OF CHEMICAL PHYSICS 140, 084109 (2014)

Recovering position-dependent diffusion from biased molecular dynamics simulations Ajasja Ljubetiˇc, Iztok Urbanˇciˇc, and Janez Štrancara) Laboratory of Biophysics, Condensed Matter Physics Department, “Jožef Stefan” Institute, 1000 Ljubljana, Slovenia

(Received 29 November 2013; accepted 6 February 2014; published online 27 February 2014) All atom molecular dynamics (MD) models provide valuable insight into the dynamics of biophysical systems, but are limited in size or length by the high computational demands. The latter can be reduced by simulating long term diffusive dynamics (also known as Langevin dynamics or Brownian motion) of the most interesting and important user-defined parts of the studied system, termed collective variables (colvars). A few hundred nanosecond-long biased MD trajectory can therefore be extended to millisecond lengths in the colvars subspace at a very small additional computational cost. In this work, we develop a method for determining multidimensional anisotropic position- and timescale-dependent diffusion coefficients (D) by analysing the changes of colvars in an existing MD trajectory. As a test case, we obtained D for dihedral angles of the alanine dipepR package, capable of determining and visualizing D in one tide. An open source Mathematica or two dimensions, is available at https://github.com/lbf-ijs/DiffusiveDynamics. Given known free energy and D, the package can also generate diffusive trajectories. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4866448] I. INTRODUCTION

Using an all-atom classical molecular dynamics (MD) model, the time evolution of large systems with millions of atoms can be simulated.1 Such an approach offers valuable insight into molecular processes, but is limited by simulation length and insufficient sampling duration. Classical MD requires femtosecond time steps to integrate fast atomic motions, which increases the computational cost and currently limits simulation length to microseconds. Additionally, high energy barriers are crossed infrequently, making it difficult to sample all relevant configurations. To avoid the limits of classical MD, the system’s time evolution can be simulated along a few chosen collective variables, usually denoted as colvars. Colvars are explicit functions of microscopic degrees of freedom in the system that capture the most relevant long-term dynamic behavior, for example, an angle between three atoms, a dihedral angle or a distance between the center of mass of two molecules. The reduction in the degrees of freedom drastically lowers the computational cost and extends the duration of trajectories to milliseconds. The long-term time evolution of the colvars is approximated by Brownian motion, also known as diffusive dynamics.2 Brownian motion describes the time evolution of slowly changing degrees of freedom that are coupled to fast changing degrees of freedom. The archetypical example is the jittery motion of a grain of pollen in water,2 where the motion of the grain represents the slow degrees of freedom and the motion of water represents the fast degrees of freea) Author to whom correspondence should be addressed. Electronic mail:

[email protected]

0021-9606/2014/140(8)/084109/11/$30.00

dom, coupled to the slow degrees by collisions. The model of Brownian motion can be used to describe a large variety of other phenomena: reaction-rate theories,3, 4 including sitespecific enzyme kinetics,5 formation of micelles,6 and even laser cooling;7 further applications are reviewed by Frey and Kroy.8 In MD, the colvars represent the slow degrees (the “grain”) and the neighboring atoms represent the fast degrees (the “water”). Brownian motion is described by system and random forces. System forces originate from potentials present in the system. In MD, the most common potential is free energy; therefore several methods are readily available for obtaining the free energy profile along a few chosen colvars. The methods apply a bias potential to traverse large energy barriers, and are thus referred to as biased MD. Among the most established methods are: Meta-dynamics,9, 10 Adaptive Biasing Force (ABF)11, 12 calculations, steered molecular dynamics (SMD),13, 14 and umbrella sampling.15 The bias perturbs system forces and thereby the normal dynamics of the system, meaning that it is possible to obtain either the dynamics or the free energy. Random forces are caused by the coupling between the slow and fast degrees of freedom and are described by the diffusion coefficient (D). In the simplest case, D is positionindependent and isotropic (a scalar). In general, however, D is a position-dependent anisotropic quantity (a tensor).The position dependence is a consequence of local environment variability and of the projection of several degrees of freedom onto fewer colvars.16 Consequently, D is more difficult to estimate than the free energy. A few methods for the positiondependent determination of D are available,17–20 but were at present only applied to one dimensional examples. There is also no freely available software implementation.

140, 084109-1

© 2014 AIP Publishing LLC


084109-2

Ljubetiˇc, Urbanˇciˇc, and Štrancar

J. Chem. Phys. 140, 084109 (2014)

Herein we present a method for determining multidimensional anisotropic timescale- and position-dependent diffusion coefficients by analyzing the motion of the colvars in an existing trajectory. To assess the ability of the algorithm to reconstruct the D, we first analyzed the trajectory generated by a simple two dimensional synthetic model system with known diffusion and free energy surfaces. Then we analyzed the rotation of the ϕ and ψ dihedral angles of the alanine dipeptide, which provides a well-studied and convenient model for thermodynamic and kinetic studies.21–23 We have used ABF MD to improve sampling and shown that ABF does not modify D, making it possible to obtain both the free energy and D from the same set of simulations. II. THEORY A. Brownian motion

Brownian motion is described by the Langevin stochastic equation24 mi x¨i = − γij x˙j + fiS + fiR , (1) j

where xi denotes the ith component of particle’s position, and mi is the mass of the particle. The dots denote the time derivatives and γ ij is the friction tensor. The symbol fi S represents the system forces, for example due to the electric field or another potential. fi R describes the random forces representing thermal fluctuations from the coupling of the slow with fast degrees of freedom. In the case of MD, the system forces originate from a position dependent free energy potential, fiS = −

∂F ( x) . ∂xi

(2)

F can be obtained by using any of the methods described in the Introduction. The random forces (fi R ) are assumed to follow an uncorrelated Gaussian distribution with zero mean and 2kB Tγ ij covariance matrix as given by Eqs. (3) and (4), R (3) fi (t) t = 0,

fiR (t )fjR (t) t = 2kB T γij δ(t − t ),

(4)

where kB is the Boltzmann constant, T the absolute temperature, and δ(t) the delta function. The 2kB Tγ ij covariance matrix is defined in terms of temperature and friction by the fluctuation-dissipation theorem.25 On longer timescales, random forces dominate the dynamics causing the particle to behave as though it has a very low mass or is experiencing very high friction. This low mass/high friction limit, often called the diffusion regime, can be approximated by the overdamped Langevin equation γij x˙j = fiS + fiR . (5) j

A numerical scheme for simulating the diffusive Brownian motion, described by Eq. (5), was first presented by Ermak and McCammon26 (see Eq. (15) in the mentioned paper; a more intuitive derivation is presented in section three of

Grassia et al.16 ), x) xi,t+t = xi,t + Ri,t ( ∂Dki ( xt ) 1 ∂F ( xt ) t, − Dik ( xt ) + k k ∂xk kB T ∂xk (6)

Ri,t ( x ) t = 0, x )Rj,t ( x ) t = 2Dij ( x )δ(t − t ), Ri,t (

(7)

where xi,t an xi,t+t are the positions at time t and t+t, Dij = kB T /γ ij is the diffusion tensor, and Ri,t the displacements steps due to random forces. The displacements follow an uncorrelated Gaussian distribution with zero mean and 2Dij covariance matrix as shown in Eq. (7). The steps (xi,t+t − xi,t ) also follow a Gaussian distribution with the same covariance (and a possibly non-zero mean), therefore the diffusion coefficient can be obtained from the covariance matrix of the steps’ distribution. However, when Dij is position dependent, additional care should be taken not to average out Dij over too large an area. B. Determination of position dependent diffusion coefficient

To determine position dependent Dij , a trajectory of points in colvars space sampled at discrete time intervals (t) is needed. The trajectory can be obtained by methods such as molecular dynamics simulations or fluorescence microscopy. To account for the position dependence of Dij , the colvars space is divided into bins as shown in Fig. 1(a). All continuous segments in the bin are identified (depicted with different colors in Fig. 1(b)), and from the segments the steps are obtained. From the distribution of step sizes, Dij is calculated (Fig. 1(d)) as described below. For a truly diffusive process, Dij does not depend on the timestep (t) at which points are sampled. Thus, the timescale dependence of Dij can be used to check the validity of approximating the long-term dynamics of the observed system by diffusive overdamped Brownian motion. In our algorithm, the timescale dependence is introduced by defining the step as a distance between every sth point instead of between consecutive points, xi (n, s, k) = xi (n + s, k) − xi (n, k),

(8)

where n is the point index (i.e., first, second. . . ), k is the number of the bin to which the step belongs and s is the stride parameter that determines how many consecutive points are taken as one step, thereby setting the timescale to st. At large sampling timescales, the size of the steps can become similar to the chosen size of the bin. In such case, very few steps remain in that bin, with larger ones being unproportionally scarcer. Consequently, Dij is poorly and improperly estimated. We have developed two binning methods to ensure unbiased sampling; both effectively work by taking into account points outside the bin (Fig. 1(c)). In the first method, binned continuous trajectories that contain fewer points than the current stride (s) value are


084109-3


J. Chem. Phys. 140, 084109 (2014)

FIG. 1. The basic outline of the method: (a) First we divide the colvars space into discrete bins and (b) then we find all the contiguous trajectories crossing each bin. (c) Optionally if the bins are very small in comparison to steps, additional care must be taken to sample the distribution of steps correctly and points from outside the bin (as indicated by the dashed outline) must be taken into account. (d) From the distribution of step’s lengths (differences between two positions in colvars space), the diffusion coefficient for each bin can be calculated. In two dimensions, the tensor can be represented with an ellipse or, equivalently, with a rotation (α D ) of the colvar orthogonal coordinate system around the origin and a rescaling of the rotated vectors by Dx˜ and Dy˜ as given by Eqs. (16) and (17).

padded by adding s points to the start and end of the trajectory, thereby ensuring that each trajectory has at least 2s+1 points and therefore contains at least two steps at stride s. If two trajectories overlap as a result of padding, they are merged into a single continuous trajectory. In this way, padding reduces the number of continuous trajectories and increases their length. The second method is conceptually orthogonal. Instead of binning points, steps are binned. Ideally, steps would be assigned to the bin with which they have most overlap. For the sake of algorithmic simplicity, each step is assigned to the bin that contains the midpoint of the segment. In this way, each step is assigned to exactly one bin, which is not the case with padding. This method is also much easier to implement numerically. The size of the bin should be chosen so that the diffusion and free energy are nearly constant inside each bin. At the same time, the bins should be large compared to the typical diffusion length scale, so that a sufficient number of steps remain in the bin to ensure efficient statistics. If trajectories are padded, smaller bin sizes can be used, since large steps are not lost. For systems such as backbone dihedral angles, about 10 bins per colvar is a reasonable amount, but the quantity required can vary depending on the roughness of the diffusion and energy surface. From the steps obtained in each bin (k), the corresponding diffusion coefficient (Dij ) is calculated by estimating the covariance matrix (σ ij 2 ) of the Gaussian distribution

σij2 (s,k)

N N 1 1 = xi (n,s,k) xi (n,s,k) − N n N n N 1 × xj (n,s,k) − xj (n,s,k) , (9) N n

where the stride (s) determines the timescale dependence, the bin index (k) the positional dependence, and n iterates over all steps in bin k. The diffusion coefficient is then given by Dij (s,k) =

σij2 (s,k) 2st

,

(10)

where st is the timescale at which steps are sampled.

C. Checking the validity of the diffusive approximation

The basic assumptions of Brownian motion are: steps are uncorrelated in time and follow a Gaussian (normal) distribution. None of the assumptions are enforced by the analysis method and can therefore be used to check the validity of the diffusive approximation. To test the correlation of steps, the autocorrelation function is used, N/s−1 D 1 1 xi (ns,s) · xi (ns + τ ,s) , C(τ ,s) = I i=1 N n=1 (11) where I is the number of components and N the number of steps xi . Usually the autocorrelation function is normalized by dividing the function by C(0,s). Getting positiondependent autocorrelations is difficult, as at large strides (even with padding) long trajectories inside a bin are very uncommon. But for testing correlation, position dependence is not necessary and the autocorrelation function can easily be calculated for the entire trajectory at different strides. To test the normality of the distribution, we have chosen the Anderson-Darling27 test due to its higher sensitivity for deviations in the tails of the normal distribution as compared to other tests, for example, the Cramér–von Mises test. At


084109-4


each stride in each bin, we obtain the p-value. A low p-value indicates that the data are likely not normally distributed. A similar strategy is used by Micheletti et al,19 except that we use the Anderson-Darling test instead of kurtosis to estimate normality.

J. Chem. Phys. 140, 084109 (2014)

relation, which states that the diffusion coefficient is the time integral of the VAF (the derivation is given in Frenkel and Smit,30 p. 89),

τ (t)dt, (14) D(τ ) = 0

D. A model for the timescale dependence of the diffusion coefficient

A model for the timescale dependence of the diffusion coefficient can be constructed by using the velocity autocorrelation function (VAF, ), which measures the timescale on which the velocity is no longer self-similar, t v(t)v(t + τ )dt , (12) (τ ) = 0 t 2 0 v(0) dt where v(t) is the velocity at time t. A single particle in vacuum has a constant VAF, as its velocity remains unchanged. In a gas collisions between particles decorrelate velocities, causing their VAF to decay exponentially. Atoms in solids fluctuate around their fixed positions, resulting in long-term oscillations of the VAF. The most complex VAFs are observed for particles in a liquid state, where collisions cause the VAF to decay over long time periods, but on a shorter timescale a transient “cage of solvent” forms, from which the observed particle can rebound, causing oscillations in the short time VAF behavior.28 A simple model for a VAF of particles in a liquid can be obtained by considering a Brownian particle trapped in a harmonic potential,29 τ ω τ − ττ 0 cos + sin , (13) (τ ) = e ω0 τ0 ω0 where τ 0 is the decay time and ω0 the oscillation frequency, related to the strength of the harmonic potential. Diffusion coefficients are connected to the VAF by the Green-Kubo

where D(τ ) is the diffusion coefficient at timescale τ . Using the Green–Kubo relation, Eq. (13) can be integrated to obtain a model for the timescale dependence of D, τ −τ D(τ ) = D ∞ − D ∞ e τ0 cos ω0 2 τ − ω02 − ττ τ , (15) + D∞ 0 e 0 sin 2τ0 ω0 ω0 where D∞ is the diffusion at an infinite timescale. If one wants to estimate the long term diffusion dynamics of their system, D∞ can be obtained by fitting expression (15) to the D(τ ) points obtained by our algorithm. The inverse square of the errors were used as weights in the nonlinear fitting procedure. E. The two-dimensional case

In two dimensions, the diffusion tensor can be intuitively presented by an ellipse whose semi-axes (Dx˜ and Dy˜ ) are rotated with respect to the colvar coordinate frame (x, y) by an angle α D (Fig. 1(d)). The semi-axis lengths Dx˜ and Dy˜ are also the eigenvalues of the diffusion tensor,

1 2 + 4D 2 − 2D D + D 2 ) Dx˜ = (Dxx + Dyy − Dxx xx yy xy yy 2 (16)

1 2 2 2 Dy˜ = (Dxx + Dyy + Dxx + 4Dxy − 2Dxx Dyy + Dyy ), 2 where the subscripts x and y signify the first and the second dimension in the colvar coordinates, respectively. For simplicity, the stride dependence is not explicitly marked. The angle of rotation (αD ) is given by

⎛

⎞

2 ⎜ ⎟ 2 −D Dxx − Dxx − Dyy + 4Dxy ⎜ ⎟ yy −1 ⎜ ⎟. αD = cos ⎜ ⎟

2 ⎝ ⎠ 2 2 8Dxy − 2 Dxx − Dyy −Dxx + Dxx − Dyy + 4Dxy + Dyy

III. METHODS A. Synthetic model system

In order to investigate how the size of the bin and the number of steps in the trajectory affect the fitted diffusion coefficients, we created a trajectory of points using a synthetic model. The diffusion was defined by the following expressions, which mimic the expected roughness and range of diffusion values of the backbone dihedral

(17)

angles, x + 10, Dx˜ (x,y) = 5 Sin 2π Lx y + 10 Dy˜ (x,y) = 5 Sin 2π Ly x 2 y 2 αD (x,y) = 180◦ + 180◦ , 4Lx 4Ly

(18)


084109-5


J. Chem. Phys. 140, 084109 (2014)

FIG. 2. Evaluation of the method with a synthetic explicitly diffusive model. Fitted diffusion tensors, where different binning methods were used, are compared to model values (black). The methods are: “w/o padding” (orange), only the points are binned; “padding” (blue), short trajectories are padded with additional points; and “midpoint” (green), steps are binned based on the midpoint of the step. The model values used to create the trajectory are given by Eq. (18). The unit cell, 4000 × 8000 space units wide, was divided into 10 × 10 bins. The timescale is unitless in this synthetic models system. Error bars in panels (a)–(d) represent standard deviation, obtained from ten independent runs with 5 × 106 points. (a) Root mean square deviation (RMSD) of difference between fitted and model diffusion surfaces as defined by Eq. (20) versus sampling timescale. Using padding (blue line), our method successfully reconstructs (fits) diffusion coefficients over timescales differing by three orders of magnitude as indicated by low RMSD. Binning steps by their midpoint (green line) compares very favorably with padding, while being easier to implement. Binning only points (orange line) gives unsatisfactory results. (b)–(d) The components of the diffusion tensor vs. timescale in the bin marked with a gray rectangle in panels (e)-(f). With longer timescales, the steps (differences between two points) increase in size. Steps comparable to the size of the bin have a lower probability of being inside the bin completely, and are therefore under sampled; consequently, diffusion is underestimated (orange line). Padding continuous trajectories with points from outside the bin (Fig. 1(c)), ensures large steps are correctly sampled and significantly improves the fitted diffusion coefficients (blue and green lines). (e)-(f) Comparison of model and fitted position-dependent diffusion coefficients at short (s = 1) and long (s = 1024) timescales, respectively. Diffusion coefficients are depicted by an ellipse in each bin as in Fig. 1(d). At short timescales, the fits match very well. At longer timescales, care must be taken to correctly bin points and ensure unbiased sampling.

where Lx and Ly are the sizes of the periodic cell and Dx˜ , Dy˜ , and α D are the components of the diffusion tensor as defined by Eqs. (16) and (17). In addition to the diffusion coefficient, free energy is needed to simulate diffusive dynamics. In the synthetic model (Figs. 2 and 3), we have used an expression which approximates the residual free energy after an ABF MD simulation, F (x,y) = 2kB T Sin(20x/Lx )Cos(20y/Ly )/3.

(19)

The quality of the reconstructed diffusion surface was evaluated by comparing the model diffusion to the fitted values in all bins by a root mean square deviation (RMSD), 2 2 K 2 1 model Dij RMSD = (k) − Dij fit (k) , K k i j

(20)

where K is the number of bins, Dij model (k) the model and Dij fit (k) the fitted diffusion tensor in bin k. B. Rotation of dihedral angles of alanine dipeptide in water

We analyzed the diffusive dynamics of dihedral angles ϕ and ψ of the alanine dipeptide.21–23 Strictly speaking, “ala-

nine dipeptide” is an incorrect, yet common, designation for an alanine residue blocked by an acetily and a NH-metyl moiety, or N-acetly-N -methy-L-alanylamide (NANMA). The alanine dipeptide was inserted into a 25 Å cube of water and equilibrated for 1 ns at pressure of 1 bar, at a temperature of 300 K using a Langevin thermostat with 1/ps damping. To obtain a good estimate for free energy, four independent 100 ns ABF11, 12 MD simulations at 300 K and constant volume were calculated. After 100 ns the free energy converged, meaning that the bias applied by ABF no longer changed. The timestep of the integrator was 1 fs and colvars were written to a file each 1 ps. Because dihedral angles rotate on sub-picosecond timescales, colvars needed to be sampled very frequently in order to determine diffusion coefficients. To this end, we performed ten independent 1 ns ABF MD simulations, using the previously obtained ABF bias and starting coordinates taken randomly from the previous 400 ns of ABF simulation. The colvars were sampled each 1 fs. Ten 1 ns simulations were performed, since one 10 ns simulation does not efficiently sample the entire ϕ–ψ colvars space. To investigate the effect of ABF MD on diffusion, we also simulated ten independent 1 ns trajectories of unbiased MD, with starting coordinates randomly taken from the available 400 ns of ABF MD simulation. The colvars were also sampled each 1 fs. The system


084109-6


FIG. 3. RMSD between model (actual) and fitted diffusion coefficients versus the number of bins (a) and the number of steps in one bin (b). RMSD is defined by Eq. (20). Error bars represent one standard deviation from ten independent runs. (a) At short timescales (dashed lines) or, equivalently, small steps compared to bin size, increasing the number of bins improves the fits as indicated by the decreasing RMSD. The choice of the binning method is irrelevant (all dashed lines overlap). At long timescales (solid lines) or large steps, binning only points (orange line) into many small bins increases the RMSD as large steps are not properly sampled. Here the binning method becomes important. Padding the trajectories in bins (blue lines) or binning steps by their midpoint (green line) drastically improves the fits at long timescales. Five million points were used in each trajectory. (b) Increasing the number of points in the trajectory improves obtained diffusion coefficients at short and long timescales (dashed and solid lines, respectively). At long timescales, padding (solid blue line) or binning by midpoint of steps (solid green line) significantly improves the fitted diffusion coefficients (10 × 10 bins were used in the calculation).

was prepared with VMD,31 the simulations were calculated with NAMD 2.932 using CHARMM 2233 all atom force field with CMAP34 correction. IV. RESULTS AND DISCUSSION

First we test the presented method on an explicitly diffusive model system. We then apply our method to a more realistic system: the rotation of the alanine dipeptide dihedral angles in explicit water. We then discuss the advantage of the diffusive approximation and compare our approach with existing methods for calculating position-dependent diffusion.

A. Synthetic model system

To test the accuracy of the presented method, we calculated diffusion coefficients from trajectories generated by a synthetic explicitly diffusive model system. The root mean square deviation (RMSD), defined by Eq. (20), between the actual (model) and calculated (fitted) value, assesses the quality of the fitted diffusion coefficients. Small RMSD values (below 1 in our examples) indicate that the diffusion coefficient was accurately fitted (Fig. 2(a)).

J. Chem. Phys. 140, 084109 (2014)

In an existing trajectory, the time step (t) at which colvars are sampled is usually fixed. The colvars can be sampled at larger timescales by sampling every sth point, thus increasing the timescales to multiples of t. For regular diffusion processes, the diffusion coefficient is independent of the sampling timescale (st). When analyzing trajectories at longer timescales, steps increase in size. Steps that are larger than the size of the bin are not taken into account if only points are binned, which leads to underestimating the diffusion (Figs. 2(b)–2(d) orange line). To solve this problem, we have developed two binning methods. Padding the trajectory with extra points outside the bin improves the fitted diffusion coefficients (Figs. 2(b)–2(d) blue line), but as a consequence, the diffusion coefficient is averaged over a slightly larger area. Binning steps by their midpoint also ensures unbiased sampling (Figs. 2(b)–2(d) green line). The method is simple to implement and ensures each steps is assigned to only one bin. In theory, the midpoint methods averages over a smaller area than padding. Padding can reach at most a full step away from the bin border, while the midpoint method reaches at most half a step. Since the sampling timescale is important, fitted and model position dependent diffusion coefficients are compared at short (s = 1) and long (s = 1024) sampling timescales in Fig. 2(e) and Fig. 2(f), respectively. The diffusion tensor in each bin is represented by an ellipse as in Fig. 1(d). At short sampling timescales, the model (black ellipses) and fitted Dij are in excellent agreement (RMSD is below 0.5) regardless of the binning method used. At longer timescales, diffusion is underestimated if only points are binned (orange ellipses). Padding (blue ellipses) or the midpoint method (green ellipses) improves the fit, but averages outside the bin, as expected. To check the effect of the energy surface on diffusion determination, we have tested several different energy surfaces and found that barriers smaller than 3 kB T do not affect the reconstructed diffusions (see Figs. S1(a) and S1(b) of the supplementary material42 ). Thus, the exact shape of the energy landscape does not play a pivotal role in reconstructing the diffusion coefficient. Higher barriers were not examined, as ABF normally reduces the energy differences to below 3 kB T. In general, insufficient sampling of higher energy regions could hampers diffusion determination. The energy surface can also influence the steps distribution; however, given that the bins are small enough that ∂F(x)/∂xi is nearly constant, the free energy only adds on offset to the average of the steps distribution, but does not affect the variance (and therefore does not perturb the diffusion) as can be seen from Eq. (6). We also checked if the model system satisfies the assumptions of Brownian motion. As expected, the steps autocorrelation function falls to zero after the first step for all strides (data not shown) and the p-value is very high for all strides in all bins (see Fig. S1(c) of the supplementary material42 ). In order to determine optimal binning for further analysis, we also explored the effect of the number of bins on the RMSD at short (s = 1) and long (s = 1024) sampling timescales (Fig. 3(a) dashed and solid lines, respectively). A small number of large bins average over a large part of the


084109-7


J. Chem. Phys. 140, 084109 (2014)

FIG. 4. Diffusion coefficients calculated from ABF MD and unbiased MD trajectories of ϕ and ψ dihedral angles of alanine dipeptide. (a) Diffusion coefficients obtained from ABF MD and unbiased MD (blue and orange ellipses, respectively) at sampling timescale 46 fs. (b) Free energy surface along ϕ and ψ dihedral angles calculated by ABF MD. Low energy regions are depicted with blue, high energy regions with red. Since diffusion coefficients obtained from ABF MD and unbiased MD trajectories match in corresponding low energy regions (the ellipses in panel (a) overlap in the regions that are blue in panel (b)), we conclude that the applied ABF bias does not change diffusive behavior. Unbiased MD (orange in panel (a)) has severe sampling problems in higher energy regions (red in panel (b)) and in some bins it cannot be determined. (c) and (d) The timescale dependence of Dx˜ and Dy˜ in the bin marked in panel (b). Dashed vertical line indicates the sampling timescale 46 fs at which diffusion coefficients in panel (a) are presented. This timescale is presented as it is long enough to be in the diffusive regime but short enough not to be contaminated by neighboring bins.

diffusion surface. In the limit of one bin, only an average value over the entire surface would be obtained and the RMSD for non-constant diffusion would be non-zero. On the contrary, when small bins are used at large timescales, care must be taken to preserve unbiased sampling of long steps. Binning only points (orange solid line) underestimates diffusion, as large steps are unlikely to fit entirely into a bin, resulting in large RMSD. Padding (blue solid line) or binning by midpoint (green solid line) ensure that large steps are correctly sampled. Padding results in lower RMSD than the midpoint method at high number of bins. When many bins are present, few steps are present in each bin and assigning a step to multiple bins and averaging over a larger area becomes beneficial. As expected the algorithm is not very sensitive to the number of points per bin (Fig. 3(b)), after the sampling exceeds a certain threshold, in our case being about 50 thousand steps per bin. Based on the findings above, we have chosen the parameters for further analysis of MD data: 10 bins per colvar were used and padding was applied as the binning method. B. Rotation of dihedral angles of alanine dipeptide in water

In the synthetic model mentioned above, the dynamics were predefined as diffusive by Eq. (6). In MD, the dynamics

are governed by Newton’s laws of motion. To explore Brownian motion in such systems, we performed biased (ABF) and unbiased MD simulations of the alanine peptide21–23 in water and monitored the time evolution of the ϕ and ψ dihedral angles. It is assumed that the (ABF) bias potential only affects the free energy potential along the chosen colvars while all the other fast degrees of freedom, which are the source of Brownian motion, are not affected by the applied bias potential. We have compared the diffusion coefficients obtained from biased (ABF) and unbiased MD simulations (Fig. 4(a); orange and blue, respectively). As expected, the applied ABF bias does not change the diffusion coefficients. In low energy regions (Fig. 4(b); blue areas) of the ϕ−ψ space, the diffusion coefficients from biased and unbiased MD are nearly identical. In high energy regions (Fig. 4(b); red areas), the advantage of determining diffusion coefficients from biased MD becomes clear. Unbiased MD simulations cannot effectively sample high energy regions, so the diffusion coefficients in those regions are very poorly determined, if they can even be determined at all. Alternatively, biased MD samples all regions of the colvars space equally, so that diffusion coefficients in all regions of the colvars space can be reliably determined. These results agree nicely with the previous work of Comer et al.,18 who have also concluded that ABF does not perturb diffusion


084109-8


and is a useful tool for obtaining the sampling necessary for an accurate diffusivity calculation in a short simulation. As MD systems are governed by Newton’s laws of motion, the dynamics of such systems is diffusive only at longer timescales, to allow for enough collisions between the atoms describing the fast and the atoms describing the slow degrees of freedom to occur. At very short timescales, the motion of the colvars is ballistic, i.e., they follow smooth trajectories (see Fig. S2 of the supplementary material42 ). The transition from ballistic to diffusive movement was thought by Einstein to be impossible to measure experimentally,35 but has recently been measured for particles in air36 and water,37 with results that are in excellent agreement with theoretical predictions. The transition from the ballistic to the diffusive regime can also be observed in our MD simulations (Figs. 4(c) and 4(d)). At short sampling timescales, the colvars follow a smooth trajectory, meaning that their velocity does not change much between several consecutive time steps. This is reflected in the linear dependence of D on the timescale τ at short timescales. At longer timescales, when the fast degrees of freedom have sufficient time to interact with the colvars, the motion becomes diffusive. In this regime D is constant, i.e., it does not depend on the timescale. The oscillations in D(τ ), at intermediate timescales between the linear and the constant regime appear because the atoms that define the colvars rebound from their first shell neighbors. This effect, also called the “transient cage” effect, has been observed in several studies of velocity autocorrelation functions in liquids.28, 29, 38 The applicability of the diffusive approximation can be checked after D has been calculated. In order to fulfill the assumptions of Brownian motion, the steps must be uncorrelated in time and normally distributed. The normalized autocorrelation function of steps for the ABF MD trajectory of the previously discussed alanine dipeptide falls rapidly to zero (Fig. 5(a)), indicating that the steps are indeed uncorrelated. The p-value of the Anderson-Darling27 normality test is low at short timescales, where the movement is ballistic. At longer

FIG. 5. Checking the applicability of the diffusive approximation. (a) Normalized autocorrelation function of steps for the MD ABF trajectory at various timescales (t equals 1 fs). (b) The average p-value in all bins versus the timescale. The error bars indicate on standard deviation of 10 independent trajectories.

J. Chem. Phys. 140, 084109 (2014)

timescales, the p-value is very close to 1, indicating that steps are very likely normally distributed (Fig. 5(b)), therefore we conclude that the system can be approximated by Brownian motion. The aim of this work is to obtain D at longer sampling timescales, where the dynamics is diffusive, not ballistic. Here the position dependence of Dij poses a challenge: as the sampling timescale is increased, Dij is averaged over larger and larger areas (due to the need to pad trajectories in bins to ensure unbiased sampling of large steps), so Dij in a particular bin becomes contaminated by the Di j of neighboring bins. Figure 6(a) demonstrates how the steps size, defined as four standard deviations of the steps distributions, increases with the timescale. If all steps started at the center of the bin, 4σ would be the size of the bin that contains 95% of all steps. We see that at long timescales, the steps sizes are large in comparison to the bin size chosen for analysis (36◦ , Fig. 6(a) horizontal gray line). If the process is diffusive and the diffusion is not too rapid, there should be a timescale window where the dynamics is no longer ballistic and averaging over neighboring bins is not yet excessive. In order to estimate the most optimal value of D for further kinetic studies and to minimize averaging over bins, we use the value of D∞ ij at “infinite” timescale as obtained from a timescale-dependent model for Brownian particles in a harmonic potential,29 given by Eq. (15). By using such a model, the human bias that could be introduced by choosing Dij at an arbitrary timescale manually is also eliminated. An example is presented in Fig. 6(b). First the Dψψ component (blue): the model nicely describes the ballistic behavior at short timescales (∼1-15 fs) and reasonably approximates the transient-cage effect (∼15-45 fs). A short plateau follows (∼45-70 fs) where a reasonable value of diffusion for further work can be chosen. At longer timescales (>70 fs), there is some discrepancy between the model and data. The “measured” diffusion increases due to averaging over neighboring bins, whereas the model predicts that the value from the short plateau continues to infinity, eliminating the contamination of averaging over neighboring bins. Sometimes, as in the case of Dϕϕ (green), the plateau is absent and the model averages diffusion values at the transient-cage region and at longer timescales, nevertheless the estimate of D∞ ϕϕ remains reasonable. In the framework of Kramers theory, kinetic constants are exponentially sensitive to the barrier height and only linearly to diffusion.17 Therefore, we conclude that the agreement between D∞ ij and diffusion data is acceptable and that values D∞ ij are suitable for further kinetic studies. Figure 6(c) shows that the timescale-dependent model nicely fits the expected plateau value (at 46 fs) of diffusion in all bins. We have also explored the effect of the Langevin thermostat39 on intrinsic diffusive dynamics, i.e., dynamics resulting from the coupling of fast degrees of freedom to the colvars. A Langevin thermostat introduces random forces and viscous drag on each atom to control temperature and could therefore affect the intrinsic diffusive dynamics. The damping parameter γ , also known as the friction parameter or coupling constant, represents the inverse decay time for inertia of atoms in the system. At higher γ values (smaller decay time), the dynamics of atoms becomes diffusive. Normally a value


084109-9


J. Chem. Phys. 140, 084109 (2014)

FIG. 7. Diffusion vs. timescale at different values of Langevin thermostat damping parameter (γ ) as obtained from ABF MD simulations of alanine dipeptide. Padding was used in the calculation of Dx˜ . The γ values are 1, 10, 100, and 1000 ps−1 (from red to orange lines, respectively). Damping values of 1 ps−1 does not significantly perturb the dynamics of the system. Error bars represent the standard deviation of 10 independent runs.

C. Advantage of diffusive approximation

The main application of determining diffusion coefficients is the lower computational cost. As the average dynamics of fast changing degrees of freedom are reflected in D, the degrees of freedom are drastically reduced. To illustrate: the alanine dipeptide system contains ∼ 2700 degrees of freedom, while the diffusive approximation of ϕ and ψ dynamics requires just 2 degrees of freedom. This results in approximately a 105 reduction in computational time. The trajectories of colvars generated by diffusive dynamics can be used for determining transition rates, transition probabilities, and lifetimes of metastable states.40

D. Comparison to other methods FIG. 6. Estimating diffusion coefficients at “infinite” timescale by a timescale dependent model for Brownian particles in a harmonic potential from MD ABF data. (a) Steps size, defined as four standard deviations of the steps distributions, in ϕ (green) and ψ (blue) directions vs. timescale. If all steps started at the center of the bin, 4σ would be the size of the bin that contains 95% of all steps. The width of the bins used during the analysis (36◦ ) is marked with a grey horizontal line. (b) Components of the diffusion tensor Dϕϕ , Dϕψ , and Dψψ (dashed-dotted green, orange and blue, respectively) are fitted by Eq. (14) (solid lines). The horizontal dashed lines represent the value of diffusion at infinity (D∞ ij ) as given by the model fit. The dashed vertical line denotes the timescale at which the short plateau region starts (46 fs). Diffusions in panel (d) (blue ellipses) are also displayed at 46 fs. The bin to which Dϕϕ , Dϕψ , and Dψψ belong is marked in panel (c) with a grey rectangle. Error bars represent the standard deviation of 10 independent runs. (d) Comparison of diffusion tensors given by the timescale model (black dashed ellipses) and diffusion values (blue ellipses) at the start of the plateau region (46 fs).

between one and ten inverse picoseconds is used. A damping of 10 ps−1 already has a noticeable effect on the diffusion; therefore, we recommend that researchers wishing to preserve realistic dynamics do not use damping parameters larger than 1 ps−1 (Fig. 7).

In this work, we present anisotropic 2D positiondependent example cases, develop the idea of padding trajectories in bins to improve sampling, explicitly treat the transition from ballistic to diffusive regime with a time dependent model and show that calculation of diffusion coefficients can be combined with biased MD. A few other studies for determining position-dependent diffusion coefficients have been published, but present only 1D examples, while overlooking the transition from ballistic to diffusive regime. There is also no freely available software for determining 1D/2D positiondependent diffusion. Hummer17 estimates the 1D position dependent diffusion coefficients and free energy using a Markovian rate equation and Bayesian inference. Very long term diffusive dynamics can be estimated, as the transitions between bins are analyzed instead of step sizes. Hummer reconstructs the 1D free energy and diffusion surface of the ψ dihedral of alanine dipeptide, which matches well with our results (see Fig. S3 of the supplementary material42 ); however, poor sampling in high energy regions is still noticed. In contrast to our method, the implementation and the generalization to higher dimensions are demanding.


084109-10


Comer et al.18 combine the Bayesian inference approach with ABF MD. They reconstruct a 1D diffusion surface of a water molecule, where the position dependence steams from a position-dependent thermostat. Comer et al. note that the sampling interval (termed frame interval in their paper) must be long enough that inertial effect are dispersed and short enough that the diffusion and energy change relatively little during one step. But the transition from the inertial to the diffusive regime is not explicitly treated. They conclude that ABF greatly improves sampling while not affecting the diffusive behavior. Micheletti et al.19 obtain D(x,τ ) in the presence of an external driving force, which can overcome sampling problems in high energy regions. A simple formula for D(x,τ ), expressed as an explicit average over the colvars trajectory, is given by their Eqs. (8) and (9). In the case of no external driving force, the expression simplifies to the Kramers-Moyal coefficients and is equivalent to our Eqs. (9) and (10). Micheletti et al. present a one dimensional test case of a model polypeptide chain in vacuum. They indicate that the longest timescale at which the step’s lengths still resemble a Gaussian distribution results is the best estimate of the diffusion coefficients. Since no explicit solvent is present in their study, the velocity autocorrelation follows a simple exponential decay and the treatment of oscillatory VAFs is not mentioned. Their model system uses very small integration steps, so diffusion steps are always much smaller than the size of the bin. Consequently, there is no need to pad the trajectories. Hegger and Stock20 use the Langevin equation to model the behavior of the slowest dihedral principle components. Their Eqs. (11) and (12) are analogous to Eqs. (9) and (10) in this paper. Hegger and Stock do not mention the possible timescale dependence of the diffusion coefficients. The diffusions are calculated by averaging over a “neighborhood” of points, similarly to the idea of padding, except that the minimal neighborhood size is not dynamically adjusted. Their sampling could be improved using biased MD. Hegger and Stock show that Langevin modeling can reproduce the density distribution of the principal dihedral components as well as the lifetimes of the metastable states as obtained by MD. Finally, we should mention that some processes are not diffusive. The diffusion length scale might be too large (i.e., the diffusion too rapid), which would limit resolution of position dependence. Memory effects, (i.e., random forces that are not normally distributed), can also occur, but this is not accounted for in this work. Therefore, these phenomena could potentially constrain the applicability of the presented approach. An interesting future direction would be to extend this work to analysis of anomalous diffusion,41 since anomalous diffusion is common in crowded biological systems, such as cells, gels, networks, etc. V. CONCLUSION

The computational costs of all atom MD can be reduced by studying long term diffusive dynamics of only the most important collective variables (colvars) in the system. Precise position-dependent diffusion coefficients are necessary for this approach. We have therefore developed a method

J. Chem. Phys. 140, 084109 (2014)

for determining multidimensional anisotropic position- and timescale-dependent diffusion coefficients by analysing the changes of colvars in an MD trajectory. First, we were able to successfully reconstruct the diffusive surface from a trajectory generated by a synthetic 2D test model. Second, we analyzed the dynamics of the ϕ and ψ dihedral angles of the alanine dipeptide (using both biased ABF MD and unbiased MD simulations). We successfully obtained the long-term diffusion coefficient by applying a timescale model for Brownian particles in a harmonic potential.29 The latter was needed since the dependence of diffusion coefficient on the sampling timescale revealed a transition from ballistic to diffusive regime due to Newton’s laws of motion governing MD. The applied ABF bias greatly enhanced sampling in high energy regions, while not affecting the diffusion coefficients. Using the diffusive model, a few hundred nanosecondlong ABF MD trajectories can be extended to millisecond lengths in the colvars subspace at a very small additional computational cost. From the extended trajectory kinetic data such as transition rates can be extracted. To increase the availability of our approach to the research community, we implemented R package availthe algorithm as an open source Mathematica able at https://github.com/lbf-ijs/DiffusiveDynamics (see the supplementary material42 ). The package is capable of obtaining and visualizing D for one and two dimensional cases. If the free energy is known, diffusive trajectories can be generated as well.

ACKNOWLEDGMENTS

The authors would like to thank Dr. Andrej Vilfan for helpful discussions and for reading the paper, William MacAskill for proofreading, and the Slovenian Research Agency (ARRS, Grant Nos. P1-0060, PR 03089-1, and PR 03090-1) for funding. 1 P.

L. Freddolino, A. S. Arkhipov, S. B. Larson, A. McPherson, and K. Schulten, Structure 14, 437 (2006). 2 A. Einstein, Ann. Phys. 324, 371 (1906). 3 P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990). 4 H. A. Kramers, Physica 7, 284 (1940). 5 B. A. Luty, R. C. Wade, J. D. Madura, M. E. Davis, J. M. Briggs, and J. A. McCammon, J. Phys. Chem. 97, 233 (1993). 6 D. I. Kopelevich, A. Z. Panagiotopoulos, and I. G. Kevrekidis, J. Chem. Phys. 122, 044908 (2005). 7 K. Molmer, J. Phys. B: At., Mol. Opt. Phys. 27, 1889 (1994). 8 E. Frey and K. Kroy, Ann. Phys. 14, 20 (2005). 9 A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. 99, 12562 (2002). 10 A. Barducci, M. Bonomi, and M. Parrinello, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 1, 826 (2011). 11 E. Darve, D. Rodríguez-Gómez, and A. Pohorille, J. Chem. Phys. 128, 144120 (2008). 12 J. Hénin, G. Fiorin, C. Chipot, and M. L. Klein, J. Chem. Theory Comput. 6, 35 (2010). 13 G. Hummer and A. Szabo, Proc. Natl. Acad. Sci. U.S.A. 98, 3658 (2001). 14 S. Park, F. Khalili-Araghi, E. Tajkhorshid, and K. Schulten, J. Chem. Phys. 119, 3559 (2003). 15 J. Kästner, Wiley Interdiscip. Rev.: Comput. Mol. Sci. 1, 932 (2011). 16 P. S. Grassia, E. J. Hinch, and L. C. Nitsche, J. Fluid Mech. 282, 373 (1995). 17 G. Hummer, New J. Phys. 7, 34 (2005). 18 J. Comer, C. Chipot, and F. D. González-Nilo, J. Chem. Theory Comput. 9, 876 (2013).


084109-11 19 C.


Micheletti, G. Bussi, and A. Laio, J. Chem. Phys. 129, 074105 (2008). 20 R. Hegger and G. Stock, J. Chem. Phys. 130, 034106 (2009). 21 P. E. Smith, J. Chem. Phys. 111, 5568 (1999). 22 J. Hermans, Proc. Natl. Acad. Sci. U.S.A. 108, 3095 (2011). 23 H. Jang and T. B. Woolf, J. Comput. Chem. 27, 1136 (2006). 24 D. S. Lemons, Am. J. Phys. 65, 1079 (1997). 25 R. Kubo, Rep. Prog. Phys. 29, 255 (1966). 26 D. L. Ermak and J. A. McCammon, J. Chem. Phys. 69, 1352 (1978). 27 T. W. Anderson and D. A. Darling, Ann. Math. Stat. 23, 193 (1952). 28 U. Balucani, J. P. Brodholt, and R. Vallauri, J. Phys.: Condens. Matter 8, 6139 (1996). 29 L. Glass and S. A. Rice, Phys. Rev. 176, 239 (1968). 30 D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, 2nd ed. (Academic Press, 2001). 31 W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph. 14, 33 (1996).

J. Chem. Phys. 140, 084109 (2014) 32 J.

C. Phillips, R. Braun, W. Wang, J. Gumbart, E. Tajkhorshid, E. Villa, C. Chipot, R. D. Skeel, L. Kalé, and K. Schulten, J. Comput. Chem. 26, 1781 (2005). 33 A. D. Mackerell, Jr., M. Feig, and C. L. Brooks III, J. Comput. Chem. 25, 1400 (2004). 34 A. D. Mackerell, J. Comput. Chem. 25, 1584 (2004). 35 T. Li and M. G. Raizen, Ann. Phys. 525, 281 (2013). 36 T. Li, S. Kheifets, D. Medellin, and M. G. Raizen, Science 328, 1673 (2010). 37 R. Huang, I. Chavez, K. M. Taute, B. Luki´ c, S. Jeney, M. G. Raizen, and E.-L. Florin, Nat. Phys. 7, 576 (2011). 38 A. Morita and B. Bagchi, J. Chem. Phys. 110, 8643 (1999). 39 D. Quigley and M. I. J. Probert, J. Chem. Phys. 120, 11432 (2004). 40 W.-N. Du, K. A. Marino, and P. G. Bolhuis, J. Chem. Phys. 135, 145102 (2011). 41 I. M. Sokolov, Soft Matter 8, 9043 (2012). 42 See supplementary material at http://dx.doi.org/10.1063/1.4866448 for three additional figures and a demonstration of the DiffusiveDynamics package.


Diffusion of molecules in the bulk of a low density amorphous ice from molecular dynamics simulations.

Molecular dynamics simulations in biology.

Floating orbital molecular dynamics simulations.

Nucleation rate analysis of methane hydrate from molecular dynamics simulations.

Molecular dynamics simulations: from structure function relationships to drug discovery.

Thermodynamics of writhe in DNA minicircles from molecular dynamics simulations.

Distance-Based Configurational Entropy of Proteins from Molecular Dynamics Simulations.

Molecular Dynamics Simulations Reveal that Water Diffusion between Graphene Oxide Layers is Slow.

Effect of correlation on viscosity and diffusion in molecular-dynamics simulations.

Microsecond molecular dynamics simulations of lipid mixing.

Molecular dynamics simulations of classical stopping power.

Extracting the diffusion tensor from molecular dynamics simulation with Milestoning.

Molecular dynamics simulations: advances and applications.

AceCloud: Molecular Dynamics Simulations in the Cloud.

Molecular dynamics simulations of large macromolecular complexes.

Accelerated molecular dynamics simulations of protein folding.

Molecular dynamics simulations of helix denaturation.

New ways to boost molecular dynamics simulations.

Molecular dynamics simulations of glycoproteins using CHARMM.

Simulations of pattern dynamics for reaction-diffusion systems via SIMULINK.

Molecular dynamics simulations: insight into molecular phenomena at interfaces.

Dynamics of lipids, cholesterol, and transmembrane α-helices from microsecond molecular dynamics simulations.

Structure and dynamics of [PF6][P(1,2,2,4)] from molecular dynamics simulations.

mineral interfaces from ab initio molecular dynamics simulations.