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

Note: This tutorial will be running in OpenMiChroM version 1.1.0 or greater. Please ensure you have the correct version installed before proceeding.

The first step is to import the OpenMiChroM module

from OpenMiChroM.ChromDynamics import MiChroM

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, timeStep=0.01)
sim_chr10.setup(platform="opencl")
sim_chr10.saveFolder('output_nucleus')
chr10 = sim_chr10.initStructure(mode='random',ChromSeq='Tutorials/Chromosomes_simulations/inputs/chr10_beads.txt', isRing=False)
sim_chr10.loadStructure(chr10, center=True)

sim_chr10.addHomoPolymerForces()
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)

sim_chr10.createSimulation()

print("Performing chr10 simulation...")
sim_chr10.run(nsteps=10**6, report=True, interval=10**5)

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

Chromosome 11 Collapse Simulation

sim_chr11 = MiChroM(name="chr11", temperature=1.0, timeStep=0.01)
sim_chr11.setup(platform="opencl")
sim_chr11.saveFolder('output_nucleus')
chr11 = sim_chr11.initStructure(mode='random',ChromSeq='Tutorials/Chromosomes_simulations/inputs/chr11_beads.txt', isRing=False)
sim_chr11.loadStructure(chr11, center=True)

sim_chr11.addHomoPolymerForces()
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)

sim_chr11.createSimulation()

print("Performing chr11 simulation...")
sim_chr11.run(nsteps=10**6, report=True, interval=10**5)

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, timeStep=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.ndb',
                'output_nucleus/chr11.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.

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()
sim_nucleus.addAngles()
sim_nucleus.addSelfAvoidance()
sim_nucleus.addFlatBottomHarmonic(nRad=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)

sim_nucleus.createSimulation()

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:

collapse_steps = 10**6
production_steps = 10**6

sim_nucleus.run(nsteps=collapse_steps, report=True, interval=10**5, totalSteps=(collapse_steps+production_steps))

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.

Now we create the reporters to save the simulation infos. There are 3 types of reporters:

statistics: Attaches a reporter to collect simulation statistics such as step number, radius of gyration (RG), total energy, potential energy, kinetic energy, and temperature.

trajectory: Attaches a reporter to save trajectory data (xyz per bead per chain) during the simulation. The file format to save the trajectory data. Options are ‘cndb’, ‘swb’,‘ndb’, ‘pdb’, ‘gro’, ‘xyz’. (Default: ‘cndb’)

energy components: Saves energy components per force group to a separate file named ‘energyComponents.txt’ in the simulation folder. Requires that statistics is True

set the number of steps interval we will save this information, here I chose 1000 steps

sim_nucleus.createReporters(statistics=True, traj=True, trajFormat="cndb", energyComponents=True, interval=1000)

sim_nucleus.run(nsteps=production_steps, report=True, interval=10**5, totalSteps=(collapse_steps+production_steps))

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")
sim_nucleus.saveStructure(mode="pdb")