Multiple Chromosome Simulation

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

Chromatin Dynamics Simulations on Chromosome 10 and Chromosome 11 of GM12878 Cell Line

The first step is to import the OpenMiChroM module

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

OpenMiChroM allows the simulation of multiples chromosomes. In this tutorial, the multiple chain simulation will be performed using the chromosomes 10 and 11 of the human GM12878 cell line. The system is generated based on the collapsed structure of each chromosome. This collapse step is similar to the one presented in the single chromosome simulation tutorial. For each chromosome, the collapse simulation should be performed individually.

Chromosome 10 Collapse Simulation

sim_chr10 = MiChroM(name="chr10", temperature=1.0, time_step=0.01)
sim_chr10.setup(platform="opencl")
sim_chr10.saveFolder('output_nucleus')
chr10 = sim_chr10.createSpringSpiral(ChromSeq='inputs/chr10_beads.txt')
sim_chr10.loadStructure(chr10, center=True)

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

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

print("Performing chr10 simulation...")
for _ in range(1000):
    sim_chr10.runSimBlock(500)

sim_chr10.saveStructure(filename="chr10", mode="ndb")
del sim_chr10

Chromosome 11 Collapse Simulation

sim_chr11 = MiChroM(name="chr11", temperature=1.0, time_step=0.01)
sim_chr11.setup(platform="opencl")
sim_chr11.saveFolder('output_nucleus')
chr11 = sim_chr11.createSpringSpiral(ChromSeq='inputs/chr11_beads.txt')
sim_chr11.loadStructure(chr11, center=True)

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

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

print("Performing chr11 simulation...")
for _ in range(1000):
    sim_chr11.runSimBlock(500)

sim_chr11.saveStructure(filename="chr11", mode="ndb")
del sim_chr11

Multiple Chromosome Simulation

Once the collapse simulations are done, the collapsed structures of each chromosome should be included in the same simulation system.

sim_nucleus = MiChroM(name="nucleus", temperature=1.0, time_step=0.01)

sim_nucleus.setup(platform="opencl")

sim_nucleus.saveFolder('output_nucleus')

The function initStructure() receives a list of files and saves the positions of the collapsed chromosomes in the variable initial_conf.

initial_conf = sim_nucleus.initStructure(
    CoordFiles=['output_nucleus/chr10_0_block1000.ndb',
                'output_nucleus/chr11_0_block1000.ndb']
)

initial_conf contains the coordinates of all beads for each chromosome. When used to load coordinate’s files (.ndb,.pdb, or .gro), the function initStructure prints the variable chains, which shows the information of each chromosome.

For each chain, we have (start,end,is_ring). start is the index of the first bead of the chain, end is the index for the last bead and is_ring indicates whether the first and last bead are connect or not. is_ring=0 represents a open chromosome chain and is_ring=1 represents a circular polymer (used for simulating bacteria genome, for example).

sim_nucleus.chains

Before loading the structures in the simulation context, it is necessary to spatially distribute the chromosome chains. This step is essential to guarantee no overlap between the chains in the initial configuration. It also helps randomizing the initial condition when simulating different replicas.

initial_conf = sim_nucleus.setFibPosition(initial_conf, factor=1.5)

This function distributes the center of mass of each chain in a spherical shell, according to the Fibonacci Sphere Algorithm. The argument factor sets the radius of the spherical shell, in comparison to the radius of the nucleus (see function documentation). For simulations with just a few chromosomes, factor=1.5 should be fine.

Now we can load the chromosomes with adjusted positions in the simulation context.

sim_nucleus.loadStructure(initial_conf, center=True)

The initial configuration of the chromosome structures can be saved in the .ndb file format. The file is stored in the path set in saveFolder. We advise saving the structure and double-checking if there is no overlap between the chains.

When having multiple chains in the simulation context, the saveStructure function will save each chromosome structure in a different file. Each chain is associated with an index starting from 0. This index follows the order of addition of the chains in the initStructure function. In this tutorial, chromosome 10 receives index 0, and chromosome 11 index 1.

sim_nucleus.saveStructure(mode='ndb')

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

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

MiChroM Homopolymer (Bonded) Potentials

sim_nucleus.addFENEBonds(kfb=30.0)
sim_nucleus.addAngles(ka=2.0)
sim_nucleus.addRepulsiveSoftCore(Ecut=4.0)
sim_nucleus.addFlatBottomHarmonic(n_rad=20)

MiChroM Non-Bonded Potentials

The addTypetoType interaction is independent of the chromosome chain and only depends on the chromatin subcompartment annotation for the interacting loci. This potential is added for both chains simultaneously.

sim_nucleus.addTypetoType(mu=3.22, rc=1.78)

On the other hand, for the Ideal Chromossome potential, we add the forces for each chromosome separately. The function addMultiChainIC receives the chromosome chain index and adds the IC potential.

The chromosome chain information can be obtained from the aforementioned variable chains. Note that the Ideal Chromosome potential was applied from the genomic distance \(d\) = 3 to \(d\) = 500. These cutoff values can be adjusted based on the chromosome length and user needs.

sim_nucleus.addMultiChainIC(chainIndex=0, mu=3.22, rc=1.78, dinit=3, dend=500)
sim_nucleus.addMultiChainIC(chainIndex=1, mu=3.22, rc=1.78, dinit=3, dend=500)

The simulation setup is complete!

Before running the production simulation, it is necessary to run a collapse and equilibration simulation as the chains are initially apart. The addFlatBottomHarmonic function used above includes in the simulation a harmonic potential to drive the chains together.

Set the parameters of the collapse simulation:

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

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

We can save the structure and check if the collapse was successful, with the two chains now interacting with each other.

sim_nucleus.saveStructure(mode='ndb')
sim_nucleus.saveStructure(mode='gro')

After the collapse, we should remove the harmonic potential and add the nucleus confinement.

sim_nucleus.removeFlatBottomHarmonic()
sim_nucleus.addAdditionalForce(sim_nucleus.addSphericalConfinementLJ)

The function addAdditionalForce adds a force to the system after the system has been already initialized. The initialization happens when the function runSimBlock is executed for the first time after the addition of the forces. In this case, we are using the function to add the Lennard-Jones spherical confinement (addSphericalConfinementLJ) in the system after removing the harmonic potential.

To run the production simulation, it is necessary to initialize the .cndb file to save the chromatin dynamics trajectory. The files will be saved in the output folder set using the saveFolder function. Each chromosome will be stored separately according to the chain index, as already mentioned.

sim_nucleus.initStorage('traj_nucleus')

Set the parameters of the production simulations:

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

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

Once the simulation is performed, it is necessary to close the .cndb files.

sim_nucleus.storage[0].close()
sim_nucleus.storage[1].close()

To visualize the chromosome’s 3D structures in the standard visualization softwares for macromolecules, there are available scripts for converting the ndb/cndb file formats 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.

sim_nucleus.saveStructure(mode="ndb")
sim_nucleus.saveStructure(mode="gro")

See the single chromosome simulation tutorial for examples on how to use our simulation analysis class `cndbTools <https://open-michrom.readthedocs.io/en/latest/OpenMiChroM.html#module-OpenMiChroM.CndbTools>`__