4.8.3. Water dynamics analysis — MDAnalysis.analysis.waterdynamics

Author:Alejandro Bernardin
Year:2014-2015
Copyright:GNU Public License v3

New in version 0.11.0.

This module provides functions to analize water dynamics trajectories and water interactions with other molecules. The functions in this module are: water orientational relaxation (WOR) [Yeh1999], hydrogen bond lifetimes (HBL) [Rapaport1983], angular distribution (AD) [Grigera1995], mean square displacement (MSD) [Brodka1994] and survival probability (SP) [Liu2004].

For more information about this type of analysis please refer to [Araya-Secchi2014] (water in a protein cavity) and [Milischuk2011] (water in a nanopore).

References

[Rapaport1983](1, 2) D.C. Rapaport (1983): Hydrogen bonds in water, Molecular Physics: An International Journal at the Interface Between Chemistry and Physics, 50:5, 1151-1162.
[Yeh1999]Yu-ling Yeh and Chung-Yuan Mou (1999). Orientational Relaxation Dynamics of Liquid Water Studied by Molecular Dynamics Simulation, J. Phys. Chem. B 1999, 103, 3699-3705.
[Grigera1995]Raul Grigera, Susana G. Kalko and Jorge Fischbarg (1995). Wall-Water Interface. A Molecular Dynamics Study, Langmuir 1996,12,154-158
[Liu2004]Pu Liu, Edward Harder, and B. J. Berne (2004).On the Calculation of Diffusion Coefficients in Confined Fluids and Interfaces with an Application to the Liquid-Vapor Interface of Water, J. Phys. Chem. B 2004, 108, 6595-6602.
[Brodka1994]Aleksander Brodka (1994). Diffusion in restricted volume, Molecular Physics, 1994, Vol. 82, No. 5, 1075-1078.
[Araya-Secchi2014]Araya-Secchi, R., Tomas Perez-Acle, Seung-gu Kang, Tien Huynh, Alejandro Bernardin, Yerko Escalona, Jose-Antonio Garate, Agustin D. Martinez, Isaac E. Garcia, Juan C. Saez, Ruhong Zhou (2014). Characterization of a novel water pocket inside the human Cx26 hemichannel structure. Biophysical journal, 107(3), 599-612.
[Milischuk2011]Anatoli A. Milischuk and Branka M. Ladanyi. Structure and dynamics of water confined in silica nanopores. J. Chem. Phys. 135, 174709 (2011); doi: 10.1063/1.3657408

4.8.3.1. Example use of the analysis classes

4.8.3.1.1. HydrogenBondLifetimes

Analyzing hydrogen bond lifetimes (HBL) HydrogenBondLifetimes, both continuos and intermittent. In this case we are analyzing how residue 38 interact with a water sphere of radius 6.0 centered on the geometric center of protein and residue 42. If the hydrogen bond lifetimes are very stable, we can assume that residue 38 is hydrophilic, on the other hand, if the are very unstable, we can assume that residue 38 is hydrophobic:

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import HydrogenBondLifetimes as HBL

u = MDAnalysis.Universe(pdb, trajectory)
selection1 = "byres name OH2 and sphzone 6.0 protein and resid 42"
selection2 = "resid 38"
HBL_analysis = HBL(universe, selection1, selection2, 0, 2000, 30)
HBL_analysis.run()
time = 0
#now we print the data ready to plot. The first two columns are the HBLc vs t
#plot and the second two columns are the HBLi vs t graph
for HBLc, HBLi in HBL_analysis.timeseries:
    print("{time} {HBLc} {time} {HBLi}".format(time=time, HBLc=HBLc, HBLi=HBLi))
    time += 1

#we can also plot our data
plt.figure(1,figsize=(18, 6))

#HBL continuos
plt.subplot(121)
plt.xlabel('time')
plt.ylabel('HBLc')
plt.title('HBL Continuos')
plt.plot(range(0,time),[column[0] for column in HBL_analysis.timeseries])

