Fully Anisotropic Rotational Diffusion Tensor from Molecular Dynamics

Subscriber access provided by Universitaetsbibliothek | Johann Christian Senckenberg

Article

Fully Anisotropic Rotational Diffusion Tensor from Molecular Dynamics Simulations Max Linke, Juergen Koefinger, and Gerhard Hummer J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b11988 • Publication Date (Web): 30 Jan 2018 Downloaded from http://pubs.acs.org on January 31, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Fully Anisotropic Rotational Diffusion Tensor from Molecular Dynamics Simulations Max Linke,† Jürgen Köfinger,† and Gerhard Hummer∗,† †Max Planck Institute of Biophysics, Max-von-Laue-Str. 3, 60438 Frankfurt am Main, Germany ‡Department of Physics, Goethe University Frankfurt, Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany E-mail: [email protected]

1

ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We present a method to calculate the fully anisotropic rotational diffusion tensor from molecular dynamics simulations. Our approach is based on fitting the timedependent covariance matrix of the quaternions that describe the rigid-body rotational dynamics. Explicit analytical expressions have been derived for the covariances by Favro, which are valid irrespective of the degree of anisotropy. We use these expressions to determine an optimal rotational diffusion tensor from trajectory data. The molecular structures are aligned against a reference by optimal rigid-body superposition. The quaternion covariances can then be obtained directly from the rotation matrices used in the alignment. The rotational diffusion tensor is determined by a fit to the time-dependent quaternion covariances, or directly by Laplace transformation and matrix diagonalization. To quantify uncertainties in the fit, we derive analytical expressions and compare them with the results of Brownian dynamics simulations of anisotropic rotational diffusion. We apply the method to microsecond long trajectories of the Dickerson-Drew B-DNA dodecamer and of horse heart myoglobin. The anisotropic rotational diffusion tensors calculated from simulations agree well with predictions from hydrodynamics.

INTRODUCTION The rotational dynamics of macromolecules in solution is a key factor in the kinetics of macromolecular binding and assembly. 1–5 Rotation also features prominently in theoretical descriptions of many experiments, such as light scattering, 6 nuclear magnetic resonance (NMR), 7–16 fluorescence anisotropy decay and fluorescence resonance energy transfer (FRET), 17–20 or dielectric spectroscopy. 21 The theoretical descriptions of the rotational dynamics used in the analysis of these experiments typically build on the model of rotational diffusion, dating back to the pioneering work of Debye 22 and Perrin. 23,24 Perrin’s hydrodynamic formulation of rotational diffusion generalizes the Stokes-Einstein formula for translational diffusion, and

2


Page 2 of 35

Page 3 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


makes it possible to estimate rotational diffusion tensors of arbitrarily shaped molecules in a neat solvent. 25 However, despite its remarkable success, at the molecular level the hydrodynamic formulation remains an approximation that is not easily generalized from neat to complex fluids. Indeed, recent NMR experiments on high density solutions of proteins 26,27 show a larger rotational diffusion anisotropy compared to proteins in dilute solutions, consistent with the behavior expected from steric clashes seen in molecular simulations. 28–30 Molecular dynamics (MD) provide a powerful tool to probe the rotational dynamics in simple and complex fluids. In practice, determining the rotational diffusion tensor of a molecule or molecular assembly from MD trajectories is far from trivial. Current methods to calculate rotational diffusion tensors using the second order rotational correlation function have focused on isotropic macromolecules or macromolecules with small anisotropy. 31,32 However, Wang and Case’s method 32 has been extended to fully asymmetric rotational diffusion by Roe and Cheatham. 33 The focus of recent quaternion-based algorithms has also been on small anisotropy or on estimates of the full anisotropic diffusion tensor from the rotational dynamics at short times only. 34,35 Here, we present an algorithm to extract from an MD simulation a fully anisotropic rotational diffusion tensor that accounts for the rotational dynamics over the entire time domain. We describe anisotropic rigid-body rotational diffusion using quaternions. Quaternions often lead to compact mathematical formulations of problems associated with three-dimensional rotations, e.g., in an elegant and powerful algorithm for the optimal superposition of macromolecular structures. 36 Favro 37 derived explicit analytical expressions for the time-dependent covariances of the Cayley-Klein parameters in the principal coordinate system (PCS). Because the Cayley-Klein parameters are closely related to quaternions, Favro’s expressions immediately give us the covariances of the quaternion coefficients in the PCS. By transforming the quaternion covariance matrix into an arbitrary body-fixed coordinate system, we arrive at a formulation that allows us to calculate rotational diffusion tensors from molec-

3



