How to calculate Global Elasticity?

To calculate global elasticity, dnaMD globalElasticity command can be used. However, it takes HDF5 file as input. Following steps can be performed to generate HDF5 files. The tutorial file can be downloaded here. We will prepare HDF5 file for both free and bound DNA.

Calculate stretching twisting and bending motions

Both free and bound DNA is superimposed on to the same DNA structure. Careful that bending calculation is fitting dependent. Therefore, at first we aligned both free and bound DNA to a common DNA structure as follows

In [ ]:
%%bash

# Align bound DNA
echo 6 0 | gmx trjconv -f inputs/F1_complex_DNA.xtc -s inputs/dna.tpr -n inputs/dna.ndx -o complex_dna_aligned.xtc -fit rot+trans

# Align free DNA
echo 6 0 | gmx trjconv -f inputs/F1_free_DNA.xtc -s inputs/dna.tpr -n inputs/dna.ndx -o free_dna_aligned.xtc -fit rot+trans

1. Run do_x3dna on DNA trajectory, -nofit is used because DNA is already superimposed to a common DNA structure using trjconv.

In [ ]:
%%bash

# For free DNA
echo 2 | do_x3dna -f free_dna_aligned.xtc -s inputs/dna.tpr -n inputs/dna.ndx -ref -noavg -nofit -name free
mv *_free.dat outputs/.

# For bound DNA
echo 2 | do_x3dna -f complex_dna_aligned.xtc -s inputs/dna.tpr -n inputs/dna.ndx -ref -noavg -nofit -name bound
mv *_bound.dat outputs/.

2. Run dnaMD to extract the parameters from do_x3dna output files and save as HDF5 file. Also calculate global axis, curvature and tangents.

In [8]:
%%bash

# For free DNA
dnaMD saveAsH5 -tbp 27 -i outputs/L-BPS_free.dat,outputs/L-BPH_free.dat,outputs/HelAxis_free.dat -o free_dna.h5
dnaMD axisCurv -tbp 27 -bs 2 -be 25 -ctan -scp 100 -s 1000 -cta 30 -io free_dna.h5 -ap free_dna_axis.pdb

# For bound DNA
dnaMD saveAsH5 -tbp 27 -i outputs/L-BPS_bound.dat,outputs/L-BPH_bound.dat,outputs/HelAxis_bound.dat -o bound_dna.h5
dnaMD axisCurv -tbp 27 -bs 2 -be 25 -ctan -scp 100 -s 1000 -cta 30 -io bound_dna.h5 -ap bound_dna_axis.pdb

Loading parameters: ['shift', 'slide', 'rise', 'tilt', 'roll', 'twist']
Reading file : L-BPS_free.dat
Reading frame 5000
Finished reading.... Total number of frame read =  5001

Loading parameters: ['x-disp', 'y-disp', 'h-rise', 'inclination', 'tip', 'h-twist']
Reading file : L-BPH_free.dat
Reading frame 5000
Finished reading.... Total number of frame read =  5001

Loading parameters: ['helical x-axis', 'helical y-axis', 'helical z-axis']
Reading file : HelAxis_free.dat
Reading frame 5000
Finished reading.... Total number of frame read =  5001
Fitting spline curve on helical axis of frame 5000 out of 5001 frames
Finished spline curve fitting...
 ... Calculating curvature and tangents
                              ... Finished

Loading parameters: ['shift', 'slide', 'rise', 'tilt', 'roll', 'twist']
Reading file : L-BPS_bound.dat
Reading frame 10000
Finished reading.... Total number of frame read =  12001

Loading parameters: ['x-disp', 'y-disp', 'h-rise', 'inclination', 'tip', 'h-twist']
Reading file : L-BPH_bound.dat
Reading frame 10000
Finished reading.... Total number of frame read =  12001

Loading parameters: ['helical x-axis', 'helical y-axis', 'helical z-axis']
Reading file : HelAxis_bound.dat
Reading frame 10000
Finished reading.... Total number of frame read =  12001
Fitting spline curve on helical axis of frame 12000 out of 12001 frames
Finished spline curve fitting...
 ... Calculating curvature and tangents
                              ... Finished