#HBL intermitent
plt.subplot(122)
plt.xlabel('time')
plt.ylabel('HBLi')
plt.title('HBL Intermitent')
plt.plot(range(0,time),[column[1] for column in HBL_analysis.timeseries])

plt.show()

where HBLc is the value for the continuos hydrogen bond lifetimes and HBLi is the value for the intermittent hydrogen bond lifetime, t0 = 0, tf = 2000 and dtmax = 30. In this way we create 30 windows timestep (30 values in x axis). The continuos hydrogen bond lifetimes should decay faster than intermittent.

4.8.3.1.2. WaterOrientationalRelaxation

Analyzing water orientational relaxation (WOR) WaterOrientationalRelaxation. In this case we are analyzing “how fast” water molecules are rotating/changing direction. If WOR is very stable we can assume that water molecules are rotating/changing direction very slow, on the other hand, if WOR decay very fast, we can assume that water molecules are rotating/changing direction very fast:

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import WaterOrientationalRelaxation as WOR

u = MDAnalysis.Universe(pdb, trajectory)
selection = "byres name OH2 and sphzone 6.0 protein and resid 42"
WOR_analysis = WOR(universe, selection, 0, 1000, 20)
WOR_analysis.run()
time = 0
#now we print the data ready to plot. The first two columns are WOR_OH vs t plot,
#the second two columns are WOR_HH vs t graph and the third two columns are WOR_dip vs t graph
for WOR_OH, WOR_HH, WOR_dip in WOR_analysis.timeseries:
      print("{time} {WOR_OH} {time} {WOR_HH} {time} {WOR_dip}".format(time=time, WOR_OH=WOR_OH, WOR_HH=WOR_HH,WOR_dip=WOR_dip))
      time += 1

#now, if we want, we can plot our data
plt.figure(1,figsize=(18, 6))

#WOR OH
plt.subplot(131)
plt.xlabel('time')
plt.ylabel('WOR')
plt.title('WOR OH')
plt.plot(range(0,time),[column[0] for column in WOR_analysis.timeseries])

#WOR HH
plt.subplot(132)
plt.xlabel('time')
plt.ylabel('WOR')
plt.title('WOR HH')
plt.plot(range(0,time),[column[1] for column in WOR_analysis.timeseries])

#WOR dip
plt.subplot(133)
plt.xlabel('time')
plt.ylabel('WOR')
plt.title('WOR dip')
plt.plot(range(0,time),[column[2] for column in WOR_analysis.timeseries])

plt.show()

where t0 = 0, tf = 1000 and dtmax = 20. In this way we create 20 windows timesteps (20 values in the x axis), the first window is created with 1000 timestep average (1000/1), the second window is created with 500 timestep average(1000/2), the third window is created with 333 timestep average (1000/3) and so on.

4.8.3.1.3. AngularDistribution

Analyzing angular distribution (AD) AngularDistribution for OH vector, HH vector and dipole vector. It returns a line histogram with vector orientation preference. A straight line in the output plot means no preferential orientation in water molecules. In this case we are analyzing if water molecules have some orientational preference, in this way we can see if water molecules are under an electric field or if they are interacting with something (residue, protein, etc):

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import AngularDistribution as AD

u = MDAnalysis.Universe(pdb, trajectory)
selection = "byres name OH2 and sphzone 6.0 (protein and (resid 42 or resid 26) )"
bins = 30
AD_analysis = AD(universe,selection,bins)
AD_analysis.run()
#now we print data ready to graph. The first two columns are P(cos(theta)) vs cos(theta) for OH vector ,
#the seconds two columns are P(cos(theta)) vs cos(theta) for HH vector and thirds two columns
#are P(cos(theta)) vs cos(theta) for dipole vector
for bin in range(bins):
      print("{AD_analysisOH} {AD_analysisHH} {AD_analysisDip}".format(AD_analysis.graph0=AD_analysis.graph[0][bin], AD_analysis.graph1=AD_analysis.graph[1][bin],AD_analysis.graph2=AD_analysis.graph[2][bin]))

#and if we want to graph our results
plt.figure(1,figsize=(18, 6))

#AD OH
plt.subplot(131)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for OH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[0][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[0][:-1]])

