This section is dealing with specific types of analysis performed by MDANSE. If you are not sure where these fit into the general workflow of data analysis, please read Workflow of the Analysis.

10. Analysis: Dynamics

This section contains the following Plugins:

10.1. Angular Correlation

  • Available for trajectories only

Theory and implementation

The angular correlation analysis computes the autocorrelation of a set of vectors describing the extent of a molecule in three orthogonal directions. This kind of analysis can be useful when trying to highlight the fact that a molecule is constrained in a given direction.

For a given triplet of non-colinear atoms g=(a1,a2,a3), one can derive an orthonormal set of three vectors v1, v2, v3 using the following scheme:

  • v_{1} = \frac{n_{1} + n_{2}}{\left| \left| {n_{1} + n_{2}} \right| \right|} where n1 and n2 are respectively the normalized vectors along (a1,a2) and (a1,a3) directions.

  • v2 is defined as the clockwise normal vector orthogonal to v1 that belongs to the plane defined by a1, a2 and a3 atoms

  • {\overrightarrow{v_{3}} = \overrightarrow{v_{1}}}\times\overrightarrow{v_{2}}

Thus, one can define the following autocorrelation functions for the vectors v1, v2 and v3 defined on triplet t :

(10.1){AC_{g,i}{(t) = \left\langle {v_{t,i}(0)\cdot v_{t,i}(t)} \right\rangle},{i = 1,2,3}}

And the angular correlation averaged over all triplets is:

(10.2){AC_{i}{(t) = {\sum\limits_{g = 1}^{N_{\mathit{triplets}}}{AC_{g,i}(t)}}},{i = 1,2,3}}

where Ntriplets is the number of selected triplets.

GUI

../_images/10000001000003220000021BA346CFB73717C1BE.png

Parameters:

10.2. Density Of States

Theory and implementation

MDANSE calculates the power spectrum of the VACF, which in case of the mass-weighted VACF defines the phonon discrete DOS, (see the section on Velocity AutoCorrelation Function) defined as:

(10.3){\mathit{DOS}\left( {n\cdot\mathit{\Delta\nu}} \right)\doteq{\sum\limits_{\alpha}\omega_{\alpha}}{\overset{\sim}{C}}_{\mathit{vv};\mathit{\alpha\alpha}}\left( {n\cdot\mathit{\Delta\nu}} \right),{n = 0}\ldots{N_{t} - 1.}}

Nt is the total number of time steps and {\mathit{\Delta\nu} = 1}\text{/}\left( {2N_{t}\Delta t} \right) is the frequency step. {\mathit{DOS}\left( {n\cdot\mathit{\Delta\nu}} \right)} can be computed either for the isotropic case or with respect to a user-defined axis. The spectrum {\mathit{DOS}\left( {n\cdot\Delta\nu} \right)} is computed from the unnormalized VACF, such that DOS(0) gives an approximate value for the diffusion constant {D = {\sum\limits_{\alpha}D_{\alpha}}} (see Eqs. (10.12) and (10.13)).

{\mathit{DOS}\left( {n\cdot\Delta\nu} \right)} is smoothed by applying a Gaussian window in the time domain [Ref10] (see the section The FCA algorithm). Its width in the time domain is {\sigma_{t} = \alpha}\text{/}T , where T is the length of the simulation. We remark that the diffusion constant obtained from DOS is biased due to the spectral smoothing procedure since the VACF is weighted by this window Gaussian function. MDANSE computes the density of states starting from both atomic velocities and atomic coordinates. In this case the velocities are computed by numerical differentiation of the coordinate trajectories correcting first for possible jumps due to periodic boundary conditions.

GUI

  • available for trajectories only

image14 image15

10.3. General AutoCorrelation Function

  • available for trajectories only

image16 image17

Format: drop-down

Default: configuration

Description: determines the variable for which the autocorrelation function is calculated. Therefore, if the selected variable is ‘configuration’, essentially position autocorrelation function is calculated.

10.4. Mean Square Displacement

Theory and implementation

Molecules in liquids and gases do not stay in the same place but move constantly. This process is called diffusion and it happens quite naturally in fluids at equilibrium. During this process, the motion of an individual molecule does not follow a simple path. As it travels, the molecule undergoes some collisions with other molecules which prevent it from following a straight line. If the path is examined in close detail, it will be seen to be a good approximation to a random walk. Mathematically, a random walk is a series of steps where each step is taken in a completely random direction from the one before. This kind of path was famously analysed by Albert Einstein in a study of Brownian motion. He showed that the Mean-Square Displacement (MSD) of a particle following a random walk is proportional to the time elapsed. This relationship can be written as