Now, we have HDF5 files of both free and bounds DNA. It can be used for the calculation of elastic properties. These files can be used with either dnaMD Python module or dnaMD globalElasticity.

Bending Stretching Twist Modulus

Following command calculate Bending Stretching Twist modulus matrix. Output matrix will be stored in csv file. Elastic modulus matrix is printed as output and average values of contour length and cumulative twist angle is also printed.

In [1]:
%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -paxis X -o elastic_modulus_BST.csv
=========== Elastic Modulus Matrix ===========

   337.233      28.506      -6.984     -20.710
    28.506     373.863       7.977      39.509
    -6.984       7.977    1080.722     -96.198
   -20.710      39.509     -96.198     448.349

=========== ====================== ===========

=========================  Average Values  ==========================

Bending-1 Angle    Bending-2 Angle    Contour Length    Sum. Twist
          0.021            0.068            5.358            9.543

=========== ====================== ====================== ===========

The above modulus matrix is in this form:

\[\begin{split}\text{modulus matrix} = \begin{bmatrix} M_{Bx} & M_{Bx,By} & M_{Bx,S} & M_{Bx,T} \\ M_{Bx,By} & M_{By} & M_{By,S} & M_{By,T} \\ M_{Bx,S} & M_{By,S} & M_{S} & M_{S,T} \\ M_{Bx,T} & M_{Bx,T} & M_{S,T} & M_{T} \end{bmatrix}\end{split}\]

Where:

  • \(M_{Bx}\) - Bending-1 stiffness in one plane
  • \(M_{By}\) - Bending-2 stiffness in another orthogonal plane
  • \(M_{S}\) - Stretch Modulus
  • \(M_{T}\) - Twist rigidity
  • \(M_{Bx,By}\) - Bending-1 and Bending-2 coupling
  • \(M_{By,S}\) - Bending-2 and stretching coupling
  • \(M_{S,T}\) - Stretching Twsiting coupling
  • \(M_{Bx,S}\) - Bending-1 Stretching coupling
  • \(M_{By,T}\) - Bending-2 Twisting coupling
  • \(M_{Bx,T}\) - Bending-1 Twisting coupling

Stretching Twist Modulus

Following command calculate Bending Stretching Twist modulus matrix. Output matrix will be stored in csv file. Elastic modulus matrix is printed as output and average values of contour length and cumulative twist angle is also printed.

In [2]:
%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype ST -o elastic_modulus_ST.csv
=========== Elastic Modulus Matrix ===========

1080.380  -97.578
-97.578  442.494

=========== ====================== ===========

=========================  Average Values  ==========================

Contour Length   Sum. Twist
       5.358         9.543

=========== ====================== ====================== ===========

The above modulus matrix is in this form:

\[\begin{split}\text{modulus matrix} = \begin{bmatrix} M_{S} & M_{S,T} \\ M_{S,T} & M_{T} \end{bmatrix}\end{split}\]

Where:

  • \(M_{S}\) - Stretch Modulus
  • \(M_{T}\) - Twist rigidity
  • \(M_{S,T}\) - Stretching Twsiting coupling

Convergence in Modulus

Same command can be used to calculate elasticity as a function of time with option -ot/--output-time and save it in csv format file. This result can beused to check their convergence.

  • If this option is used, -fgap/--frame-gap is an essential option.
  • This options also gives final value and error as the output on display.
  • The output file is in csv format and can be opened as spreadsheet.

NOTE: Elastic properties cannot be calculated using a single frame because fluctuations are required. Therefore, here time means trajectory between zero time to given time.

In [3]:
%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -paxis X -fgap 100 -em block -ot modulus_time_BST.csv

# Print first  and last 3 line of output file
echo "====================================="
echo "Elastic modulus as a function of time"
echo "====================================="
head -4 modulus_time_BST.csv
printf ".\n.\n.\n"
tail -3 modulus_time_BST.csv
==============================================
Elasticity               Value         Error
----------------------------------------------
bend-1                 337.233         3.932
bend-2                 373.863         4.935
stretch               1080.722        24.015
twist                  448.349        29.356
bend-1-bend-2           28.506         4.685
bend-2-stretch           7.977         7.251
stretch-twist          -96.198        13.450
bend-1-stretch          -6.984         7.131
bend-2-twist            39.509        35.576
bend-1-twist           -20.710         5.081
==============================================