#AD HH
plt.subplot(132)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for HH')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[1][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[1][:-1]])

#AD dip
plt.subplot(133)
plt.xlabel('cos theta')
plt.ylabel('P(cos theta)')
plt.title('PDF cos theta for dipole')
plt.plot([float(column.split()[0]) for column in AD_analysis.graph[2][:-1]],[float(column.split()[1]) for column in AD_analysis.graph[2][:-1]])

plt.show()

where P(cos(theta)) is the angular distribution or angular probabilities.

4.8.3.1.4. MeanSquareDisplacement

Analyzing mean square displacement (MSD) MeanSquareDisplacement for water molecules. In this case we are analyzing the average distance that water molecules travels inside protein in XYZ direction (cylindric zone of radius 11[nm], Zmax 4.0[nm] and Zmin -8.0[nm]). A strong rise mean a fast movement of water molecules, a weak rise mean slow movement of particles:

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import MeanSquareDisplacement as MSD

u = MDAnalysis.Universe(pdb, trajectory)
selection = "byres name OH2 and cyzone 11.0 4.0 -8.0 protein"
MSD_analysis = MSD(universe, selection, 0, 1000, 20)
MSD_analysis.run()
#now we print data ready to graph. The graph
#represents MSD vs t
time = 0
for msd in MSD_analysis.timeseries:
      print("{time} {msd}".format(time=time, msd=msd))
      time += 1

#Plot
plt.xlabel('time')
plt.ylabel('MSD')
plt.title('MSD')
plt.plot(range(0,time),MSD_analysis.timeseries)
plt.show()

4.8.3.1.5. SurvivalProbability

Analyzing survival probability (SP) SurvivalProbability for water molecules. In this case we are analyzing how long water molecules remain in a sphere of radius 12.3 centered in the geometrical center of resid 42, 26, 34 and 80. A slow decay of SP means a long permanence time of water molecules in the zone, on the other hand, a fast decay means a short permanence time:

import MDAnalysis
from MDAnalysis.analysis.waterdynamics import SurvivalProbability as SP
import matplotlib.pyplot as plt

universe = MDAnalysis.Universe(pdb, trajectory)
selection = "byres name OH2 and sphzone 12.3 (resid 42 or resid 26 or resid 34 or resid 80) "
sp = SP(universe, selection, verbose=True)
sp.run(start=0, stop=100, tau_max=20)
tau_timeseries = sp.tau_timeseries
sp_timeseries = sp.sp_timeseries

# print in console
for tau, sp in zip(tau_timeseries, sp_timeseries):
      print("{time} {sp}".format(time=tau, sp=sp))

# plot
plt.xlabel('Time')
plt.ylabel('SP')
plt.title('Survival Probability')
plt.plot(taus, sp_timeseries)
plt.show()

4.8.3.2. Output

4.8.3.2.1. HydrogenBondLifetimes

