Nuclear Deformation Simulation¶
Import the necessary packages¶
[1]:
import parmed as pmd
import json
import sys
from sys import platform
import mdtraj as md
import simtk.openmm.app as mmapp
import mdtraj.reporters
import simtk.unit as u
import random
from openNucleome import OpenNucleome
from openNucleome.utils import coor_transformation, final_frame
import warnings
import numpy as np
import pandas as pd
Important parameters¶
We defined the important parameters in the next block. First, we set the transition probability between dP particles and P particles as 0.2, and the transition frequency as 4000. When creating the system, we set type 6, 7 as the dP and P particles, respectively, so we kept using 6, 7 here. In this example, we ran a simulation with total length of 100,000 steps, and output one configuration and the energy every 2000 steps.
[2]:
prob_P_dP = 0.2 # Transition probability from P to dP
prob_dP_P = 0.2 # Transition probability from dP to P
transition_freq = 4000
sampling_freq = 2000
dP_type = 6
P_type = 7
total_steps = 100000
Initialize the system¶
We first set up an example “model” of class “OpenNucleome” with the conserved temperature, damping coefficient, timestep, and the mass scale. In this folder, we also included the initial configuration “human.pdb” used for the simulation and created a system according to the initial configuration.
In this example, we moved all the lamina beads, and would consider the dynamics of the membrane, and that is why we set “True” (on) for membrane dynamics. Consequently, we need to designate the bonds between specific lamina beads, so we included a text file, which logs all the pairs linked by the bonds, for the variable “membrane_bond”.
[3]:
model = OpenNucleome(1.0, 0.1, 0.005, 1.0) # 1.0: temperature (LJ reduced unit);
# 0.1: damping coefficient (LJ reduced unit);
# 0.005: timestep (LJ reduced unit);
# 1.0: mass_scale
PDB_file = "strength_one/human.pdb" #The initial configuration
# Generate new elements and construct topology as well
# flag_membrane: True for including lamina dynamics, False for excluding lamina dynamics;
# lam_bond: A file contains the lamina bond when membrane_dynamics is on.
warnings.filterwarnings("ignore")
model.create_system(PDB_file, flag_membrane = True, lam_bond = 'input_params/lamina_bond.txt')
Add the force field¶
Different from the sphere nucleus situation, we used a csv file to log all the flags of interactions between chromosomes and chromosomes, or chromosomes and nuclear landmarks.
The column from “bond” to “lam_squeeze” is boolen, which means the switches of the corresponding interactions. For example, if the spec-chrom is True, it means in the simulation, the interaction between speckle and chromosome will be present; if the inter is False, it means we exclude the interchromosomal interactions between chromosome. Remember some potentials require the corresponding parameter files, such as ideal requires ideal_param_file and spec-chrom requires the chr_spec_param, so if you turn those potentials on, you should also include the corresponding parameter files.
In this example, we turn every interaction on, and set the strength of force squeezing the nucleus as 1.0 (k) when loading the customized force field. The users can set their specific strengths.
Because during the simulations, we would transit the dP speckle and P speckle particles, here, we logged the start index and end index of speckle particle, and computed the number of speckle particles.
[4]:
# Add the default force field
# Add the customized force field
# In this customized setting, I logged all the flags for different interactions
# (True means on, False means off, and NaN means N/A)
# and all the corresponding files required for interactions.
force_field = pd.read_csv('input_params/input.csv', sep=' ', header=0, index_col=0)
model.load_customized_settings(force_field, 'input_params', k = 1.0)
index_spec_spec_potential = 6 #The 6-th potential is speckle-speckle interaction in the default settings
start_spec_index = model.N_chr_nuc+1
end_spec_index = model.N_chr_nuc_spec+1
N_spec = end_spec_index-start_spec_index
force_field
[4]:
| bond | angle | softcore | ideal | compt | inter | spec-spec | spec-chrom | nuc-nuc | nuc-spec | ... | lam-chrom | hard-wall | lam-lam | lam_squeeze | ideal_param_file | compt_param_file | interchr_param_file | chr_spec_param | chr_nuc_param | chr_lam_param | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| chromosome | True | True | True | True | True | True | NaN | NaN | NaN | NaN | ... | NaN | NaN | NaN | NaN | ideal_param_file.txt | compt_param_file.txt | interchr_param_file.txt | NaN | NaN | NaN |
| speckle | NaN | NaN | NaN | NaN | NaN | NaN | True | True | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | chr_spec_param.txt | NaN | NaN |
| nucleolus | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | True | True | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | chr_nuc_param.txt | NaN |
| lamina | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | ... | True | True | True | True | NaN | NaN | NaN | NaN | NaN | chr_lam_param.txt |
4 rows × 21 columns
Perform the simulation¶
We first created the simulation with a specific Platform, and here, we used “CUDA” but users can also use “CPU”, “Reference”, and “OpenCL” according to their hardware. Before the simulation, we minimized the energy to make the system much more reasonable and stable. After randomly setting velocity, we started our simulation with a total length of 100,000 steps and output the configuration and energy every 2000 steps, and change the speckle types every 4000 steps as we mentioned in the previous blocks.
[5]:
#model.save_system("model_before_simulation_0.xml")
simulation = model.create_simulation(platform_type = "CUDA") # Users can also use CPU, Reference, OpenCL.
simulation.context.setPositions(model.chr_positions)
simulation.minimizeEnergy()
simulation.reporters.append(mdtraj.reporters.DCDReporter('results/step_100000.dcd', sampling_freq))
def setVelocity(context):
sigma = u.sqrt(1.0*u.kilojoule_per_mole / model.chr_system.getParticleMass(1))
velocs = u.Quantity(1.0 * np.random.normal(size=(model.chr_system.getNumParticles(), 3)),
u.meter) * (sigma / u.meter)
context.setVelocities(velocs)
setVelocity(simulation.context)
simulation.reporters.append(mmapp.statedatareporter.StateDataReporter(sys.stdout, sampling_freq, step=True,
potentialEnergy=True, kineticEnergy=True, temperature=True, progress=True,
remainingTime=True, separator='\t', totalSteps = total_steps))
for i in range(total_steps//transition_freq):
simulation.step(transition_freq)
# Change the type of speckles every 4000 steps, non-equilibrium scheme.
# Do the chemical modification, and change the spec-spec potential on the fly.
for j in np.random.randint(start_spec_index-1, end_spec_index-1, N_spec):
if model.compart_type[j] == dP_type-1:
model.compart_type[j] = P_type-1 if random.random() < prob_dP_P else dP_type-1
else:
model.compart_type[j] = dP_type-1 if random.random() < prob_P_dP else P_type-1
# Update the context after changing the type of speckles.
for m in range(model.chr_system.getNumParticles()):
model.chr_system.getForce(index_spec_spec_potential).setParticleParameters(m, [model.compart_type[m]])
model.chr_system.getForce(index_spec_spec_potential).updateParametersInContext(simulation.context)
# Keep the final result of bead types in case constructing the configuration for the continuous simulation.
np.savetxt('results/compt_final_frame.txt', (np.array(model.compart_type)+1).reshape((-1,1)), fmt='%d')
#"Progress (%)" "Step" "Potential Energy (kJ/mole)" "Kinetic Energy (kJ/mole)" "Temperature (K)" "Time Remaining"
2.0% 2000 123136.44445905171 100752.69234185581 114.52055593216552 --
4.0% 4000 138145.93252829043 99854.70038122084 113.4998532971868 9:06
6.0% 6000 142575.6182367836 104465.90880812544 118.741183730089 8:38
8.0% 8000 141094.42452352808 107327.50563145665 121.99381799170929 8:21
10.0% 10000 139697.16552712512 107342.47070395519 122.01082804258003 8:13
12.0% 12000 137654.04464475974 106994.70033531352 121.61553482482218 8:03
14.0% 14000 138386.55393252496 105782.62134612219 120.2378251433289 7:53
16.0% 16000 137409.19082230702 106199.98460544148 120.71222112592687 7:42
18.0% 18000 138944.10917304258 106093.39691125305 120.59106821467195 7:32
20.0% 20000 137593.22835381792 106262.0393691768 120.78275567816587 7:21
22.0% 22000 138695.92867257126 105989.21260070454 120.47264710962654 7:10
24.0% 24000 138082.96422167047 105636.71725733916 120.07198324896527 6:59
26.0% 26000 138577.68014618923 105942.89960543014 120.42000543978695 6:48
28.0% 28000 138547.97569732778 105614.2894846135 120.04649071927538 6:37
30.0% 30000 139048.79707005448 105953.13915726854 120.43164422721559 6:26
32.0% 32000 139369.4478256774 106096.93739319986 120.5950925037643 6:15
34.0% 34000 139349.70734011347 105551.45826302844 119.9750736061584 6:04
36.0% 36000 138652.69678783504 106307.39016938012 120.83430366889978 5:53
38.0% 38000 138994.0196299976 105924.91785091767 120.39956647706012 5:42
40.0% 40000 138711.95389602106 106045.77361075555 120.53693718629414 5:31
42.0% 42000 139755.55680584465 106157.32762501454 120.66373506566056 5:20
44.0% 44000 139418.1188406263 105861.9173488944 120.32795695131931 5:09
46.0% 46000 139318.77307771594 106084.84309645879 120.58134552039144 4:58
48.0% 48000 139156.1582248191 106000.44885852029 120.48541880298527 4:47
50.0% 50000 139394.24871664512 106154.93144347834 120.66101144572286 4:36
52.0% 52000 138976.02488041663 105917.33440921981 120.39094675743841 4:24
54.0% 54000 139636.7751483859 105824.53815483821 120.2854698873599 4:13
56.0% 56000 139109.49119912702 105962.1255783129 120.44185864342542 4:02
58.0% 58000 139591.83060822415 105548.44699685683 119.9716508502199 3:51
60.0% 60000 139158.00653002597 105773.10164969457 120.22700457961139 3:40
62.0% 62000 139141.23690284294 106340.19938493858 120.87159626642863 3:29
64.0% 64000 138657.07116640714 106162.69666414094 120.66983778441957 3:18
66.0% 66000 138872.56066384647 105709.55038510424 120.154769029554 3:07
68.0% 68000 139552.46901895298 106237.45014421445 120.75480633361563 2:56
70.0% 70000 140314.6318439555 105998.06802734527 120.48271263103683 2:45
72.0% 72000 139878.4964976344 105824.37540777247 120.28528490089441 2:34
74.0% 74000 140049.54976571468 105978.09999204957 120.46001596209643 2:23
76.0% 76000 139040.28804717027 105698.80343463195 120.14255350742695 2:12
78.0% 78000 139441.61604379857 105796.50827126215 120.25360971791694 2:01
80.0% 80000 139103.81283076043 105733.9975154245 120.1825568622178 1:50
82.0% 82000 139735.49375639576 105520.89870734174 119.94033808470546 1:39
84.0% 84000 139625.44875567633 105857.33796241746 120.32275179130842 1:28
86.0% 86000 139302.36985220472 106414.25964148348 120.95577686304458 1:17
88.0% 88000 139197.9629622838 105156.448308775 119.52608550957876 1:06
90.0% 90000 138969.3590577207 105676.13902677708 120.1167920063982 0:55
92.0% 92000 138261.92996336846 105622.26065619843 120.05555115210468 0:44
94.0% 94000 138923.98136101887 106060.5908358773 120.55377918643204 0:33
96.0% 96000 138863.37451388055 105655.18511185612 120.0929747373017 0:22
98.0% 98000 138649.21396001056 106102.68049923926 120.60162040571038 0:11
100.0% 100000 138531.91066385998 106315.62779116184 120.84366696236692 0:00
Post-process the output files¶
We converted the trajectory with the default OpenMM unit to the nex trajectory with the reduced unit, and saved the final configuration according to the last frame.
[6]:
coor_transformation('results/step_100000.dcd', 'results/reduced_traj.dcd')
final_frame('results/step_100000.dcd', 'results/compt_final_frame.txt', 'results/final_config.pdb')