=====================================
Elastic modulus as a function of time
=====================================
#Time, bend-1, bend-2, stretch, twist, bend-1-bend-2, bend-2-stretch, stretch-twist, bend-1-stretch, bend-2-twist, bend-1-twist
1000.000, 584.19567, 363.51886, 1121.88830, 1506.57584, 41.49043, -101.57973, 112.06246, -254.22188, -74.92882, 252.68516
2000.000, 354.70978, 367.93829, 1032.49032, 662.58414, -18.43006, -53.55808, 9.82719, -81.25370, -207.81222, 111.70108
3000.000, 380.26015, 359.63694, 967.63012, 543.68667, 5.70038, -76.61296, -8.92719, -88.41420, -95.81366, 56.24528
.
.
.
118000.000, 335.33614, 371.13759, 1082.41850, 449.84105, 27.64596, 10.92246, -97.17016, -7.92023, 36.92839, -17.34603
119000.000, 335.79993, 373.00012, 1081.59948, 449.85747, 28.30290, 9.58966, -96.91853, -6.18165, 37.45497, -19.36140
120000.000, 337.23300, 373.86276, 1080.72201, 448.34909, 28.50566, 7.97658, -96.19760, -6.98413, 39.50860, -20.71008

Some of the plots from above data can be found here


Same as above but only for stretching and twisting motions.

In [4]:
%%bash

dnaMD globalElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -estype ST -paxis X -fgap 100 -gt "gmx analyze" -em "block" -ot modulus_time_ST.csv

# Print first  and last 3 line of output file
echo "====================================="
echo "Elastic modulus as a function of time"
echo "====================================="
head -4 modulus_time_ST.csv
printf ".\n.\n.\n"
tail -3 modulus_time_ST.csv
==============================================
Elasticity               Value         Error
----------------------------------------------
stretch               1080.380        22.847
twist                  442.494        19.169
stretch-twist          -97.578        14.389
==============================================

=====================================
Elastic modulus as a function of time
=====================================
#Time, stretch, twist, stretch-twist
1000.000, 991.91168, 1373.35825, 200.50841
2000.000, 1004.78021, 516.21083, 3.60912
3000.000, 931.30826, 509.38341, -16.08887
.
.
.
118000.000, 1081.86818, 444.95480, -98.78361
119000.000, 1081.20986, 444.61885, -98.34127
120000.000, 1080.37979, 442.49438, -97.57808

Global deformation free energy

To caluclate global deformation enrgy, dnaMD globalEnergy can be used. At first, elastic matrix from reference DNA (most often free or unbound DNA) is calculated and subsequently this matrix is used to calculate deformation free energy of probe DNA (most often bound DNA).

The deformation free energy is calculated using elastic matrix as follows

\[G = \frac{1}{2L_0}\mathbf{xKx^T}\]
\[\mathbf{x} = \begin{bmatrix} (\theta^{x} - \theta^{x}_0) & (\theta^{y} - \theta^{y}_0) & (L - L_0) & (\phi - \phi_0) \end{bmatrix}\]

Where, \(\mathbf{K}\), \(\theta^{x}_0\), \(\theta^{y}_0\), \(L_0\) and \(\phi_0\) is calculated from reference DNA while \(\theta^{x}\), \(\theta^{y}\), \(L\) and \(\phi\) is calculated for probe DNA from each frame.

This command gives output energy as a function of time in csv file and also average energies with error.

In [5]:
%%bash

dnaMD globalEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -estype BST -et "all" -paxis X -gt "gmx analyze" -em "block" -o energy_all_BST.csv