(10.4){\left\langle r^{2} \right\rangle = 6}{\mathit{Dt} + C}

where \langle r^{2}\rangle is the MSD and t is the time. D and C are constants. The constant D defines the so-called diffusion coefficient.

The Figure Fig. 10.1 shows an example of an MSD analysis performed on a water box of 768 water molecules. To get the diffusion coefficient out of this plot, the slope of the linear part of the plot should be calculated.

../_images/10000000000001BC00000163C18A769B32940652.png

Fig. 10.1 MSD calculated for a 100 ps MD simulation of 256 water molecules using NPT condition at 1 bar and 300 K.

Defining

(10.5){d_{\alpha}\left( {t,t_{0}} \right)\doteq R_{\alpha}{\left( {t_{0} + t} \right) - R_{\alpha}}\left( t_{0} \right),}

the MSD of particle \alpha can be defined as:

(10.6)\mathrm{\Delta}_{\alpha}^{2}{(t) = \left\langle {d_{\alpha}^{2}\left( {t,t_{0}} \right)} \right\rangle_{t_{0}}}

where R_{\alpha}(t_0) and R_{\alpha}(t_0 + t) are respectively the position of particle \alpha at times t_0 and t_0 + t. One can introduce an MSD with respect to a given axis n:

(10.7){\mathrm{\Delta}_{\alpha}^{2}\left( {t,t_{0};n} \right)\doteq\left\langle {d_{\alpha}^{2}\left( {t,\tau;n} \right)} \right\rangle_{t_{0}}}

with

(10.8){d_{\alpha}^{}\left( {t,\tau;n} \right)\doteq n\bullet d_{\alpha}^{}\left( {t,t_{0}} \right).}

The calculation of MSD is the standard way to obtain diffusion coefficients from Molecular Dynamics (MD) simulations. Assuming Einstein-diffusion in the long time limit one has for isotropic systems

(10.9){D_{\alpha} = {\lim\limits_{t\rightarrow\infty}{\frac{1}{6t}\mathrm{\Delta}_{\alpha}^{2}(t)}}}.

There exists also a well-known relation between the MSD and the velocity autocorrelation function. Writing

(10.10){d_{\alpha}{(t) = {\int\limits_{0}^{t}{\mathit{d\tau}v_{\alpha}(\tau)}}}}

in Eq. (10.6) one can show (see e.g. [Ref11]) that

(10.11){\mathrm{\Delta}_{\alpha}^{2}{(t) = 6}{\int\limits_{0}^{t}{\mathit{d\tau}\left( {t - \tau} \right)C_{\mathit{\upsilon\upsilon};\mathit{\alpha\alpha}}(t)}}.}

Using now the definition (10.9) of the diffusion coefficient one obtains the relation

(10.12){{D_{\alpha} = {\int\limits_{0}^{t}{\mathit{d\tau}C_{\mathit{\upsilon\upsilon};\mathit{\alpha\alpha}}(t)}}}.}

With Eq. (10.42) this can also be written as

(10.13){{D_{\alpha} = \pi}{\overset{\sim}{C}}_{\mathit{\upsilon\upsilon};\mathit{\alpha\alpha}}(0).}

Computationally, the MSD is calculated using the Fast Correlation Algorithm (FCA) [Ref12]. In this framework, in the discrete case, the mean-square displacement of a particle is given by

(10.14){\mathrm{\Delta}^{2}{(m) = \frac{1}{N_{t} - m}}{\sum\limits_{k = 0}^{N_{t} - m - 1}\left\lbrack {r{\left( {k + m} \right) - r}(k)} \right\rbrack^{2}},{m = 0.}..{N_{t} - 1}}

where r(k) is the particle trajectory and N_{t} is the number of frames of the trajectory. We define now the auxiliary function

(10.15){S(m)\doteq{\sum\limits_{k = 0}^{N_{t} - m - 1}\left\lbrack {r{\left( {k + m} \right) - r}(k)} \right\rbrack^{2}},{m = 0}...N{t - 1},}

which is split as follows:

(10.16){S{(m) = S_{\mathit{AA} + \mathit{BB}}}{(m) - 2}S_{\mathit{AB}}(m),}

(10.17){S_{\mathit{AA} + \mathit{BB}}{(m) = \sum\limits_{k = 0}^{N_{t} - m - 1}}\left\lbrack {r^{2}{\left( {k + m} \right) + r^{2}}(k)} \right\rbrack,}

(10.18){S_{\mathit{AB}}{(m) = {\sum\limits_{k = 0}^{N_{t} - m - 1}{r(k)}}}\cdot r\left( {k + m} \right).}

