Single Chromosome Simulation

This tutorial should take between 20 to 30 minutes of reading and performing simulations.

Chromatin Dynamics Simulations on Chromosome 10 of GM12878 Cell Line

The first step is to import the OpenMiChroM module

from OpenMiChroM.ChromDynamics import MiChroM
from OpenMiChroM.CndbTools import cndbTools

MiChroM class sets the initial parameters of the simulation:

• time_step=0.01: set the simulation time step to perfom the integration

• temperature=1.0: set the temperature of your simulation

sim = MiChroM(temperature=1.0, time_step=0.01)

There are four hardware platform options to run the simulations:

platform="cuda"
platform="opencl"
platform="hip"
platform="cpu"

Choose accordingly.

sim.setup(platform="opencl")

Set the directory name in which the output of the simulation is saved:

sim.saveFolder('output_chr10')

The next step is to load the chromatin compartment sequence for chromosome 10 and generate an initial 3D structure to start the simulation. We can use the createSpringSpiral function to set the initial configuration of the polymer based in the sequence file.

The first column of the sequence file should contain the locus index. The second should have the locus type annotation. A template file of the chromatin sequence of types can be found here.

The loci positions are stored in the variable chr10 as a NumPy array \([N:3]\), where \(N\) is the number of beads.

chr10 = sim.createSpringSpiral(ChromSeq='inputs/chr10_beads.txt', isRing=False)

We can check the position of the first five beads:

print(chr10[:5])

The initial structure should then be loaded into the sim object.

The option center=True moves your system to the origin.

sim.loadStructure(chr10, center=True)

The initial 3D chromosome structure can be saved in .ndb file format. The file is stored in the path given in saveFolder.

sim.saveStructure(mode='ndb')

The next step is to add the force field in the simulation object sim.

In this tutorial, the forces can be divided into two sets:

MiChroM Homopolymer (Bonded) Potentials

sim.addFENEBonds(kfb=30.0)
sim.addAngles(ka=2.0)
sim.addRepulsiveSoftCore(Ecut=4.0)

MiChroM Non-Bonded Potentials

sim.addTypetoType(mu=3.22, rc=1.78)
sim.addIdealChromosome(mu=3.22, rc=1.78, dinit=3, dend=500)

The last potential adds a spherical constrain to collapse the initial structure.

sim.addFlatBottomHarmonic(kr=5*10**-3, n_rad=15.0)

Run a short simulation to generate a collapsed structure.

block = 3*10**2
n_blocks = 2*10**3

Two variables control the chromatin dynamics simulation steps:

block: The number of time-steps performed in each cycle (or block) n_blocks: The number of cycles (or blocks) simulated.

The initial collapse simulation will run for \(3\times10^2 \times 2\times10^3 = 6\times10^5\) time-steps.

for _ in range(n_blocks):
    sim.runSimBlock(block, increment=False)

Details about the output of each simulation block:

• bl=0: index number of the simulated block. The parameter increment=False is used to ignore the steps counting.

• pos[1]=[X,Y,Z]: spatial position for the locus 1.

• dr=1.26: average of the loci displacements in each block (in units of sigma).

• t=0: current simulation time.

• kin=1.5: kinetic energy of the system (reduced units).

• pot=19.90: total potential energy of the system (reduced units).

• RG=7.654: radius of gyration at the end of the simulated block.

• SPS=12312: steps per second of each block.

The radius of gyration is a good parameter to check the performance of the collapse. If the chromosome polymer is not collapsed, it is necessary to rerun the initial collapse steps. We can also save the structure for inspection.

print(sim.chromRG())
sim.saveStructure(mode='ndb')

The structure can also be saved using stardard file formats used for macromolecules, as the pdb and gro formats.

sim.saveStructure(mode='gro')
sim.saveStructure(mode='pdb')

The next step is to remove the spherical constrain force to run the production simulation.

sim.removeFlatBottomHarmonic()

If necessary, one could remove any of the forces applied in the system. To see the forces in the system:

sim.forceDict