# Print first  and last 3 line of output file
echo "==========================================================="
echo "Deformation free energy of bound DNA as a function of time"
echo "==========================================================="
head -4 energy_all_BST.csv
printf ".\n.\n.\n"
tail -3 energy_all_BST.csv
==============================================
Energy(kJ/mol)         Average         Error
----------------------------------------------
full                    27.714         0.724
diag                    27.012         0.751
stretch                  2.468         0.157
twist                   15.061         0.478
st_coupling             -0.743         0.040
b1                       1.477         0.180
b2                       8.005         0.266
bend                     9.718         0.386
bs_coupling              0.034         0.002
bt_coupling              0.825         0.020
bb_coupling              0.235         0.020
st                      26.269         0.713
bs                      11.984         0.499
bt                      25.369         0.634
==============================================

===========================================================
Deformation free energy of bound DNA as a function of time
===========================================================
#Time, full, diag, stretch, twist, st_coupling, b1, b2, bend, bs_coupling, bt_coupling, bb_coupling, st, bs, bt
0.000, 29.76525, 29.12752, 3.05586, 16.01787, -0.96687, 0.27243, 9.78136, 10.18484, 0.05805, 1.09664, 0.13105, 28.16064, 13.16770, 27.16829
10.000, 27.69000, 26.82145, 1.78113, 13.65491, -0.68154, 3.00373, 8.38168, 11.78823, 0.02173, 0.69127, 0.40282, 26.13991, 13.18827, 25.73159
20.000, 32.65326, 31.73848, 3.36328, 15.22752, -0.98900, 0.93046, 12.21722, 13.41835, 0.05997, 1.11574, 0.27067, 30.74949, 16.57094, 29.49094
.
.
.
49980.000, 28.48093, 28.59762, 4.40662, 16.38502, -1.17429, 0.45068, 7.35530, 7.95214, 0.05514, 0.91465, 0.14616, 27.42332, 12.26774, 25.10564
49990.000, 33.18350, 32.91519, 3.56658, 21.03477, -1.19700, 0.05849, 8.25534, 8.36962, 0.06281, 1.21256, 0.05579, 31.71819, 11.94323, 30.56117
50000.000, 30.81292, 29.94787, 1.86541, 19.01796, -0.82313, 0.67654, 8.38796, 9.25575, 0.03664, 1.02777, 0.19124, 29.12474, 10.96655, 29.11023

Some of the plots from above data can be found here


Deformation free energy can be calculated as the following terms:

  • full : Use entire elastic matrix – all motions with their coupling
  • diag : Use diagonal of elastic matrix – all motions but no coupling
  • b1 : Only bending-1 motion
  • b2 : Only bending-2 motion
  • stretch : Only stretching motion
  • twist : Only Twisting motions
  • st_coupling : Only stretch-twist coupling motion
  • bs_coupling : Only Bending and stretching coupling
  • bt_coupling : Only Bending and Twisting coupling
  • bb_coupling : Only bending-1 and bending-2 coupling
  • bend : Both bending motions with their coupling
  • st : Stretching and twisting motions with their coupling
  • bs : Bending (b1, b2) and stretching motions with their coupling
  • bt : Bending (b1, b2) and twisting motions with their coupling

When all is used, all above terms were calculated.

In [6]:
%%bash

dnaMD globalEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 24 -estype ST -et "all" -gt "gmx analyze" -em "block"  -o energy_all_ST.csv

# Print first  and last 3 line of output file
echo "==========================================================="
echo "Deformation free energy of bound DNA as a function of time"
echo "==========================================================="
head -4 energy_all_ST.csv
printf ".\n.\n.\n"
tail -3 energy_all_ST.csv
==============================================
Energy(kJ/mol)         Average         Error
----------------------------------------------
full                    15.398         0.353
diag                    17.211         0.408
stretch                  2.581         0.137
twist                   14.630         0.362
st_coupling             -0.906         0.037
==============================================

===========================================================
Deformation free energy of bound DNA as a function of time
===========================================================
#Time, full, diag, stretch, twist, st_coupling
0.000, 24.65120, 27.93891, 4.01772, 23.92119, -1.64386
10.000, 16.45031, 18.42131, 2.11880, 16.30251, -0.98550
20.000, 22.50793, 25.46474, 3.54666, 21.91808, -1.47841
.
.
.
49980.000, 13.37563, 15.94228, 5.74292, 10.19936, -1.28332
49990.000, 17.67070, 20.02852, 2.88301, 17.14551, -1.17891
50000.000, 16.73218, 18.89026, 2.53138, 16.35888, -1.07904

