Python script to compute rotational correlation function of a protein backbone N-H vector
NMR and fluorescence spectroscopy are examples of experimental methods which are sensitive to inter- and intramolecular dynamics across range of timescales, from picoseconds to hours. However, only all-atom MD simulations can provide atomistic insight into these motions. Scripts within this repository serve for calculation of rotational correlation function, C(τc), of a vector of interest, using trajectories obtained from all-atom MD simulations. Reorientation of a vector is determined by both overall rotational diffusion of a protein and internal motions. Since global rotational diffusion is much slower compared to internal motions of backbone or side-chain vectors, rotational correlation function C(τc) can be separated into correlation function for global and internal rotational motions1:
C(τc) = Cinternal(τc)Cglobal(τc)
If MD trajectory is aligned, overall rotational diffusion of a protein will be removed, and only timescales of internal rotational motions can be determined from the analysis of correlation function.
Rotational correlation function can be expressed in terms of second Legendre polynomial of the dot product of the unit vector e at time t and t+τc The scalar product between two unit vectors equals to the cosine of the angle θ between them, which further simplifies the equation:
C(τc) = ⟨P2(e(t)⋅e(t+τc))⟩ = ⟨½(3cos2θ(t+τc)-1)⟩
Within this repository, two scripts are provided, namely RotCorrFun_full_traj.py and RotCorrFun_subtraj.py. In the first case, rotational correlation function is calculated for the entire trajectory. In the latter case, trajectory is divided in N subtrajectories of subtraj_length. The idea is that by subsequently averaging over several subtrajectories, less noisy C(τc) can be obtained. Length of subtrajectory should be fewfold longer than the overal rotational diffusion time of the protein (if its estimate is available from other source).
Particularly convenient is to calculate rotational correlation function for the backbone N-H bond vector, because this allows to directly compare the results with NMR data of {1H}-15N labelled protein:
vector_origin = "resname LYS and name N"
vector_end = "resname LYS and name H"
In this particular example, I studied reorientation of N-H bond vectors in lysine residues. Script calculates and saves as .txt file individual lysine correlation functions for each of the (sub)trajectories, as well as the average correlation function across all lysine residues within a (sub)trajectory. Depending on the system under the study, if local environments for different residues vary, such averaging will likely not be meaningful. Therefore, it is left to the user to decide if C(τc) should be averaged over residues and/or trajectories.
Using mdtraj text-based atom selection language, reorientation of any other vector of choice can be monitored by defining its origin and end, as given in the example above.
C(τc) can be modelled as weighted sum of exponential decays, whose characteristic times correspond to rotational correlation times of internal and global rotational motions. Furthermore, amplitudes of exponential terms relate to so-called order parameter, S2, which is a measure of bond flexibility, and in the case of backbone vectors, it relates to secondary structure. Full framework related to analysis of rotational correlation function using fluorescence, NMR and MD data is presented in work by Möckel et al, 2019.1
- time-ordered trajectory obtained from all-atom MD simulations, in any of the mdtraj compatible formats (dcd, xtc, nc,..). If global rotational diffusion of the protein is to be quantified, trajectory should not be aligned. Otherwise, if only internal motions of backbone and side-chain vectors are of interest, trajectory should be superimposed on the first frame in trajectory.
- topology as .pdb file
RotCorrFun_full_traj.py and RotCorrFun_subtraj.py are Python scripts developed on Python 3.8.8. Scripts were tested under the following configuration:
- Windows 10
- Python 3.8.8
- numpy 1.23.0
- mdtraj 1.9.4
- Mockel, C.; Kubiak, J.; Schillinger, O.; Kuhnemuth, R.; Della Corte, D.; Schroder, G. F.; Willbold, D.; Strodel, B.; Seidel, C. A. M.; Neudecker, P., Integrated NMR, Fluorescence, and Molecular Dynamics Benchmark Study of Protein Mechanics and Hydrodynamics. J Phys Chem B 2019, 123 (7), 1453-1480.
- Milana Popara