Scripts.DatabankLib.analyze_nmrpca module

This code aims to calculate the relaxation time based on the PCAlipids analysis.

For details check:

Principal Component Analysis of Lipid Molecule Conformational Changes in Molecular Dynamics Simulations Pavel Buslaev, Valentin Gordeliy, Sergei Grudinin, and Ivan Gushchin. Journal of Chemical Theory and Computation (2016), 12(3), 1019–1028. DOI: 10.1021/acs.jctc.5b01106

and

Principal component analysis highlights the influence of temperature, curvature and cholesterol on conformational dynamics of lipids Pavel Buslaev, Khalid Mustafin, and Ivan Gushchin Biochimica et Biophysica Acta (BBA) - Biomembranes Volume 1862, Issue 7, 1 July 2020, 183253 DOI: 10.1016/j.bbamem.2020.183253

The main idea is to run PCA on concatenated lipid trajectory, calculate the autocorrelation times, which are linearly related to equilibration times of individual lipids. The parameters of linear transform are calculated based on the trajectories from NMRlipids databank.

The initial code was developed by Dr. Pavel Buslaev (pavel.i.buslaev@jyu.fi), Dr. Samuli Olilla, Patrik Kula, Alexander Kuzmin

MEM_CUTOFF

Environmental variable which sets the trajectory size cutoff (in Bytes) after which the algorithm will start skipping frames for in-memory representation for computing ACF.

Default: 1,500,000,000.

class Scripts.DatabankLib.analyze_nmrpca.Parser(system: System, eq_time_fname='eq_times.json', path=None, v=True)

Bases: object

Class Parser is a basic class to work with the trajectory. It has basic utility for preparing trajectories:

1. Checking if simulation has a correct name

2. Checking if trajectory is downloaded. Downloading, if not.

3. Calling AmberMerger to merge head groups and tails for Amber trajectories

4. Concatenate trajectories

Parameters:

• system – simulation data

• eq_time_fname – name of the output file

• path – name of a particular trajectory

• v – verbosity for testing

validate_path()

Path validation. Behaviour depends on the input.

1. If ther were any errors, validation fails. Currently the only error tested is that .xtc and .tpr files are not present. This is the case for non-GROMACS trajectories.

2. If path is not provided, parser is iterating over all trajectories.

3. If path is provided and current path is equal to the provided one, parser reports that it finds the trajectory for the analysis.

TODO: must be removed. Substitute path-validation part.

download_traj()

TODO: must be removed. We don’t download trajectory in the analysis part

Basic trajectory and TPR download.

prepare_gmx_traj()

Preparing trajectory. If centered trajectory is found, use it. If whole trajectory is found, use it. Otherwise, call gmx trjconv to make whole trajectory. The selected trajectory is loaded into Universe.

prepare_OpenMM_traj()

prepare_traj()

concatenate_traj()

Create Concatenator and corresponding concatenated trajectories for all lipids available for the trajectory.

dump_data(data)

Write data to json file

class Scripts.DatabankLib.analyze_nmrpca.Topology(traj, lipid_resname: str, mapping_dict: dict)

Bases: object

Class Topology is a class needed to extract lipid specific data and also to merge lipids from Amber trajectories, where lipid is often represented as 3 residue types. It has the following methods:

1. Simple constructor, which sets the force field,residue names, trajectory, and loads mapping_file

2. Mapping file loader

3. Function that outputs atomNames

4. Checker if the merge is needed. Currently defunct

5. Runner, that returns lists for head, tail1, tail2 orders for merging

Currently defunct.

Parameters:

• traj – MDAnalysis trajectory

• lipid_resname – resname of the lipid (str)

• mapping_dict – preloaded mapping_dict

atom_names()

Extract all names of heavy atoms from the mapping

is_merge_needed()

Checker for merge. Currently it checks if the RESIDUE key is in the mapping file.

NOTE: This has to be changed if the content of mapping file changes

get_lipid_resnames()

Helper function that gets the residue names for lipid, if merge is needed

assign_resnames(resnames)

Helper function that finds head, sn-1 and sn-2 tails

NOTE: currently there are only lipids with 2 different tails available in databank: e.g. POPC or POPG. This leads to different names of tails. This won’t be the case for DOPC. Currently the algorithm is using that all three groups differ

run_merger()

Find lipid tails that correspond to a particular head-group.