Same as above but only for stretching and twisting motions.

Local elastic properties

Local elastic properties can be caluclated using either local base-step parameters or local helical base-step parameters.

In case of base-step parameters: Shift (\(Dx\)), Slide (\(Dy\)), Rise (\(Dz\)), Tilt (\(\tau\)), Roll (\(\rho\)) and Twist (\(\omega\)), following elastic matrix is calculated.

\[\begin{split}\mathbf{K}_{base-step} = \begin{bmatrix} K_{Dx} & K_{Dx,Dy} & K_{Dx,Dz} & K_{Dx,\tau} & K_{Dx,\rho} & K_{Dx,\omega} \\ K_{Dx,Dy} & K_{Dy} & K_{Dy,Dz} & K_{Dy,\tau} & K_{Dy,\rho} & K_{Dy,\omega} \\ K_{Dx,Dz} & K_{Dy,Dz} & K_{Dz} & K_{Dz,\tau} & K_{Dz,\rho} & K_{Dz,\omega} \\ K_{Dx,\tau} & K_{Dy,\tau} & K_{Dz,\tau} & K_{\tau} & K_{\tau, \rho} & K_{\tau,\omega} \\ K_{Dx,\rho} & K_{Dy,\rho} & K_{Dz,\rho} & K_{\tau, \rho} & K_{\rho} & K_{\rho,\omega} \\ K_{Dx,\omega} & K_{Dy,\omega} & K_{Dz,\omega} & K_{\tau, \omega} & K_{\rho, \omega} & K_{\omega} \\ \end{bmatrix}\end{split}\]

In case of helical-base-step parameters: x-displacement (\(dx\)), y-displacement (\(dy\)), h-rise (\(h\)), inclination (\(\eta\)), tip (\(\theta\)) and twist (\(\Omega\)), following elastic matrix is calculated.

\[\begin{split}\mathbf{K}_{helical-base-step} = \begin{bmatrix} K_{dx} & K_{dx,dy} & K_{dx,h} & K_{dx,\eta} & K_{dx,\theta} & K_{dx,\Omega} \\ K_{dx,dy} & K_{dy} & K_{dy,h} & K_{dy,\eta} & K_{dy,\theta} & K_{dy,\Omega} \\ K_{dx,h} & K_{dy,h} & K_{h} & K_{h,\eta} & K_{h,\theta} & K_{h,\Omega} \\ K_{dx,\eta} & K_{dy,\eta} & K_{h,\eta} & K_{\eta} & K_{\eta, \theta} & K_{\eta,\Omega} \\ K_{dx,\theta} & K_{dy,\theta} & K_{h,\theta} & K_{\eta, \theta} & K_{\theta} & K_{\theta,\Omega} \\ K_{dx,\Omega} & K_{dy,\Omega} & K_{h,\Omega} & K_{\eta, \Omega} & K_{\theta, \Omega} & K_{\Omega} \\ \end{bmatrix}\end{split}\]
In [7]:
%%bash

dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 7 -o local_elasticity_4-7bps.csv
=========== ============== Elastic Matrix =============== ===========

 150.32688   -12.74081    -4.99484    -0.86502    -0.08939    -0.16720
 -12.74081   129.99578    30.42056     0.27454    -0.58620    -2.06912
  -4.99484    30.42056   504.89222     0.65083    -0.06531     0.21751
  -0.86502     0.27454     0.65083     0.03433    -0.00104     0.00052
  -0.08939    -0.58620    -0.06531    -0.00104     0.01295     0.01427
  -0.16720    -2.06912     0.21751     0.00052     0.01427     0.05907

=========== ====================== ====================== ===========

Local elastic properties as a function of time

Same command can be used to calculate elasticity as a function of time with option -ot/--output-time and save it in csv format file. This result can be used to check their convergence.

  • If this option is used, -fgap/--frame-gap is an essential option.
  • The output file is in csv format and can be opened as spreadsheet.

NOTE: Elastic properties cannot be calculated using a single frame because fluctuations are required. Therefore, here time means trajectory between zero time to given time.

In [8]:
%%bash

dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 7 -fgap 200 -ot local_elasticity_time_4-7bps.csv

# Print first 3 and last 3 line of first seven columns from output file
echo "====================================="
echo "Elastic matrix as a function of time"
echo "====================================="
head -4 local_elasticity_time_4-8bps.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'
printf ".\n.\n.\n"
tail -3 local_elasticity_time_4-8bps.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'
=====================================
Elastic matrix as a function of time
=====================================
#Time, shift, slide, rise, tilt, roll, twist, ... ... ...
2000.000, 158.96347, 127.71194, 503.26390, 0.03143, 0.01542, 0.05937, ... ... ...
4000.000, 150.93805, 126.16929, 467.98254, 0.03050, 0.01409, 0.06100, ... ... ...
6000.000, 147.43936, 132.39932, 480.47463, 0.03256, 0.01358, 0.06143, ... ... ...
.
.
.
116000.000, 151.11253, 129.64049, 504.93262, 0.03417, 0.01293, 0.05897, ... ... ...
118000.000, 150.77496, 129.62910, 504.85789, 0.03426, 0.01296, 0.05910, ... ... ...
120000.000, 150.32688, 129.99578, 504.89222, 0.03433, 0.01295, 0.05907, ... ... ...

Local elastic properties of the consecutive overlapped DNA segments

Above method gives local elasticities of a small local segment of the DNA. However, we mostly interested in large segment of the DNA. This large segment can be further divided into smaller local segments. For these smaller segments local elasticities can be calculated. Here these segments overlapped with each other.

Same command can be used to calculate elasticity as a function of time with option -os/--output-segments and save it in csv format file.

  • If this option is used, -fgap/--frame-gap is an essential option.
  • The output file is in csv format and can be opened as spreadsheet.
In [9]:
%%bash

dnaMD localElasticity -i free_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -fgap 200 -gt "gmx analyze" -em "acf" -os local_elasticity_segments.csv

# Print first 3 and last 3 line of first seven columns from output file
echo "========================================"
echo "Elasticity as a function of DNA segments"
echo "========================================"
head -4 local_elasticity_segments.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'
printf ".\n.\n.\n"
tail -3 local_elasticity_segments.csv | awk '{print $1, $2, $3, $4, $5, $6 ,$7, "... ... ..."}'
========================================
Elasticity as a function of DNA segments
========================================
#bps, shift, shift-error, slide, slide-error, rise, rise-error, ... ... ...
4-7, 150.32688, 0.25654, 129.99578, 0.27262, 504.89222, 1.61893, ... ... ...
5-8, 162.67911, 0.46671, 122.65855, 0.47513, 510.18741, 1.36958, ... ... ...
6-9, 172.50587, 2.36981, 119.84961, 1.46698, 510.68380, 2.57983, ... ... ...
.
.
.
15-18, 153.98831, 2.02295, 137.21947, 0.29044, 546.15664, 1.31991, ... ... ...
16-19, 188.51349, 1.91188, 148.97261, 0.46969, 521.74643, 1.14698, ... ... ...
17-20, 174.28796, 2.07481, 134.04039, 0.67123, 573.20011, 2.54414, ... ... ...

Local deformation energy of a local small segment

At first, elastic matrix from reference DNA (most often free or unbound DNA) is calculated and subsequently this matrix is used to calculate deformation free energy of probe DNA (most often bound DNA).

\[G = \frac{1}{2}\mathbf{xKx^T}\]

When helical='False'

\[\mathbf{K} = \mathbf{K}_{base-step}\]
\[\mathbf{x} = \begin{bmatrix} (Dx_{i}-Dx_0) & (Dy_i - Dy_0) & (Dz_i - Dz_0) & (\tau_i - \tau_0) & (\rho_i - \rho_0) & (\omega_i - \omega_0) \end{bmatrix}\]

When helical='True'

\[\mathbf{K} = \mathbf{K}_{helical-base-step}\]
\[\mathbf{x} = \begin{bmatrix} (dx_{i}-dx_0) & (dy_i - dy_0) & (h_i - h_0) & (\eta_i - \eta_0) & (\theta_i - \theta_0) & (\Omega_i - \Omega_0) \end{bmatrix}\]
In [10]:
%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 10 -be 13 -et all -gt "gmx analyze" -em "block" -o local_energy_time_4-7bps.csv

