Sphere Model: Langevin Dynamics¶
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
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 freezed all the lamina beads, and would not consider the dynamics of the membrane, and that is why we set “False” (off) for membrane dynamics. Consequently, there was also no need to set the bond between lamina beads, so we set “None” for the variable “lam_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 = "initial_configs/human_1.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 = False, lam_bond = None)
Add the force field¶
In this example, we loaded default settings for the interactions between chromosomes and chromosomes, and chromosomes and nuclear landmarks. All the types of potential can be found in “Section: Energy function of the whole nucleus model” in Supporting Information. According to the order of added potential, the index of speckle-speckle potential is 6, and this index is needed when transiting the speckle particle between type 6 and type 7.
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
model.load_default_settings()
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 #Currently, the number of speckle is 1600.
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 95731.85929571482 86345.68982083745 110.69893791415754 --
4.0% 4000 99532.94729061538 91125.1333946349 116.82639289816234 12:18
6.0% 6000 102489.94966086837 92875.1373802856 119.0699742854793 12:07
8.0% 8000 101970.62104604443 93946.25490129643 120.44319363439149 11:51
10.0% 10000 104080.58569641429 93895.10736796792 120.37762026728836 11:36
12.0% 12000 102469.63879127387 94266.58652759668 120.85387274165592 11:20
14.0% 14000 103107.25368497109 94476.86194657334 121.12345499409341 11:05
16.0% 16000 101997.27652061332 93615.4220947793 120.01905155637233 10:49
18.0% 18000 102695.14047996342 94293.17782819153 120.88796394805793 10:34
20.0% 20000 102284.79145205846 94000.65665731127 120.51293905682084 10:18
22.0% 22000 103181.19623978612 94092.67913100094 120.63091587908977 10:03
24.0% 24000 102471.11959102858 93940.2788685578 120.43553209994224 9:47
26.0% 26000 102967.84568309924 94193.30678616917 120.75992491908966 9:32
28.0% 28000 102104.78401933674 94638.44177388668 121.33060737551574 9:16
30.0% 30000 103145.97248266514 93784.43750129726 120.23573667443081 9:01
32.0% 32000 102815.6439428291 93162.17754139185 119.43797228328442 8:45
34.0% 34000 103375.09306618437 94517.692018082 121.17580092543459 8:30
36.0% 36000 103472.25183194957 93685.96252558767 120.1094874633261 8:15
38.0% 38000 104177.71982035965 94111.53321030273 120.65508763055053 7:59
40.0% 40000 102929.15471019791 94273.22366385003 120.86238183962666 7:44
42.0% 42000 103807.91912237357 93392.1141090396 119.73276098529199 7:28
44.0% 44000 102317.0566719379 93764.79770755101 120.2105576028052 7:13
46.0% 46000 103016.71835973949 93875.29215481882 120.35221630034223 6:57
48.0% 48000 102693.02057004237 94267.37592089028 120.85488477826573 6:42
50.0% 50000 104012.35688091657 94484.31101151256 121.13300502002912 6:26
52.0% 52000 102642.30455152449 95044.68483915059 121.85142869216766 6:11
54.0% 54000 103806.89709221238 94755.90310995904 121.48119792816715 5:55
56.0% 56000 103170.73528723426 94183.26523600474 120.74705121409615 5:40
58.0% 58000 103770.68482632717 94051.89853501321 120.57863337752292 5:24
60.0% 60000 102732.78325580872 93822.51084952969 120.2845483663741 5:09
62.0% 62000 103234.22757751837 93808.19539673388 120.26619533191987 4:53
64.0% 64000 102881.23494367936 93540.57531709578 119.92309472502322 4:38
66.0% 66000 102959.45930786872 94174.98848670874 120.73644006073863 4:23
68.0% 68000 102207.69833765701 94185.86471289878 120.75038385677814 4:07
70.0% 70000 103380.66204900964 94032.35585227858 120.55357880644698 3:52
72.0% 72000 102557.51791599748 94209.38874932617 120.78054269682411 3:36
74.0% 74000 103333.50519057558 93847.19960158567 120.31620042262259 3:21
76.0% 76000 102915.37827298888 93417.95780964197 119.76589373591789 3:05
78.0% 78000 103285.69125801143 93846.46844148931 120.31526304350987 2:50
80.0% 80000 102317.42949138349 93657.30623426713 120.07274884880665 2:34
82.0% 82000 103100.49753877813 93921.214984692 120.41109136988563 2:19
84.0% 84000 102898.78388727782 93804.70340561948 120.26171844710522 2:03
86.0% 86000 103186.26768478611 94014.1095649399 120.53018627077486 1:48
88.0% 88000 102891.62793189581 94033.03357096677 120.55444767134507 1:32
90.0% 90000 103181.6413625974 93624.9302086184 120.0312413727458 1:17
92.0% 92000 102430.7777392424 94168.32895799648 120.72790225465954 1:01
94.0% 94000 103673.4257128592 93934.5022553975 120.42812623008454 0:46
96.0% 96000 102454.2311218296 94226.1241681962 120.80199823326495 0:30
98.0% 98000 103528.39761667955 93969.1743303926 120.47257734148573 0:15
100.0% 100000 102507.23206314092 94393.40130067056 121.0164548082331 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')