Two following input files are required in this tutorial
tutorial_data/elasticity_DNA/free_dna.h5
tutorial_data/elasticity_DNA/bound_dna.h5
These two files should be present inside tutorial_data/elasticity_DNA of the present working directory.
In [1]:
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import dnaMD
%matplotlib inline
eyDNA object is initialized by using the total number of base-pairs and HDF5 file.
This class contains all the required functions to calculate the elastic properties and deformation free energy.
In [2]:
eyDNA = dnaMD.dnaEY(27, 'BST', filename='elasticity_DNA/free_dna.h5')
In [3]:
# All frames
avg, mod_matrix = eyDNA.getStretchTwistBendModulus([4,20], paxis='X')
print('Average values for all frames: ', avg)
print('Modulus matrix for all frames: \n', mod_matrix )
print(' ')
# Elastic matrix
avg, mod_matrix = eyDNA.getStretchTwistBendModulus([4,20], paxis='X', matrix=True)
print('Average values for all frames: ', avg)
print('Elastic constant matrix for all frames: \n', mod_matrix )
print(' ')
The elastic matrix is in this form:
$$\text{Elastic matrix} = \begin{bmatrix} K_{Bx} & K_{Bx,By} & K_{Bx,S} & K_{Bx,T} \\ K_{Bx,By} & K_{By} & K_{By,S} & K_{By,T} \\ K_{Bx,S} & K_{By,S} & K_{S} & K_{S,T} \\ K_{Bx,T} & K_{Bx,T} & K_{S,T} & K_{T} \end{bmatrix} $$Where:
Where:
Elasticities cannot be calcualted from an individual snapshot or frame. However, these properties can be calculated as a function of time by considering all the frames up to that time. For example, 0-50 ns, 0-100 ns, 0-150 ns etc. By this method, we can analyze the convergence and also further we can calculate error using block average method.
Elasticities over the time can be calculated using getElasticityByTime method.
If esType='BST'
, A ordered dictionary of 1D arrays of shape (nframes). The keys in dictionary are name of the elasticity in the same order as listed above..
bend-1
- Bending-1 stiffness in one planebend-2
- Bending-2 stiffness in another orthogonal planestretch
- Stretch Modulustwist
- Twist rigiditybend-1-bend-2
- Bending-1 and Bending-2 couplingbend-2-stretch
- Bending-2 and stretching couplingstretch-twist
- Stretching Twsiting couplingbend-1-stretch
- Bending-1 Stretching couplingbend-2-twist
- Bending-2 Twisting couplingbend-1-twist
- Bending-1 Twisting couplingIf esType='ST'
, 2D array with three properties of shape (3, frame) will be returned.
stretch
- Stretch Modulustwist
- Twist rigiditystretch-twist
- Stretching Twsiting couplingIn the following example, modulus as a function of time was calculated by adding 1000 frames.
In [4]:
time, modulus = eyDNA.getModulusByTime([4,20], frameGap=500, masked=True)
print('Keys in returned dictionary:\n', '\n'.join(list(modulus.keys())), '\n-----------')
# Stretching modulus
plt.plot(time, modulus['stretch'])
plt.scatter(time, modulus['stretch'])
plt.xlabel('Time (ps)')
plt.ylabel(r'Stretching Modulus (pN)')
plt.show()
# Twist rigidity
plt.plot(time, modulus['twist'])
plt.scatter(time, modulus['twist'])
plt.xlabel('Time (ps)')
plt.ylabel(r'Rigidity (pN nm$^2$)')
plt.show()
# Stretch twist coupling
plt.plot(time, modulus['stretch-twist'])
plt.scatter(time, modulus['stretch-twist'])
plt.xlabel('Time (ps)')
plt.ylabel(r'Stretch-Twist Coupling (pN nm)',)
plt.show()
Deformation energy of a probe DNA (bound DNA) can be calculated with reference to the DNA present in the current object.
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.
We already loaded the data for reference DNA above. Here, we will load data for probe DNA.
In [5]:
# Load parameters of bound DNA
boundDNA = dnaMD.DNA(27, filename='elasticity_DNA/bound_dna.h5')
Deformation free energy can be calculated for the following motions that can be used with which
option.
'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 couplingwhich
can be either 'all'
or a list of energy terms given above.
In [6]:
# Deformation free energy of bound DNA and calculate all above listed terms
time, energy = eyDNA.getGlobalDeformationEnergy([4,20], boundDNA, paxis='X', which='all', masked=True)
energyTerms=list(energy.keys())
print('Keys in returned dictionary:\n', '\n'.join(energyTerms), '\n-----------')
# Plot two energy terms
fig = plt.figure(figsize=(8,8))
fig.subplots_adjust(hspace=0.3)
ax1 = fig.add_subplot(211)
ax1.set_title('Bound DNA, entire elastic matrix')
ax1.plot(time, energy['full'])
ax1.set_xlabel('Time (ps)')
ax1.set_ylabel(r'Deformation Free Energy (kJ/mol)',)
ax2 = fig.add_subplot(212)
ax2.set_title('Bound DNA, only diagonal of elastic matrix')
ax2.plot(time, energy['diag'])
ax2.set_xlabel('Time (ps)')
ax2.set_ylabel(r'Deformation Free Energy (kJ/mol)',)
plt.show()
# Calculate average and error for each energy terms
error = dnaMD.get_error(time, list(energy.values()), len(energyTerms), err_type='block', tool='gmx analyze')
print("==============================================")
print('{0:<16}{1:>14}{2:>14}'.format('Energy(kJ/mol)', 'Average', 'Error'))
print("----------------------------------------------")
for i in range(len(energyTerms)):
print('{0:<16}{1:>14.3f}{2:>14.3f}'.format(energyTerms[i], np.mean(energy[energyTerms[i]]),error[i]))
print("==============================================\n")
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.
$$ \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} $$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.
$$ \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} $$
In [7]:
# base-step
avg, matrix = eyDNA.calculateLocalElasticity([10,13], helical=False)
# Print matrix in nice format
out = ''
mean_out = ''
for i in range(matrix.shape[0]):
for j in range(matrix.shape[0]):
if j != matrix.shape[0]-1:
out += '{0:>10.5f} '.format(matrix[i][j])
else:
out += '{0:>10.5f}\n'.format(matrix[i][j])
mean_out += '{0:>15.3f} '.format(avg[i])
print('Average values for all frames: ', mean_out)
print('=========== ============== Elastic Matrix =============== ===========\n')
print(out)
print('=========== ====================== ====================== ===========')
# helical base-step
avg, matrix = eyDNA.calculateLocalElasticity([10,13], helical=True)
# Print matrix in nice format
out = ''
mean_out = ''
for i in range(matrix.shape[0]):
for j in range(matrix.shape[0]):
if j != matrix.shape[0]-1:
out += '{0:>10.5f} '.format(matrix[i][j])
else:
out += '{0:>10.5f}\n'.format(matrix[i][j])
mean_out += '{0:>15.3f} '.format(avg[i])
print('\n\nAverage values for all frames: ', mean_out)
print('=========== ============== Elastic Matrix =============== ===========\n')
print(out)
print('=========== ====================== ====================== ===========')
In [8]:
# Here calculate energy for one base-step
time, energy = eyDNA.getLocalDeformationEnergy([10,13], boundDNA, helical=False, which='all')
energyTerms=list(energy.keys())
print('Keys in returned dictionary:\n', '\n'.join(energyTerms), '\n-----------')
# Plot two energy terms
fig = plt.figure(figsize=(8,8))
fig.subplots_adjust(hspace=0.3)
ax1 = fig.add_subplot(211)
ax1.set_title('Bound DNA, entire elastic matrix')
ax1.plot(time, energy['full'])
ax1.set_xlabel('Time (ps)')
ax1.set_ylabel(r'Local Deformation Energy (kJ/mol)',)
ax2 = fig.add_subplot(212)
ax2.set_title('Bound DNA, only diagonal of elastic matrix')
ax2.plot(time, energy['diag'])
ax2.set_xlabel('Time (ps)')
ax2.set_ylabel(r'Local Deformation Energy (kJ/mol)',)
plt.show()
# Calculate average and error for each energy terms
error = dnaMD.get_error(time, list(energy.values()), len(energyTerms), err_type='block', tool='gmx analyze')
print("==============================================")
print('{0:<16}{1:>14}{2:>14}'.format('Energy(kJ/mol)', 'Average', 'Error'))
print("----------------------------------------------")
for i in range(len(energyTerms)):
print('{0:<16}{1:>14.3f}{2:>14.3f}'.format(energyTerms[i], np.mean(energy[energyTerms[i]]),error[i]))
print("==============================================\n")
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 [9]:
# First calculation for local base-step parameters
segments, energies, error = eyDNA.getLocalDeformationEnergySegments([4,20], boundDNA, span=4,
helical=False, which='all',
err_type='block',
tool='gmx analyze')
energyTerms=list(energies.keys())
print('Keys in returned dictionary:\n', '\n'.join(energyTerms), '\n-----------')
# Now plot the data
fig = plt.figure(figsize=(14,8))
fig.subplots_adjust(hspace=0.3)
mpl.rcParams.update({'font.size': 16})
xticks = range(len(segments))
ax1 = fig.add_subplot(111)
ax1.set_title('Local base-step parameters')
for term in energyTerms:
ax1.errorbar(xticks, energies[term], yerr=error[term], ms=10, elinewidth=3, fmt='-o', label=term)
ax1.set_xticks(xticks)
ax1.set_xticklabels(segments, rotation='vertical')
ax1.set_xlabel('base-step number')
ax1.set_ylabel(r'Deformation Energy (kJ/mol)',)
plt.legend()
plt.show()
Same as the above but energy is calculated using helical base-step parameters
In [10]:
# Secind calculation for local base-step parameters
segments, energies, error = eyDNA.getLocalDeformationEnergySegments([4,20], boundDNA, span=4,
helical=True, which='all',
err_type='block',
tool='gmx analyze')
energyTerms=list(energies.keys())
print('Keys in returned dictionary:\n', '\n'.join(energyTerms), '\n-----------')
# Now plot the data
fig = plt.figure(figsize=(14,8))
fig.subplots_adjust(hspace=0.3)
mpl.rcParams.update({'font.size': 16})
xticks = range(len(segments))
ax1 = fig.add_subplot(111)
ax1.set_title('Local base-step parameters')
for term in energyTerms:
ax1.errorbar(xticks, energies[term], yerr=error[term], ms=10, elinewidth=3, fmt='-o', label=term)
ax1.set_xticks(xticks)
ax1.set_xticklabels(segments, rotation='vertical')
ax1.set_xlabel('base-step number')
ax1.set_ylabel(r'Deformation Energy (kJ/mol)',)
plt.legend()
plt.show()