ular simulations using common root-mean-square distance (RMSD) fitting algorithms and standard linear-algebra algorithms. To assess uncertainties in the estimated diffusion tensor, we perform Brownian dynamics simulations of anisotropic rotational diffusion. We apply our method to ten 1-µs long trajectories of the Dickerson-Drew dodecamer, a B-DNA double helix, and five 800 ns long trajectories of myoglobin, and compare the resulting rotational diffusion tensors to estimates from hydrodynamic theory and to experimental results. For these relatively rigid molecules in neat aqueous solution, we expect that the hydrodynamic formulation is essentially exact, which allows us to test our algorithm to calculate rotational diffusion tensors from MD simulation trajectories.

THEORY Rotational Diffusion in the Principal Coordinate System. To describe the current orientation of a rigid body we use quaternions or, more precisely, column vectors of the socalled Euler parameters, 38,39 q = (q0 , q1 , q2 , q3 )T , with q02 + q12 + q22 + q32 = 1 and T indicating the matrix transpose. Quaternions have proven useful in computer modeling of rigid-body motions as an alternative to Euler angles. The quaternion coefficients are related to the rotation axes ~v and angle φ through q = (cos(φ/2), sin(φ/2)~v ). A quaternion q = (1, , 0, 0) with → 0 thus corresponds to an infinitesimal rotation about the x axis by an angle /2. See eq 7 below for a representation of the 3 × 3 rotation matrix in terms of the coefficients of a “quaternion” u = (u0 , u1 , u2 , u3 )T . In this work, we will use results of Favro 37 who derived explicit expression for the rotational dynamics of rigid bodies using Cayley-Klein parameters. 38 Because of the close relationship of the Cayley-Klein parameters and quaternions, Favro’s results are readily adapted to quaternions. The quaternion q(t) describes the orientation of the rigid body in the PCS as a function of time t. In the PCS, the rotational diffusion tensor is diagonal, D = diag(D1 , D2 , D3 ). For rotational diffusion of a rigid body that is aligned with the PCS at time 0, q(0) = (1, 0, 0, 0)T ,

4


Page 4 of 35

Page 5 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Favro 37 derived explicit expressions for the average covariance coefficients hqi (t)qj (t)i as a function of time, 1 4 1 hq12 (t)i = 4 1 hq22 (t)i = 4 1 hq32 (t)i = 4

hq02 (t)i =

1 + e−3Dt eD1 t + eD2 t + eD3 t

(1)

1 + e−3Dt eD1 t − eD3 t − eD2 t

(2)

1 + e−3Dt eD2 t − eD1 t − eD3 t

(3)

1 + e−3Dt eD3 t − eD2 t − eD1 t

(4)

hqi (t)qj (t)i = 0 for i 6= j,

(5)

where D = TrD/3 = (D1 + D2 + D3 )/3. The average h· · · i is over repeated initializations of the stochastic trajectories or, equivalently, different starting points on a long equilibrium trajectory, appropriately rotated into the PCS.

Rotational Diffusion in the Reference Coordinate System. In practice, we do not a priori know the PCS in which D is diagonal. Therefore, the time-dependent covariance matrix with elements Cij (t) = [hu(t)uT(t)i]ij = hui (t)uj (t)i for quaternions u(t) describing the orientation in an arbitrarily chosen reference coordinate system (RCS) is in general not diagonal. However, an appropriately chosen rigid body rotation R transforms the system from the RCS into the PCS, so that

RC(t)RT = hdiag( q12 (t), q22 (t), q32 (t) )i ,

(6)

where RT is the matrix of orthonormal eigenvectors of C(t). For ideal Brownian rotations, the eigenvalues follow eqs 2 to 4, and R will be independent of time t. Since RT rotates a

5



Page 6 of 35

vector from the PCS into the RCS, 



 1    T  RT   0 ,R   0







 0   0     1  , and RT  0       0 1

     

are the principal axes in the RCS. Therefore, the row vectors in the matrix R of eigenvectors of C(t) give the principal axes.

Calculating Orientational Covariance Matrix from Simulation Trajectories. To calculate hu(t)uT(t)i from simulations, it is convenient to express the products of quaternion coefficients ui (t)uj (t) in terms of rotation matrix elements. Common tools provide a rotation matrix for an RMSD superposition of two molecular structures. Let U be the RMSD superposition matrix that rotates the macromolecule of interest into the (predefined) reference orientation. Note that this RCS will in general differ from the initially unknown PCS. The matrix U is the output of an RMSD alignment, a task that can itself be accomplished using quaternions. 36,40,41 We can represent this rotation in terms of the coefficients of the quaternion u = (u0 , u1 , u2 , u3 )T as 38 