Hydrogen bond lifetimes (HBL) data is returned per window timestep, which is stored in HydrogenBondLifetimes.timeseries (in all the following descriptions, # indicates comments that are not part of the output):

results = [
    [ # time t0
        <HBL_c>, <HBL_i>
    ],
    [ # time t1
        <HBL_c>, <HBL_i>
    ],
    ...
 ]

4.8.3.2.2. WaterOrientationalRelaxation

Water orientational relaxation (WOR) data is returned per window timestep, which is stored in WaterOrientationalRelaxation.timeseries:

results = [
    [ # time t0
        <WOR_OH>, <WOR_HH>, <WOR_dip>
    ],
    [ # time t1
        <WOR_OH>, <WOR_HH>, <WOR_dip>
    ],
    ...
 ]

4.8.3.2.3. AngularDistribution

Angular distribution (AD) data is returned per vector, which is stored in AngularDistribution.graph. In fact, AngularDistribution returns a histogram:

results = [
    [ # OH vector values
      # the values are order in this way: <x_axis  y_axis>
        <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ...
    ],
    [ # HH vector values
        <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ...
    ],
    [ # dip vector values
       <cos_theta0 ang_distr0>, <cos_theta1 ang_distr1>, ...
    ],
 ]

4.8.3.2.4. MeanSquareDisplacement

Mean Square Displacement (MSD) data is returned in a list, which each element represents a MSD value in its respective window timestep. Data is stored in MeanSquareDisplacement.timeseries:

results = [
     #MSD values orders by window timestep
        <MSD_t0>, <MSD_t1>, ...
 ]

4.8.3.2.5. SurvivalProbability

Survival Probability (SP) computes two lists: a list of taus (SurvivalProbability.tau_timeseries) and a list of their corresponding mean survival probabilities (SurvivalProbability.sp_timeseries). Additionally, a list SurvivalProbability.sp_timeseries_data is provided which contains a list of SPs for each tau, which can be used to compute their distribution, etc.

results = [ tau1, tau2, …, tau_n ], [ sp_tau1, sp_tau2, …, sp_tau_n]

Additionally, for each

4.8.3.3. Classes

class MDAnalysis.analysis.waterdynamics.HydrogenBondLifetimes(universe, selection1, selection2, t0, tf, dtmax, nproc=1)[source]

Hydrogen bond lifetime analysis

This is a autocorrelation function that gives the “Hydrogen Bond Lifetimes” (HBL) proposed by D.C. Rapaport [Rapaport1983]. From this function we can obtain the continuous and intermittent behavior of hydrogen bonds in time. A fast decay in these parameters indicate a fast change in HBs connectivity. A slow decay indicate very stables hydrogen bonds, like in ice. The HBL is also know as “Hydrogen Bond Population Relaxation” (HBPR). In the continuos case we have:

\[C_{HB}^c(\tau) = \frac{\sum_{ij}h_{ij}(t_0)h'_{ij}(t_0+\tau)}{\sum_{ij}h_{ij}(t_0)}\]

where \(h'_{ij}(t_0+\tau)=1\) if there is a H-bond between a pair \(ij\) during time interval \(t_0+\tau\) (continuos) and \(h'_{ij}(t_0+\tau)=0\) otherwise. In the intermittent case we have:

\[C_{HB}^i(\tau) = \frac{\sum_{ij}h_{ij}(t_0)h_{ij}(t_0+\tau)}{\sum_{ij}h_{ij}(t_0)}\]

where \(h_{ij}(t_0+\tau)=1\) if there is a H-bond between a pair \(ij\) at time \(t_0+\tau\) (intermittent) and \(h_{ij}(t_0+\tau)=0\) otherwise.

Parameters:
  • universe (Universe) – Universe object
  • selection1 (str) – Selection string for first selection [‘byres name OH2’]. It could be any selection available in MDAnalysis, not just water.
  • selection2 (str) – Selection string to analize its HBL against selection1
  • t0 (int) – frame where analysis begins
  • tf (int) – frame where analysis ends
  • dtmax (int) – Maximum dt size, dtmax < tf or it will crash.
  • nproc (int) – Number of processors to use, by default is 1.

New in version 0.11.0.

run(**kwargs)[source]

Analyze trajectory and produce timeseries

class MDAnalysis.analysis.waterdynamics.WaterOrientationalRelaxation(universe, selection, t0, tf, dtmax, nproc=1)[source]

Water orientation relaxation analysis

Function to evaluate the Water Orientational Relaxation proposed by Yu-ling Yeh and Chung-Yuan Mou [Yeh1999]. WaterOrientationalRelaxation indicates “how fast” water molecules are rotating or changing direction. This is a time correlation function given by:

\[C_{\hat u}(\tau)=\langle \mathit{P}_2[\mathbf{\hat{u}}(t_0)\cdot\mathbf{\hat{u}}(t_0+\tau)]\rangle\]

where \(P_2=(3x^2-1)/2\) is the second-order Legendre polynomial and \(\hat{u}\) is a unit vector along HH, OH or dipole vector.

Parameters:
  • universe (Universe) – Universe object
  • selection (str) – Selection string for water [‘byres name OH2’].
  • t0 (int) – frame where analysis begins
  • tf (int) – frame where analysis ends
  • dtmax (int) – Maximum dt size, dtmax < tf or it will crash.

New in version 0.11.0.

static lg2(x)[source]

Second Legendre polynomial

run(**kwargs)[source]

Analyze trajectory and produce timeseries

class MDAnalysis.analysis.waterdynamics.AngularDistribution(universe, selection_str, bins=40, nproc=1, axis='z')[source]

Angular distribution function analysis

The angular distribution function (AD) is defined as the distribution probability of the cosine of the \(\theta\) angle formed by the OH vector, HH vector or dipolar vector of water molecules and a vector \(\hat n\) parallel to chosen axis (z is the default value). The cosine is define as \(\cos \theta = \hat u \cdot \hat n\), where \(\hat u\) is OH, HH or dipole vector. It creates a histogram and returns a list of lists, see Output. The AD is also know as Angular Probability (AP).

Parameters:
  • universe (Universe) – Universe object
  • selection (str) – Selection string to evaluate its angular distribution [‘byres name OH2’]
  • bins (int (optional)) – Number of bins to create the histogram by means of numpy.histogram()
  • axis ({'x', 'y', 'z'} (optional)) – Axis to create angle with the vector (HH, OH or dipole) and calculate cosine theta [‘z’].

New in version 0.11.0.

run(**kwargs)[source]

Function to evaluate the angular distribution of cos(theta)

class MDAnalysis.analysis.waterdynamics.MeanSquareDisplacement(universe, selection, t0, tf, dtmax, nproc=1)[source]

Mean square displacement analysis

Function to evaluate the Mean Square Displacement (MSD). The MSD gives the average distance that particles travels. The MSD is given by:

\[\langle\Delta r(t)^2\rangle = 2nDt\]

where \(r(t)\) is the position of particle in time \(t\), \(\Delta r(t)\) is the displacement after time lag \(t\), \(n\) is the dimensionality, in this case \(n=3\), \(D\) is the diffusion coefficient and \(t\) is the time.

Parameters:
  • universe (Universe) – Universe object
  • selection (str) – Selection string for water [‘byres name OH2’].
  • t0 (int) – frame where analysis begins
  • tf (int) – frame where analysis ends
  • dtmax (int) – Maximum dt size, dtmax < tf or it will crash.

New in version 0.11.0.

run(**kwargs)[source]

Analyze trajectory and produce timeseries

class MDAnalysis.analysis.waterdynamics.SurvivalProbability(universe, selection, t0=None, tf=None, dtmax=None, verbose=False)[source]

Survival probability analysis

Function to evaluate the Survival Probability (SP). The SP gives the probability for a group of particles to remain in certain region. The SP is given by:

\[P(\tau) = \frac1T \sum_{t=1}^T \frac{N(t,t+\tau)}{N(t)}\]

where \(T\) is the maximum time of simulation, \(\tau\) is the timestep, \(N(t)\) the number of particles at time t, and \(N(t, t+\tau)\) is the number of particles at every frame from t to tau.

Parameters:
  • universe (Universe) – Universe object
  • selection (str) – Selection string; any selection is allowed. With this selection you define the region/zone where to analyze, e.g.: “resname SOL and around 5 (resname LIPID)” and “resname ION and around 10 (resid 20)” (see SP-examples )
  • verbose (Boolean) – If True, prints progress and comments to the console.

New in version 0.11.0.

run(tau_max=20, start=0, stop=None, step=1, verbose=False)[source]

Computes and returns the survival probability timeseries

Parameters:
  • start (int) – Zero-based index of the first frame to be analysed
  • stop (int) – Zero-based index of the last frame to be analysed (inclusive)
  • step (int) – Jump every step’th frame
  • tau_max (int) – Survival probability is calculated for the range \(1 <= au <= tau_max\)
  • verbose (Boolean) – Overwrite the constructor’s verbosity
Returns:

  • tau_timeseries (list) – tau from 1 to tau_max. Saved in the field tau_timeseries.
  • sp_timeseries (list) – survival probability for each value of tau. Saved in the field sp_timeseries.