Klauda Lab Simulation Wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
[[Analysis]]

======Brief instruction to obtain lipid wobbble correlation times======

Authors: Xiaohong Zhuang, Viviana Monje and Dr. Jeffery B. Klauda
Updated by Eric Zhang to handle any lipid type automatically

This brief instruction is used to calculate the time scale of lipid wobble, which describes the lipid mobility, and by calculating the correlation time (tau) of the cross-tail vectors of glycerol phospholipids (C22-C32) or sphingolipids (C4S-C2F). This could be adjusted to any vector of interest by changing the atom selection

Download Legacy Scripts: {{ ::wob-s.tar.gz |}}\\
Download Legacy Example: {{ :wob-e.tar.gz |}}\\
Download Modern Python Scripts: {{ :asym_lip_wobble.tar.gz |}}\\
Download Modern Python Example: {{ :asym_lip_wobble_eg.tar.gz |}}\\

If on ZT1, use the files in the following directory: /afs/shell.umd.edu/project/energybio/shared/jbklauda/scripts/lip_wobble/

Below, is the instructions to use the legacy scripts. Scroll down further to get instructions on the modern python workflow.
=====Calculate the correlation function=====

1.	Update **correl.csh** to your system

<code commend>
# Update the path of the input psf, crd file
set p = /lustre/xzhuang/simulation/lipids/dmpi/charmm-gui

# Update Name of the atoms, use C22 and C32 for phospholipids, and C4S and C2F for sphingolipids 
set atm1 = C22
set atm2 = C32

# Update Nf, number of frames per dcd file
set Nf = 1000
</code>

2.	Update **set-range.str**

<code commend>
! Update the following numbers
set rn = DMPI   ! Update the resname of lipid
set type = G   ! Update lipid type, put G for glycerol lipids, S for sphingolipids
set ti = 21      ! first trajectory
set tf = 50      ! last trajectory

! Manually update the numbers if the top and bottom in your system are not symmetric
calc Nh = @Nlip / 2   ! Calculate number of lipids on the top/bot leaflet
calc Nb = @Nh + 1    ! Calculate number of beginning bottom leaflet subset
</code>

//Note://

{{::dmpi-sa.jpg?300|}}{{::dmpi-corr.jpg?300|}}   

The example simulation has 50 dcd, with each dcd 1,000,000 time steps (i.e. 2ns), therefore total simulation time is 100 ns. Unlike other analyses, we need to use the longest possible equilibrated time range to calculate correlation function (usually 40% or more).
Here, the last 60% (last 60 ns, i.e. 60ns~100 ns) is used to calculate the correlation function (as SA/lipid indicates that the system reach equilibrium at by 60 ns).

As you can see in the graph of the example correlation function output in corr directory, the correlation is very weak in the tail, so we only use the last 65%~90% to do the exponential fitting to obtain the correlation time tau ( and also to calculate T1).
As you can see the following in **expfit/1_line_sed.scr**
<code commend>
linenew=$(( $line * 3 / 4 ))
</code>


3.	Submit the analysis job by command: ''correl.scr''

After run finished, the **cross-top.dat** and **cross-bot.dat** will be obtained in the corr directory. The first column is the time in ps, and the second 2nd is the correlation function C2(t).

4.	Go to corr directory, then calucate the top and bottom leaflet average by command ''1_comb_avg.scr''

5.	Adjust the time range (axis) in 2_correl.gnu, then plot the correlation function by command ''1_comb_avg.scr'' (if you use WinSCP), otherwise, you may plot it in Excel or any other software you like.

=====Fitting correlation function to exponential functions=====

1.	Go to expfit directory, obtain the last 80% of correlation by command ''1_line_sed.scr''

2.	Run the exponential fit by command ''2_expfit.csh''

The current scripts fitting both 2nd order exponential function since it fitting the 2nd order exponential is more likely to work than the 3rd order fitting based on my experience. The results from the two types of fitting are generally close.  In both fittings, the longest correlation time corresponds to lipid wobble movement time scale. 

In **cross-fit2.dat**, CF values in the end are, a1 plateau a2 t1 t2, the longest time scale t2 (i.e. CF(4) ) is what we need.
In **cross-fit3.dat**, CF values in the end are, a1 a2 plateau a3 t1 t2 t3, the longest time scale t3 (i.e. CF(6) ) is what we need.

If the fitting work fine, you will see the statement “STOP Normal Exit statement executed” on the screen. 