2



2

 1 − 2(u2 + u3 ) 2(u0 u3 + u1 u2 ) 2(u1 u3 − u0 u2 )  2 2 U =  2(u1 u2 − u0 u3 ) 1 − 2(u1 + u3 ) 2(u0 u1 + u2 u3 )  2(u0 u2 + u1 u3 ) 2(u2 u3 − u0 u1 ) 1 − 2(u1 2 + u2 2 )

6


  .  

(7)

Page 7 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We now invert these relations to obtain explicit expressions for all pair products ui uj of the covariance matrix coefficients in terms of the coefficients Uij of the rotation matrix,

u0 u0 = (1 + U11 + U22 + U33 )/4

(8)

u0 u1 = (U23 − U32 )/4 u0 u2 = (U31 − U13 )/4 u0 u3 = (U12 − U21 )/4 u1 u1 = (1 + U11 − U22 − U33 )/4 u1 u2 = (U12 + U21 )/4 u1 u3 = (U13 + U31 )/4 u2 u2 = (1 − U11 + U22 − U33 )/4 u2 u3 = (U23 + U32 )/4 u3 u3 = (1 − U11 − U22 + U33 )/4.

The remaining coefficients in the uuT matrix follow from its symmetry, ui uj = uj ui . We use these relations to calculate the average covariance matrix hu(t)uT(t)i from an MD trajectory of a macromolecule. To evaluate hu(t)uT(t)i at time t from a trajectory segment starting at time τ we use the 3 × 3 rotation matrix S(τ + t) that rotates the macromolecule structure at time τ + t into optimal superposition with the reference (i.e., into the RCS). S(τ + t) can be broken up into a rotation S(τ ) to move to the RCS at time τ and another rotation U (t; τ ) to then establish optimal superposition from there, S(τ + t) = U (t; τ )S(τ ). Because SS T = S TS = I, we can solve for the rotation matrix U (t; τ ) describing the orientation of the molecule at time τ + t with t ≥ 0, starting in the RCS at time τ ,

U (t; τ ) = S(τ + t)S T(τ ).

7


(9)


Page 8 of 35

The sums of coefficients Uij (t; τ ) in eq 8 provide the coefficients ui (t; τ )uj (t; τ ) at time t for the trajectory segment starting at time τ . The average time-dependent covariance matrix is obtained by repeating this procedure for all Nτ starting times τ and averaging over them,

hui (t)uj (t)i =

1 X ui (t; τ )uj (t; τ ). Nτ τ

(10)

Estimating the Diffusion Tensor from Trajectory Data. To determine the optimal rotational diffusion tensor from the covariance matrix Cij (t) = hui (t)uj (t)i (i = 1, 2, 3) we use a least-squares fit to the covariance matrix. The statistical variances s2ij (t) of the covariance matrix coefficients, hui (t)uj (t)i, are given by s2ij (t) = hu2i (t)u2j (t)i − hu2i (t)ihu2j (t)i.

(11)

Estimates of the fourth moments hu2i (t)u2j (t)i from simulations are problematic because they are themselves subject to large uncertainties. Therefore, we instead use analytic expressions hu2i (t)u2j (t)iid for the fourth moments that are valid for ideal rotational diffusion with given diffusion coefficients and PCS. For the fourth-order correlations, Favro 37 obtained

hq14 (t)i =

1 3 1 + e−3Dt eD1 t − eD2 t − eD3 t 8 2 1 + e−3Dt e−3D1 t − e−3D2 t − e−3D3 t 2 ! + e−6Dt cosh(2t∆)

8


(12)

Page 9 of 35 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and

hq12 (t)q22 (t)i =

1 1 1 −3Dt D3 t − e e − e−3D3 t 8 3 2 + e−6Dt ∆−1 (D3 − D) sinh(2t∆) ! 1 − cosh(2t∆) , 3

(13)

where

∆ = D12 + D22 + D32 − D1 D2 − D1 D3 − D2 D3

1/2

.

(14)

The other hqi2 (t)qj2 (t)i can be found by permutations of the indices. Averages over products of odd powers of the qi vanish because of symmetry. We transform the hqi2 (t)qj2 (t)i from the PCS into the RCS with R in eq 6, resulting in

hu2i (t)u2j (t)iid

=

3 X

2 2 hqi2 (t)qj2 (t)iRki Rkj +

k=1 3 X

2 2 hqi2 (t)qj2 (t)i Rmi Rni +

m

Fully Anisotropic Rotational Diffusion Tensor from Molecular Dynamics

Recommend Documents