# sim.removeForce(forceName="TypetoType")

To run the production simulation, it is necessary to initialize the .cndb file to save the chromatin dynamics trajectory.

sim.initStorage(filename="traj_chr10")

Set the parameters of the production simulation:

\(block = 5\times10^2\) \(n\_blocks = 2\times10^3\)

block = 5*10**2
n_blocks = 2*10**3

for _ in range(n_blocks):
    sim.runSimBlock(block, increment=True)
    sim.saveStructure()

Once the simulation is completed, it is necessary to close the .cndb file to avoid losing the trajectory data.

sim.storage[0].close()

The simulation should generate the traj_chr10_0.cndb trajectory file in the output_chr10 folder. This file contains 2000 frames (one snapshot per block).

Trajectory analysis using cndbTools

cndbTools is a class that allows analyses in the chromatin dynamics trajectories using the binary format .cndb (compact ndb).

cndbTools = cndbTools()

Load the cndb file in the variable chr10_traj.

chr10_traj = cndbTools.load('output_chr10/traj_chr10_0.cndb')

print(chr10_traj) # Print the information of the cndb trajectory.

Extract the loci XYZ position over the simulated 2000 frames and save in the variable chr10_xyz.

chr10_xyz = cndbTools.xyz(frames=[1,2001,1], XYZ=[0,1,2])

The variable chr10_xyz allows the cndbTools to perform several analyses. In this example, the radius of gyration can be obtained as a function of the simulated frames.

import matplotlib.pyplot as plt
import matplotlib as mpl


chr10_RG = cndbTools.compute_RG(chr10_xyz)
plt.plot(chr10_RG)
plt.ylabel(r'Radius of Gyration ($\sigma$)',fontsize=11)
plt.xlabel(r'Simulation Frames',fontsize=11)

cndbTools allows the selection of beads to compute the analyses. An example is the Radial Distribution Probability (RDP) for each chromatin subcompartments A1 and B1.

chr10_A1 = cndbTools.xyz(frames=[1,2001,1], beadSelection=chr10_traj.dictChromSeq['A1'], XYZ=[0,1,2])
chr10_B1 = cndbTools.xyz(frames=[1,2001,1], beadSelection=chr10_traj.dictChromSeq['B1'], XYZ=[0,1,2])

print("Computing RDP...")
r_A1, RDP_chr10_A1 = cndbTools.compute_RDP(chr10_A1, radius=15.0, bins=200)
r_B1, RDP_chr10_B1 = cndbTools.compute_RDP(chr10_B1, radius=15.0, bins=200)

plt.plot(r_A1, RDP_chr10_A1, color='red', label='A')
plt.plot(r_B1, RDP_chr10_B1, color='blue', label='B')
plt.xlabel(r'r ($\sigma$)', fontsize=11,fontweight='normal', color='k')
plt.ylabel(r'$\rho(r)/N_{type}$', fontsize=11,fontweight='normal', color='k')
plt.legend()
plt.gca().set_xlim([1/200,15.0])

We can also use cndbTools to generate the in silico Hi-C map (contact probability matrix).

In this tutorial, the trajectory contains 2,000 snapshots of chromosome 10 of the GM12878 cell line. For this set of structures, we expect the in silico Hi-C to not be fully converged due to inadequate sampling. To produce a converged map, it is recommended to simulate around 20 replicas with 10,000 frames on each, which generates an ensemble of 200,000 chromosome structures.

print("Generating the contact probability matrix...")
chr10_sim_HiC = cndbTools.traj2HiC(chr10_xyz)

plt.matshow(chr10_sim_HiC, norm=mpl.colors.LogNorm(vmin=0.001, vmax=chr10_sim_HiC.max()),cmap="Reds")
plt.colorbar()

To visualize the chromosome’s 3D structures in the standard visualization softwares for macromolecules, there are available scripts for converting the ndb/cndb file format to .pdb and .gro. For details, please check the Nucleome Data Bank.

The ndb plugin for visualizing the chromatin dynamics trajectories in VMD/Chimera/Pymol is under development.