If the 2nd order fitting did not go well, the following two things that you can try. But please discuss with Dr.Klauda about these.

a.	Adjust the portion of correlation function to fit by modify the **1_line_sed.scr**. 

The current script uses first 75% data in the following line, you may change it to 80% (must use integer fraction format instead of floating point (decimal) here, i.e. 4/5), or 90% (9/10), or any other fraction depends your correlation function graph.

<code commend>
linenew=$(( $line * 3 / 4 ))
</code>

b.	Use longer equilibrated time trajectory range by adjust ti value in **set-range.str**, and redo the correlation function calculation.

3.	This step is just to check the wellness of the fitting. 

You may adjust the time range (axis) in 3_expfit_cross.gnu, then plot the correlation function and its exponential fitting by command ''3_plot_expfit.scr'' (if you use WinSCP), otherwise, you may plot it in Excel or any other software you like.


=====bilayer wobble=====

Please make the following edits to adjust the code for individual lipids to all glycerol lipids. (Sphino lipids would need to be calculated separately, by adjusting the cross-chain vector atom selection in **correl.csh**)

1. **set-range.str**

a. delete or comment out the following line. No need to specify lipid if for overall bilayer wobble
<code commend>
! set rn = xxx
</code>
b. change lipid count/selection. With //type C32//, we are able to select all glycerol lipids.
<code commend>
if @type .eq. G then
        define lipid select type C32 end     ! all glycerol lipids
</code>
2. **cross-*.inp**

a. change vector atom selection line.

change this line 
<code commend>
enter v@K vect xyz MEMB @R @atm1 MEMB @R @atm2
</code>
to 
<code commend>
enter v@K vect xyz sele resid @R .and. (type @atm1 .or. type @atm2) end
</code>

3. everything else the same as above sections

======Modern Python Workflow======
This workflow mirrors the old setup, but includes scripts that dynamically set up CHARMM files based on your membrane setup. I've also included many advanced options when fitting exponentials, most of which may not be needed. A detailed outline of the scripts is also included in a dev notes file, in case you want to modify them. Feel free to message ezhang32@umd.edu for any questions/clarifications about the scripts. Examples for all steps can be found in asym_lip_wobble_eg at the start.

=====Step 1: Define System Parameters=====
Open ''config.py''. This file contains everything you need to edit based on your system and/or preferences. First, check the contents of cross_tail_dict. This dictionary/map specifies the atom pair that defines your cross-tail vector. By default, it includes selections for glycerophospholipids (GL), sphingolipids (SL), and lipid A (LIPA). You can customize it by adding your own (lipid type, atom pair list) entry, as I did with NAAA. The atoms should sit on separate acyl chains, be at roughly the same z-level, and be near the base of the chain. 

Next, modify/populate the dictionary called config_dict, which contains subdirectories defining the lipid residues, top leaflet counts, and bottom leaflet counts for each residue. If a lipid only exists in one leaflet, just set that count to 0 in the appropriate subdirectory. 

Next, edit the workspace fields by specifying your charmm/dyn directories and your equilibrated DCD range (as wide as possible). 

Finally (Optional), specify what kind of analysis you want by adjusting the fields. The default implementation used in prior publications is a 2nd-order exponential fit, an unweighted fit, with 75% of the data untruncated. There are a lot of other options, which are explained inside the file via comments.

=====Step 2: Run Scripts=====
Run ''0_setup.py''. This generates your CHARMM scripting files. Double-check that there are no errors.

Next, run ''1_correl.scr''. This submits your Slurm jobs for the CHARMM scripts. This is usually pretty quick and can take anywhere from minutes to an hour.

Next, run ''2_analysis.py.'' This creates a results folder which holds all of the generated plots and wobble_results.txt, which contains your wobble parameters. I suggest playing around with the analysis fields in config.py to see how your results change. 

If these fits fail or give unphysical results, ask Professor Klauda for advice. But in general, boundary-pinning errors indicate you don't have enough data (the graph looks like it decays to zero, but a proper correlation function should plateau). You can try adjusting the truncated window to capture more or less of the data. For more rigid systems, 3rd-order exponential fitting may be preferred, and/or you can raise the upper bounds manually by editing the "tau1_ub" variable inside src/utils/analysis_utils.py (the lower bound is 0, which should not be touched). Using weighted least squares is a fancy alternative that enables you to use most of the data with minimal, if any, truncation in scenarios where data is limited. Unweighted least squares seems to be what is reported in the literature in most cases.