Sphere Model: Brownian 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 (A larger damping coefficient for the Brownian Dynamics)¶
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, 100, 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 99669.37249189641 93484.31840286203 119.85097090889614 --
4.0% 4000 102879.7194162421 93650.98329636063 120.06464256681159 12:11
6.0% 6000 104829.51531281356 93560.6536394227 119.94883600940348 11:58
8.0% 8000 105119.10886859841 93785.34285541113 120.23689737786663 11:42
10.0% 10000 105999.76281040246 93984.31355069151 120.49198648178462 11:27
12.0% 12000 105999.28209041196 93586.10689548533 119.98146819311104 11:11
14.0% 14000 107031.57428997077 93881.41950974891 120.36007183643666 10:56
16.0% 16000 106615.33266635888 93369.33133546957 119.7035524765387 10:41
18.0% 18000 106557.47937161615 94058.32191220552 120.58686842706025 10:26
20.0% 20000 107007.09729298094 93408.50286975966 119.75377209087473 10:10
22.0% 22000 107926.75086174044 94138.0784156327 120.6891197408138 9:55
24.0% 24000 107391.15834572143 93151.26510476637 119.42398206388837 9:40
26.0% 26000 107501.08014869876 94274.11024334499 120.8635184731827 9:24
28.0% 28000 107569.04196787113 94384.1702547747 121.00462020503947 9:09
30.0% 30000 108045.97777092032 94601.02906711765 121.282642654724 8:54
32.0% 32000 107432.52383930977 93589.96460811926 119.9864139488545 8:39
34.0% 34000 108329.11488470988 93621.14913681203 120.02639387388383 8:23
36.0% 36000 107958.7920328365 93807.04607312058 120.26472184895187 8:08
38.0% 38000 108615.67340018699 93614.86099439148 120.01833220121839 7:53
40.0% 40000 108189.84506806819 93740.44910207292 120.17934162929558 7:37
42.0% 42000 107939.65143493259 93834.05573605302 120.29934941410589 7:22
44.0% 44000 107817.74457492097 94328.42213398592 120.93314868428571 7:07
46.0% 46000 107645.46570848562 93837.58325264072 120.30387184412577 6:52
48.0% 48000 108009.67809897116 93728.81099086833 120.16442106345902 6:36
50.0% 50000 108136.9096828193 93568.82887522175 119.9593170180332 6:21
52.0% 52000 107895.9760849019 93795.82888007448 120.25034091854111 6:06
54.0% 54000 107886.52785448956 94409.24909202677 121.03677236751204 5:51
56.0% 56000 108117.69439003256 93699.94336879319 120.12741151370089 5:35
58.0% 58000 108090.49330945831 93610.20515342659 120.0123632101665 5:20
60.0% 60000 107685.93326105399 94339.93374000705 120.94790706496308 5:05
62.0% 62000 108413.78749640704 93574.87401785952 119.96706715438546 4:50
64.0% 64000 108240.0185759508 94063.73526929892 120.59380858680689 4:34
66.0% 66000 108650.74974339592 93739.3187409187 120.17789245702663 4:19
68.0% 68000 108106.430386827 93157.33329748194 119.43176175144505 4:04
70.0% 70000 108253.03044933111 94041.39083856613 120.56516205267707 3:48
72.0% 72000 107594.28937622189 93723.37105782727 120.15744683218914 3:33
74.0% 74000 108512.84246920235 93776.25733124021 120.22524933986725 3:18
76.0% 76000 107968.59366912623 93890.56791089648 120.37180048546071 3:03
78.0% 78000 108024.27839576988 93862.2149239762 120.33545072040258 2:47
80.0% 80000 107697.96213950164 93773.7042630089 120.22197619511371 2:32
82.0% 82000 108108.68537563692 94185.73338566914 120.75021548954557 2:17
84.0% 84000 108250.88620647619 93754.43863044186 120.19727681442845 2:02
86.0% 86000 108259.21231578765 93841.62919263162 120.30905891554517 1:46
88.0% 88000 108737.65904237027 93828.09609682369 120.29170889685012 1:31
90.0% 90000 108590.35111284332 94000.81867313881 120.51314676817357 1:16
92.0% 92000 108451.29548567365 93910.97405537797 120.39796205213196 1:01
94.0% 94000 108650.86069591212 94028.24587136641 120.54830963177517 0:45
96.0% 96000 108038.8927258845 94171.9106591098 120.73249415160302 0:30
98.0% 98000 108029.91821528843 93447.38024299666 119.80361457791906 0:15
100.0% 100000 108042.33377743122 94088.84286342005 120.62599761683903 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')