# Print first 3 and last 3 line from output file
echo "==============================================="
echo "Local deformation energy as a function of time"
echo "==============================================="
head -4 local_energy_time_4-7bps.csv
printf ".\n.\n.\n"
tail -3 local_energy_time_4-7bps.csv
==============================================
Energy(kJ/mol)         Average         Error
----------------------------------------------
full                    19.143         0.671
diag                    43.865         1.839
shift                    1.366         0.064
slide                   15.994         1.069
rise                     1.194         0.030
tilt                     1.705         0.040
roll                     6.365         0.365
twist                   17.242         0.781
==============================================

===============================================
Local deformation energy as a function of time
===============================================
#Time, full, diag, shift, slide, rise, tilt, roll, twist
0.000, 32.99568, 97.62372, 0.85609, 43.00855, 0.55182, 0.05174, 8.64099, 44.51453
10.000, 24.91292, 63.21362, 0.18680, 31.85144, 0.31550, 2.74114, 5.80446, 22.31428
20.000, 34.19821, 78.73263, 0.77440, 31.54860, 4.31769, 0.19352, 8.74777, 33.15064
.
.
.
49980.000, 17.43224, 50.92180, 0.07917, 12.95852, 2.46789, 0.92035, 5.74851, 28.74736
49990.000, 18.86625, 58.72850, 3.25713, 25.11722, 0.01483, 1.27080, 1.81182, 27.25671
50000.000, 19.01918, 49.56285, 2.16952, 13.54590, 1.86029, 4.52988, 4.38336, 23.07391

Some of the plots from above data can be found here


Same as the above but energy is calculated using helical base-step parameters

In [11]:
%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 10 -be 13 -et all -helical -gt "gmx analyze" -em "block" -o local_helical_energy_time_4-7bps.csv

# Print first 3 and last 3 line from output file
echo "======================================================"
echo "Local helical deformation energy as a function of time"
echo "======================================================"
head -4 local_energy_time_4-7bps.csv
printf ".\n.\n.\n"
tail -3 local_energy_time_4-7bps.csv
==============================================
Energy(kJ/mol)         Average         Error
----------------------------------------------
full                    18.541         0.607
diag                    82.103         3.204
x-disp                  38.632         1.974
y-disp                   0.940         0.024
h-rise                   6.768         0.333
inclination             15.242         0.857
tip                      1.288         0.030
h-twist                 19.232         0.825
==============================================

======================================================
Local helical deformation energy as a function of time
======================================================
#Time, full, diag, shift, slide, rise, tilt, roll, twist
0.000, 32.99568, 97.62372, 0.85609, 43.00855, 0.55182, 0.05174, 8.64099, 44.51453
10.000, 24.91292, 63.21362, 0.18680, 31.85144, 0.31550, 2.74114, 5.80446, 22.31428
20.000, 34.19821, 78.73263, 0.77440, 31.54860, 4.31769, 0.19352, 8.74777, 33.15064
.
.
.
49980.000, 17.43224, 50.92180, 0.07917, 12.95852, 2.46789, 0.92035, 5.74851, 28.74736
49990.000, 18.86625, 58.72850, 3.25713, 25.11722, 0.01483, 1.27080, 1.81182, 27.25671
50000.000, 19.01918, 49.56285, 2.16952, 13.54590, 1.86029, 4.52988, 4.38336, 23.07391

Deformation energy of the consecutive overlapped DNA segments

Above method gives energy of a small local segment of the DNA. However, we mostly interested in large segment of the DNA. This large segment can be further divided into smaller local segments. For these smaller segments local deformation energy can be calculated. Here these segments overlapped with each other.

In [12]:
%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -et all -gt "gmx analyze" -em "block" -os local_energy_segments.csv