NOTE: currently there are only lipids with 2 different tails available in databank: e.g. POPC or POPG. This leads to different names of tails. This won’t be the case for DOPC. Currently the algorithm is using that all three groups differ

TODO: get the correspondence from structure

class Scripts.DatabankLib.analyze_nmrpca.Concatenator(topology: Topology, traj, lipid_resname: str)

Bases: object

Class Concatenator is a class needed to concatenate trajectory for lipid types. It has the following methods:

1. Simple constructor, which sets the topology,residue names, and trajectory

2. concatenateTraj() method to do the basic by lipid concatenation

3. alignTraj() method to perform the alignment of the concatenated trajectory

4. The enveloping concatenate method

Parameters:

• topology – topology for lipid

• traj – MDAnalysis trajectory

• lipid_resname – lipid resname in the trajectory

concatenate_traj()

concatenateTraj performs basic trajectory concatination. First, it extracts coordinates from trajectory, next, it reshapes the coordinate array, swaps time and resid axes to obtain continuous trajectories of individual lipids (this is needed for autocorrelation time analysis), and finally merges individual lipid trajectories.

concatenate_traj_with_merging()

concatenateTrajWithMerging performs basic trajectory concatination. In contrast to basic concatenateTraj it additionally merges splitted lipids. First, it creates extracts coordinates from trajectory, next, it reshapes the coordinate array, swaps time and resid axes to obtain continuous trajectories of individual lipids (this is needed for autocorrelation time analysis), and finally merges individual lipid trajectories.

align_traj(concatenated_traj)

alignTraj alignes the concatenated trajectory in two steps: (1) it computes average structure after alignment to the first frame, and (2) it alignes the structure to the calculated average structure in (1).

concatenate()

Simple enveloping function to perform concatenation

class Scripts.DatabankLib.analyze_nmrpca.PCA(aligned_traj, av_pos, n_lipid, n_frames, traj_time)

Bases: object

Class PCA is a class that actually performs PCA. It has the following methods:

1. Simple constructor, which sets the aligned trajtory, average coordinates, number of lipids, number of frames in the concatenated trajectory and trajectory length in ns

2. PCA() prepares the trajectory coordinates for analysis and calculates principal components

3. get_proj() projects the trajectory on principal components

4. get_lipid_autocorrelation() calculates the autocorrelation timeseries for individual lipid

5. get_autocorrelations() calculates the autocorrelation timeseries for trajectory

Parameters:

• aligned_traj – np.array with positions of concatenated and aligned trajectory

• av_pos – np.array with average positions for the lipid

• n_lipid – number of lipids of particular type in the system

• n_frames – number of frames in the concatenated trajectory

• traj_time – trajectory length in ns

PCA()

PCA calculates the PCA. First the data is centered and then covariance matrix is calculated manually.

get_proj(cdata)

Projecting the trajectory on the 1st principal component

get_lipid_autocorrelation(data, variance, mean)

Autocorrelation calculation for individual lipid.

get_autocorrelations()

Autocorrelation calculation for the trajectory.

class Scripts.DatabankLib.analyze_nmrpca.TimeEstimator(autocorrelation)

Bases: object

Class TimeEstimator is a class that estimates equilbration time from autocorrelation data. It includes the following methods:

1. Simple constructor, which sets the autocorrelation data

2. Helper method get_nearest_value() that finds the interval in the data where the autocorrelation falls below specific value

3. timerelax() method calculates the logarithms of autocorrelation and calculates the decay time

4. calculate_time() method calculates the estimated equilibration time

Parameters:

autocorrelation – autocorrelation data

get_nearest_value(iterable, value)

get_nearest_value method return the indices that frame the particular value As an input it gets an array, which is expected to decay. The method tries to find and index (index_A), for which the array gets below the cutoff value for the first time. Then, for index_A-1 the array was larger than the cutoff value for the last time. It means that the function, represented by the array data is equal to cutoff somewhere between timesteps that correspond to index_A-1 and index_A. It can happen, that this index_A is equal 0. Then, method searches for the first occurance, where array is smaller than array[0]. The method returns the interval inbetween the array gets below cutoff.

timerelax()

Estimates autocorrelation decay time.

calculate_time()

Basic enveloping method that estimates equilibration time from autocorrelation decay time. They are linearly connected and the coefficient is calculated experimentally.