The function SAB(m) can be computed using the FCA method described in the section The FCA algorithm. For SAA+BB(m) the following recursion relation holds:

(10.19){S_{\mathit{AA} + \mathit{BB}}{(m) = S_{\mathit{AA} + \mathit{BB}}}{\left( {m - 1} \right) - r^{2}}{\left( {m - 1} \right) - r^{2}}\left( {N_{t} - m} \right),}

(10.20){S_{\mathit{AA} + \mathit{BB}}{(0) = {\sum\limits_{k = 0}^{N_{t} - 1}{r^{2}(k)}}}.}

This allows one to construct the following efficient scheme for the computation of the MSD:

  1. Compute

    (10.21){\mathit{DSQ}{{(k)} = r}2{(k)},{k = 0}...N{t - 1}}

    ;

    (10.22){\mathit{DSQ}{{({- 1})} = \mathit{DSQ}}{{({Nt})} = 0}}

    .

  2. Compute

    (10.23){\mathit{SUMSQ} = 2}\cdot{\sum\limits_{k = 0}^{N_{t} - 1}{\mathit{DSQ}(k)}}

  3. Compute SAB(m) using the Fast Fourier Transform (FFT) method.

  4. Compute MSD(m) in the following loop:

(10.24){\mathit{SUMSQ}\leftarrow{\mathit{SUMSQ} - \mathit{DSQ}}{{({m - 1})} - \mathit{DSQ}}{({N{t - m}})}}