# Print first 3 and last 3 line from output file
echo "=================================================="
echo "Local deformation energy as a function of segments"
echo "=================================================="
head -4 local_energy_segments.csv
printf ".\n.\n.\n"
tail -3 local_energy_segments.csv
==================================================
Local deformation energy as a function of segments
==================================================
#bps, full, full-error, diag, diag-error, shift, shift-error, slide, slide-error, rise, rise-error, tilt, tilt-error, roll, roll-error, twist, twist-error
4-7, 10.44541, 0.15599, 16.11461, 0.50862, 0.95000, 0.02392, 8.97490, 0.44581, 1.17289, 0.13403, 1.71030, 0.05021, 1.08986, 0.03874, 2.21666, 0.17392
5-8, 9.49626, 0.25053, 19.32846, 0.70809, 1.14060, 0.02737, 9.59726, 0.51038, 1.21808, 0.18759, 1.75363, 0.04521, 1.11454, 0.05955, 4.50435, 0.19880
6-9, 10.52567, 0.28519, 20.09691, 0.54228, 1.22326, 0.02970, 11.45709, 0.43514, 1.51242, 0.08500, 1.90555, 0.04658, 1.31208, 0.11407, 2.68650, 0.09203
.
.
.
15-18, 13.98537, 0.52298, 41.05547, 2.63829, 1.04577, 0.02497, 17.01431, 0.86024, 2.44413, 0.14371, 1.78045, 0.04873, 1.78854, 0.20280, 16.98227, 2.09869
16-19, 14.04847, 0.37844, 39.34525, 1.49491, 1.97148, 0.09511, 14.20378, 0.94662, 2.52573, 0.08194, 1.70940, 0.04099, 0.98835, 0.09072, 17.94652, 1.15043
17-20, 7.89839, 0.32320, 14.91991, 0.78938, 1.03028, 0.02583, 5.81921, 0.49549, 1.40673, 0.06566, 1.46641, 0.03364, 1.05240, 0.14133, 4.14489, 0.38861

Some of the plots from above data can be found here


Same as the above but energy is calculated using helical base-step parameters

In [13]:
%%bash

dnaMD localEnergy -ir free_dna.h5 -ip bound_dna.h5 -tbp 27 -bs 4 -be 20 -span 4 -et all -helical -gt "gmx analyze" -em "block" -os local_helical_energy_segments.csv

# Print first 3 and last 3 line from output file
echo "=========================================================="
echo "Local helical deformation energy as a function of segments"
echo "=========================================================="
head -4 local_helical_energy_segments.csv
printf ".\n.\n.\n"
tail -3 local_helical_energy_segments.csv
==========================================================
Local helical deformation energy as a function of segments
==========================================================
#bps, full, full-error, diag, diag-error, x-disp, x-disp-error, y-disp, y-disp-error, h-rise, h-rise-error, inclination, inclination-error, tip, tip-error, h-twist, h-twist-error
4-7, 9.94039, 0.17649, 17.73868, 0.51652, 9.00892, 0.37199, 1.19946, 0.04468, 1.54174, 0.13701, 1.69527, 0.06446, 1.73424, 0.04823, 2.55906, 0.16107
5-8, 9.01980, 0.22539, 19.49642, 0.57573, 8.19696, 0.26473, 1.36335, 0.06127, 1.48339, 0.27526, 1.60444, 0.08879, 1.76889, 0.04544, 5.07939, 0.32815
6-9, 10.30985, 0.26888, 17.00383, 0.38032, 7.38157, 0.25950, 1.96202, 0.28106, 0.99682, 0.02458, 1.68845, 0.14166, 1.89859, 0.04606, 3.07639, 0.16278
.
.
.
15-18, 13.74199, 0.55739, 53.29447, 4.60218, 27.45388, 1.71961, 0.66820, 0.01620, 1.17571, 0.02924, 3.37921, 0.40760, 1.29946, 0.03306, 19.31802, 2.43431
16-19, 12.60741, 0.33640, 46.34514, 2.01539, 19.89670, 0.90439, 0.90573, 0.04403, 1.01976, 0.02504, 1.68335, 0.09711, 1.12857, 0.06776, 21.71103, 1.22043
17-20, 7.54442, 0.33087, 15.27288, 0.76005, 5.60088, 0.40791, 0.91852, 0.02377, 1.13807, 0.02707, 1.66290, 0.20744, 1.25180, 0.02884, 4.70071, 0.43934