(10.25){\mathit{MSD}{(m)}\leftarrow{} (SUMSQ - 2 S_{AB}(m)/(N_t - m)}

m running from 0 to N_t - 1

It should be noted that the efficiency of this algorithm is the same as for the FCA computation

of time correlation functions since the number of operations in step (1), (2), and (4) grows

linearly with N_t.

GUI

image18 image19

10.5. Order Parameter

Theory and implementation

Adequate and accurate cross comparison of the NMR and MD simulation data is of crucial importance in versatile studies conformational dynamics of proteins. NMR relaxation spectroscopy has proven to be a unique approach for a site-specific investigation of both global tumbling and internal motions of proteins. The molecular motions modulate the magnetic interactions between the nuclear spins and lead for each nuclear spin to a relaxation behaviour which reflects its environment. Since its first applications to the study of protein dynamics, a wide variety of experiments has been proposed to investigate backbone as well as side chain dynamics. Among them, the heteronuclear relaxation measurement of amide backbone 15N nuclei is one of the most widespread techniques. The relationship between microscopic motions and measured spin relaxation rates is given by Redfield’s theory [Ref13]. Under the hypothesis that 15N relaxation occurs through dipole-dipole interactions with the directly bonded 1H atom and chemical shift anisotropy (CSA), and assuming that the tensor describing the CSA is axially symmetric with its axis parallel to the N-H bond, the relaxation rates of the 15N nuclei are determined by a time correlation function,

(10.26){C_{\mathit{ii}}{(t) = \left\langle {P_{2}\left( {\mu_{i}(0)\cdot\mu_{i}(t)} \right)} \right\rangle}}

which describes the dynamics of a unit vector μi(t) pointing along the 15N-1H bond of the residue i in the laboratory frame. Here P2(.) is the second order Legendre polynomial. The Redfield theory shows that relaxation measurements probe the relaxation dynamics of a selected nuclear spin only at a few frequencies. Moreover, only a limited number of independent observables are accessible. Hence, to relate relaxation data to protein dynamics one has to postulate either a dynamical model for molecular motions or a functional form for Cii(t), yet depending on a limited number of adjustable parameters. Usually, the tumbling motion of proteins in solution is assumed isotropic and uncorrelated with the internal motions, such that:

(10.27){C_{\mathit{ii}}{(t) = C^{G}}(t)\cdot C_{\mathit{ii}}^{I}(t)}

where CG(t) and

(10.28){C_{\mathit{ii}}^{I}(t)}

denote the global and the internal time correlation function, respectively. Within the so-called model free approach [Ref14], [Ref15] the internal correlation function is modelled by an exponential,

(10.29){C_{\mathit{ii}}^{I}{(t) = {S_{i}^{2} + \left( {1 - S_{i}^{2}} \right)}}\exp\left( \frac{- t}{\tau_{\mathit{eff},i}} \right)}

Here the asymptotic value

(10.30){S_{i}^{2} = C_{\mathit{ii}}}\left( {+ \infty} \right)

is the so-called generalized order parameter, which indicates the degree of spatial restriction of the internal motions of a bond vector, while the characteristic time

(10.31){tau_{\mathit{eff},i}}

is an effective correlation time, setting the time scale of the internal relaxation processes.

(10.32){S_{i}^{2}}

can adopt values ranging from 0 (completely disordered) to 1 (fully ordered). So,

(10.33){S_{i}^{2}}

is the appropriate indicator of protein backbone motions in computationally feasible timescales as it describes the spatial aspects of the reorientational motion of N-H peptidic bonds vector.

When performing Order Parameter analysis, MDANSE computes for each residue i both

(10.34){C_{\mathit{ii}}(t)}

and

(10.35){S_{i}^{2}}

. It also computes a correlation function averaged over all the selected bonds defined as:

(10.36){C^{I}{(t) = {\sum\limits_{i = 1}^{N_{\mathit{bonds}}}{C_{\mathit{ii}}^{I}(t)}}}}

where Nbonds is the number of selected bonds for the analysis.

GUI

  • available for trajectories only

../_images/10000001000003220000027563B8CBFF70E2089C.png

Format: int or float

Default: 0

Description: <insert>

    • y-component

Format: int or float

Default: 0

Description: <insert>

    • z-component

Format: int or float

Default: 1

Description: <insert>

10.6. Position AutoCorrelation Function

  • available for trajectories only

image20 image21

10.7. Velocity AutoCorrelation Function

Theory and implementation

The Velocity AutoCorrelation Function (VACF) is another interesting property describing the dynamics of a molecular system. Indeed, it reveals the underlying nature of the forces acting on the system.

In a molecular system that would be made of non-interacting particles, the velocities would be constant at any time triggering the VACF to be a constant value. Now, if we think about a system with small interactions such as in a gas-phase, the magnitude and direction of the velocity of a particle will change gradually over time due to its collision with the other particles of the molecular system. In such a system, the VACF will be represented by a decaying exponential.

In the case of solid phase, the interactions are much stronger and, as a results, the atoms are bound to a given position from which they will move backwards and forwards oscillating between positive and negative values of their velocity. The oscillations will not be of equal magnitude however, but will decay in time, because there are still perturbative forces acting on the atoms to disrupt the perfection of their oscillatory motion. So, in that case the VACF will look like a damped harmonic motion.

Finally, in the case of liquid phase, the atoms have more freedom than in solid phase and because of the diffusion process, the oscillatory motion seen in solid phase will be cancelled quite rapidly depending on the density of the system. So, the VACF will just have one very damped oscillation before decaying to zero. This decaying time can be considered as the average time for a collision between two atoms to occur before they diffuse away.

Mathematically, the VACF of atom \alpha in an atomic or molecular system is usually defined as

(10.37){C_{\mathit{vv};\mathit{\alpha\alpha}}(t)\doteq\frac{1}{3}\left\langle {v_{\alpha}\left( t_{0} \right)\cdot v_{\alpha}\left( {t_{0} + t} \right)} \right\rangle_{t_{0}}.}

In some cases, e.g. for non-isotropic systems, it is useful to define VACF along a given axis,

(10.38){C_{\mathit{vv};\mathit{\alpha\alpha}}\left( {t;n} \right)\doteq\left\langle {v_{\alpha}\left( {t_{0};n} \right)v_{\alpha}\left( {{t_{0} + t};n} \right)} \right\rangle_{t_{0}},}

where v_{\alpha}(t; n) is given by

(10.39){v_{\alpha}\left( {t;n} \right)\doteq n\cdot v_{\alpha}(t).}

The vector n is a unit vector defining a space-fixed axis.

The VACF of the particles in a many body system can be related to the incoherent dynamic structure factor by the relation:

(10.40){\lim\limits_{q\rightarrow 0}\frac{\omega^{2}}{q^{2}}S{\left( {q,\omega} \right) = G}(\omega),}

where G(\omega) is the Density Of States (DOS). For an isotropic system it reads

(10.41){G{(\omega) = {\sum\limits_{\alpha}{b_{\alpha,\mathit{inc}}^{2}{\overset{\sim}{C}}_{\mathit{vv};\mathit{\alpha\alpha}}(\omega)}}},}

(10.42){{\overset{\sim}{C}}_{\mathit{vv};\mathit{\alpha\alpha}}{(\omega) = \frac{1}{2\pi}}{\int\limits_{- \infty}^{+ \infty}\mathit{dt}}\exp\left\lbrack {{- i}\omega t} \right\rbrack C_{\mathit{vv};\mathit{\alpha\alpha}}(t).}

For non-isotropic systems, relation (10.40) holds if the DOS is computed from the atomic velocity autocorrelation functions

(10.43){C_{\mathit{vv};\mathit{\alpha\alpha}}\left( {t;n_{q}} \right)}

, where n_q is the unit vector in the direction of q.

GUI

  • available for trajectories only

image22 image23