OpenMiChroM

OpenMiChroM.ChromDynamics

The ChromDynamics classes perform chromatin dynamics based on the compartment annotation sequences of chromosomes. Simulations can be performed either using the default parameters of MiChroM (Minimal Chromatin Model) or using custom values for the type-to-type and Ideal Chromosome parameters.

Details about the MiChroM energy function and default parameters are described in: Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G., & Onuchic, J.N. (2016). Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), 12168-12173.

class OpenMiChroM.ChromDynamics.MiChroM(name='OpenMiChroM', timeStep=0.01, collisionRate=0.1, temperature=1.0)

Bases: object

The MiChroM class performs chromatin dynamics simulations using the default MiChroM energy function parameters for type-to-type and Ideal Chromosome interactions.

Details about the MiChroM (Minimal Chromatin Model) energy function and the default parameters are described in: Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G., & Onuchic, J.N. (2016). Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), 12168-12173.

The MiChroM class sets up the environment to start chromatin dynamics simulations.

Parameters

• name (str, optional) – Name used in the output files. Defaults to “OpenMichrom”.

• timeStep (float, optional) – Simulation time step in units of τ. Defaults to 0.01.

• collisionRate (float, optional) – Friction/damping constant in units of reciprocal time (1/τ). Defaults to 0.1.

• temperature (float, optional) – Temperature in reduced units. Defaults to 1.0.

addAdditionalForce(forceFunction, *args, **kwargs)

Add an additional force after the system has already been initialized.

Parameters

• forceFunciton (function, required) – Force function to be added. Example: addSphericalConfinementLJ

• **args (collection of arguments, required) – Arguments of the function to add the force. Consult respective documentation.

addAngles(kA=2.0)

Adds an angular potential between bonds connecting beads i-1, i, and i+1.

This function adds angle forces to the system to enforce stiffness between sequential beads, according to the method described in: Halverson, J.D., Lee, W.B., Grest, G.S., Grosberg, A.Y., & Kremer, K. (2011). Molecular dynamics simulation study of nonconcatenated ring polymers in a melt. I. Statics. The Journal of Chemical Physics, 134(20), 204904.

Parameters

kA (float or array-like, optional) – Angle potential coefficient(s). If a single float is provided, the same coefficient is used for all angles. If an array is provided, it must have a length of N-2 (number of angles), and each angle will use the corresponding coefficient. Defaults to 2.0.

Raises

ValueError – If kA is an array and its length does not match the expected number of angles.

addCorrelatedNoise(act_seq='none', atol=1e-05)

This function assigns active force to monomers, as defined by the act_seq array. :param act_seq: This list should contain the same number of elements as the number of monomers, where the value at index i corresponds to the active force F for the bead i :type act_seq: list, required :param atol: tolerance for force values. If force is less than atol then it is ignored. :type atol: float, optional

Returns Error if act_seq is not the same length as the sequence file.

addCustomIC(mu=3.22, rc=1.78, dinit=3, dend=200, IClist=None, CutoffDistance=3.0)

Adds the Ideal Chromosome potential using custom values for interactions between beads separated by a genomic distance \(d\). The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 3.22).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 1.78).

• dinit (int, required) – The first neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential. (Default value = 3).

• dend (int, required) – The last neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential. (Default value = 200).

• IClist (file, optional) – A one-column text file containing the energy interaction values for loci i and j separated by a genomic distance \(d\). (Default value: None).

addCustomMultiChainIC(mu=3.22, rc=1.78, dinit=3, dend=1000, chainIndex=None, IClist=None, CutoffDistance=3.0)

addCustomTypes(name='CustomTypes', mu=3.22, rc=1.78, TypesTable=None, CutoffDistance=3.0)

Adds the type-to-type potential using custom values for interactions between the chromatin types. The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

The function receives a txt/TSV/CSV file containing the upper triangular matrix of the type-to-type interactions. A file example can be found here.

Parameters

• name (string, required) – Name to customType Potential. (Default value = “CustomTypes”)

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 3.22).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 1.78).

• TypesTable (file, required) – A txt/TSV/CSV file containing the upper triangular matrix of the type-to-type interactions. (Default value: None).

addCylindricalConfinement(rConf=5.0, zConf=10.0, kConf=30.0)

Adds a cylindrical confinement potential to the system.

This potential confines particles within a cylinder of radius rConf and height 2 * zConf. Particles outside this cylinder experience a harmonic restoring force pushing them back inside.

Parameters

• rConf (float, optional) – Radius of the cylindrical confinement in nanometers. Defaults to 5.0.

• zConf (float, optional) – Half-height of the cylindrical confinement along the z-axis in nanometers. Defaults to 10.0.

• kConf (float, optional) – Force constant (stiffness) of the confinement potential. Defaults to 30.0.

addFENEBonds(kFb=30.0, bonds=None)

Adds FENE (Finite Extensible Nonlinear Elastic) bonds to the system.

This function initializes the FENE bond force if it has not been added yet, and adds bonds between specified pairs of loci. By default, if no bonds are specified, it will add bonds between neighboring loci in each chain according to the chain definitions. If a chain is specified as a ring, it will also add a bond between the first and last loci of that chain.

The FENE bonds are defined according to the method described in: Halverson, J.D., Lee, W.B., Grest, G.S., Grosberg, A.Y., & Kremer, K. (2011). Molecular dynamics simulation study of nonconcatenated ring polymers in a melt. I. Statics. The Journal of Chemical Physics, 134(20), 204904.

Parameters

• kFb (float, optional) – Bond coefficient (force constant) for the FENE bonds. Defaults to 30.0.

• bonds (list of tuples, optional) – A list of tuples specifying the pairs of loci to bond. Each tuple is (i, j). If None, bonds between neighboring loci in the chains will be added. Defaults to None.

Raises

ValueError – If any of the loci indices are out of bounds.

addFlatBottomHarmonic(kR=0.005, nRad=10.0)

Adds a flat-bottom harmonic potential to confine the chromosome chain inside the nucleus wall.

The potential is defined as:

V(r) = step(r - r0) * (kR / 2) * (r - r0)^2

where:

• r is the distance from the origin (center of the nucleus)

• r0 (nRad) is the nucleus radius

• kR is the spring constant of the potential

This potential applies no force when particles are inside the nucleus (r ≤ r0) and applies a harmonic restoring force when particles are outside the nucleus (r > r0).

Parameters

• kR (float, optional) – Spring constant of the harmonic potential. Defaults to 5e-3.

• nRad (float, optional) – Nucleus radius in units of σ. Defaults to 10.0.

addHarmonicBonds(bondStiffness=30.0, equilibriumDistance=1.0)

Adds harmonic bonds to all monomers within each chain. For each chain, bonds are created between consecutive monomers. If the chain forms a ring, an additional bond is created between the first and last monomer to complete the ring.

Parameters

• bondStiffness (float, optional) – The stiffness of the bond. Default is 30.0.

• equilibriumDistance (float, optional) – The equilibrium distance for the bond. Default is 1.0.

addHomoPolymerForces(harmonic=False, softcoreLJ=False, kFb=30, kA=2.0, eCut=4.0)

Adds homopolymer forces to the system based on the specified parameters.

Depending on the harmonic flag, this function will add either harmonic bonds or FENE bonds to the polymer chains. Additionally, it sets up angle forces and repulsive soft-core interactions.

Parameters

• harmonic (bool, optional) – If True, adds harmonic bonds to the monomers. If False, adds FENE bonds. Default is False.

• KFb (float, optional) – The stiffness of the bonds. Relevant when harmonic is True. Default is 30.0.

• kA (float, optional) – The stiffness of the angles between bonded monomers. Default is 2.0.

• eCut (float, optional) – The cutoff distance for the repulsive soft-core potential. Default is 4.0.

addIdealChromosome(mu=3.22, rc=1.78, Gamma1=- 0.03, Gamma2=- 0.351, Gamma3=- 3.727, dinit=3, dend=500, CutoffDistance=3.0)

Adds the Ideal Chromosome potential for interactions between beads separated by a genomic distance \(d\) according to the MiChroM energy function parameters reported in “Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G. and Onuchic, J.N., 2016. Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), pp.12168-12173”.

The set of parameters \(\{\gamma_d\}\) of the Ideal Chromosome potential is fitted in a function: \(\gamma(d) = \frac{\gamma_1}{\log{(d)}} +\frac{\gamma_2}{d} +\frac{\gamma_3}{d^2}\).

The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 3.22).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 1.78).

• Gamma1 (float, required) – Ideal Chromosome parameter. (Default value = -0.030).

• Gamma2 (float, required) – Ideal Chromosome parameter. (Default value = -0.351).

• Gamma3 (float, required) – Ideal Chromosome parameter. (Default value = -3.727).

• dinit (int, required) – The first neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential. (Default value = 3).

• dend (int, required) – The last neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential. (Default value = 500).

addLoops(mu=3.22, rc=1.78, X=- 1.61299, looplists=None)

Adds the Loops interactions according to the MiChroM energy function parameters reported in “Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G. and Onuchic, J.N., 2016. Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), pp.12168-12173”.

The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Note

For Multi-chain simulations, the ordering of the loop list files is important! The order of the files should be the same as used in the other functions.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 3.22).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 1.78).

• X (float, required) – Loop interaction parameter. (Default value = -1.612990).

• looplists (list[str], required) – List with the names of the files containing loop information. Each file should be a two-column text file containing the index i and j of the loci pairs that forms the loop anchors. (Default value: None).

addMultiChainIC(mu=3.22, rc=1.78, Gamma1=- 0.03, Gamma2=- 0.351, Gamma3=- 3.727, dinit=3, dend=500, chainIndex=0, CutoffDistance=3.0)

Adds the Ideal Chromosome potential for multiple chromosome simulations. The interactions between beads separated by a genomic distance \(d\) is applied according to the MiChroM energy function parameters reported in “Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G. and Onuchic, J.N., 2016. Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), pp.12168-12173”.

The set of parameters \(\{\gamma_d\}\) of the Ideal Chromosome potential is fitted in a function: \(\gamma(d) = \frac{\gamma_1}{\log{(d)}} +\frac{\gamma_2}{d} +\frac{\gamma_3}{d^2}\).

The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 3.22).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 1.78).

• Gamma1 (float, required) – Ideal Chromosome parameter. (Default value = -0.030).

• Gamma2 (float, required) – Ideal Chromosome parameter. (Default value = -0.351).

• Gamma3 (float, required) – Ideal Chromosome parameter. (Default value = -3.727).

• dinit (int, required) – The first neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential. (Default value = 3).

• dend (int, required) – The last neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential. (Default value = 500).

• chainIndex (integer, required) – The index of the chain to add the Ideal Chromosome potential. All chains are stored in self.chains. (Default value: 0).

addRepulsiveSoftCore(eCut=4.0, cutoffDistance=3.0)

Adds a soft-core repulsive interaction that allows chain crossing, representing the activity of topoisomerase II.

Details can be found in the following publications:

• Oliveira Jr., A.B., Contessoto, V.G., Mello, M.F., & Onuchic, J.N. (2021). A scalable computational approach for simulating complexes of multiple chromosomes. Journal of Molecular Biology, 433(6), 166700.

• Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G., & Onuchic, J.N. (2016). Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), 12168-12173.

• Naumova, N., Imakaev, M., Fudenberg, G., Zhan, Y., Lajoie, B.R., Mirny, L.A., & Dekker, J. (2013). Organization of the mitotic chromosome. Science, 342(6161), 948-953.

Parameters

• eCut (float, optional) – Energy cost for chain crossing in units of (k_B T). Defaults to 4.0.

• cutoffDistance (float, optional) – Cutoff distance for the nonbonded interactions. Defaults to 3.0.

addSelfAvoidance(Ecut=4.0, k_rep=20.0, r0=1.0)

This adds Soft-core self avoidance between all non-bonded monomers. This force is well behaved across all distances (no diverging parts) :param Ecut: energy associated with full overlap between monomers :type Ecut: float, required :param k_rep: steepness of the sigmoid repulsive potential :type k_rep: float, required :param r0: distance from the center at which the sigmoid is half as strong :type r0: float, required

addSphericalConfinementLJ(radius='density', density=0.1)

Adds a spherical confinement potential to the system according to the MiChroM energy function.

This potential describes the interaction between the chromosome and a spherical wall, effectively confining the particles within a sphere of specified radius.

Parameters

• radius (float or str, optional) – Radius of the spherical confinement. If set to “density”, the radius is calculated based on the specified density. Defaults to “density”.

• density (float, optional) – Density of the chromosome beads inside the nucleus. Required if radius is “density”. Defaults to 0.1.

Notes

• If radius is “density”, the radius is calculated using the formula:

radius = (3 * N / (4 * π * density)) ** (1/3)

where N is the number of particles in the system. - The confinement potential is modeled using a shifted Lennard-Jones potential.

addTypetoType(mu=3.22, rc=1.78)

Adds the type-to-type interactions according to the MiChroM energy function parameters reported in “Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G. and Onuchic, J.N., 2016. Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), pp.12168-12173”.

The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 3.22).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 1.78).

buildClassicMichrom(ChromSeq=None, chromosome=None, CoordFiles=None, mode='auto')

Builds a Classic Michrom simulation setup by initializing the structure, loading it, adding polymer forces, defining interactions between types, setting up the ideal chromosome, adding flat-bottom harmonic potentials, and creating the simulation.

Parameters

• ChromSeq (str, optional) – Path to the chromosome sequence file. If None, a default path is used. Default is None.

• CoordFiles (str, optional) – Path to the coordinate files. If None, a default path is used. Default is None.

• mode (str, optional) – Mode of initialization for the structure. Supported modes include ‘spring’, ‘line’, etc. Default is ‘auto’.

Example

```python simulation = MichromSimulation() simulation.buildClassicMichrom(

ChromSeq=”/path/to/chromosome_sequence.txt”, CoordFiles=”/path/to/coordinate_files/”, mode=’spring’

createLine(ChromSeq, chromosome=None)

Creates a straight line for the initial configuration of the chromosome polymer.

Parameters

• ChromSeq (file, required) – Chromatin sequence of types file. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB).

• length_scale (float, required) – Length scale used in the distances of the system in units of reduced length \(\sigma\). (Default value = 1.0).

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

createRandomWalk(ChromSeq=None, chromosome=None)

Creates a chromosome polymer chain with beads position based on a random walk.

Args:

ChromSeq (file, required):

Chromatin sequence of types file. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB).

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

createReporters(statistics=True, traj=False, trajFormat='cndb', outputName=None, energyComponents=False, interval=1000)

Configures and attaches reporters to the simulation for data collection during simulation runs. This method sets up custom reporters for the OpenMM Simulation object to collect simulation statistics and/or save trajectory data at specified intervals. It supports saving energies per force group and various trajectory file formats.

Parameters

• statistics (bool, optional) – If True, attaches a reporter to collect simulation statistics such as step number, radius of gyration (RG), total energy, potential energy, kinetic energy, and temperature. (Default: True)

• traj (bool, optional) – If True, attaches a reporter to save trajectory data during the simulation. (Default: False)

• file_name (str, optional) – The file name for saving trajectory data. If None, defaults to self.name. (Default: None)

• trajFormat (str, optional) – The file format to save the trajectory data. Options are ‘cndb’, ‘swb’,`’ndb’, `’pdb’, ‘gro’, ‘xyz’. (Default: ‘cndb’)

• energyComponents (bool, optional) – If True, saves energy components per force group to a separate file named ‘energyComponents.txt’ in the simulation folder. Requires that statistics is True. (Default: False)

• interval (int, optional) – The interval (in time steps) at which to report data. (Default: 1000)

createSimulation()

Initializes the simulation context and adds forces to the system.

This function checks if the simulation context has already been created. If not, it loads the particles, processes any exceptions (bonds that should not be included in nonbonded interactions), adds forces to the system, and sets up the simulation context.

createSpringSpiral(ChromSeq=None, isRing=False, chromosome=None)

Creates a spring-spiral-like shape for the initial configuration of the chromosome polymer.

Parameters

• ChromSeq (file, required) – Chromatin sequence of types file. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB).

• isRing (bool, optional) – Whether the chromosome chain is circular or not (Used to simulate bacteria genome, for example). f bool(isRing) is True , the first and last particles of the chain are linked, forming a ring. (Default value = False).

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

getForces()

Gets the energy values for each force applied in the system.

getLoops(looplists)

Get the loop position (CTFC anchor points) for each chromosome.

Note

For Multi-chain simulations, the ordering of the loop list files is important! The order of the files should be the same as used in the other functions.

Parameters

looplists (list[str]) – List with the names of the files containing loop information. Each file should be a two-column text file containing the index i and j of the loci pairs that forms the loop anchors.

getPositions()

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

initPositions()

Internal function that sets the locus coordinates in the OpenMM system.

Raises

ValueError – If the simulation context has not been initialized.

initStructure(mode='auto', CoordFiles=None, ChromSeq=None, isRing=False, chromosome=None)

Creates the coordinates for the initial configuration of the chromosomal chains and sets their sequence information.

Args:

mode (str, required):

• ‘auto’ - Creates a spring-spiral-like shape when a CoordFiles is not provided. If CoordFiles is provided, it loads the respective type of coordinate files (.ndb, .gro, or .pdb). (Default value = ‘auto’).

• ‘line’ - Creates a straight line for the initial configuration of the chromosome polymer. Can only be used to create single chains.

• ‘spring’ - Creates a spring-spiral-like shape for the initial configuration of the chromosome polymer. Can only be used to create single chains.

• ‘random’ - Creates a chromosome polymeric chain with beads positions based on a random walk. Can only be used to create single chains.

• ‘ndb’ - Loads a single or multiple .ndb files and gets the position and types of the chromosome beads.

• ‘pdb’ - Loads a single or multiple .pdb files and gets the position and types of the chromosome beads.

• ‘gro’ - Loads a single or multiple .gro files and gets the position and types of the chromosome beads.

CoordFiles (list of files, optional):

List of files with xyz information for each chromosomal chain. Accepts .ndb, .pdb, and .gro files. All files provided in the list must be in the same file format.

ChromSeq (list of files, optional):

List of files with sequence information for each chromosomal chain. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB). If the chromatin types considered are different from the ones used in the original MiChroM (A1, A2, B1, B2, B3, B4, and NA), the sequence file must be provided when loading .pdb or .gro files, otherwise, all the chains will be defined with ‘NA’ type. For the .ndb files, the sequence used is the one provided in the file.

isRing (bool, optional):

Whether the chromosome chain is circular or not (used to simulate bacteria genome, for example). To be used with the option 'random'. If bool(isRing) is True , the first and last particles of the chain are linked, forming a ring. (Default value = False).

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

initVelocities()

Internal function that sets the initial velocities of the loci in the OpenMM system.

Raises

ValueError – If the simulation context has not been initialized.

loadBed(bedFile, chromosome)

Reads a BED file and extracts types for a specified chromosome.

Parameters

• bedFile (str) – Path to the BED file.

• chromosome (str) – The chromosome identifier (e.g., ‘chr1’).

Returns

A list of tuples containing start positions and types.

Return type

List[Tuple[int, str]]

loadGRO(GROfiles=None, ChromSeq=None)

Loads a single or multiple .gro files and gets position and types of the chromosome beads. Initially, the MiChroM energy function was implemented in GROMACS. Details on how to run and use these files can be found at the Nucleome Data Bank.

• Contessoto, V.G., Cheng, R.R., Hajitaheri, A., Dodero-Rojas, E., Mello, M.F., Lieberman-Aiden, E., Wolynes, P.G., Di Pierro, M. and Onuchic, J.N., 2021. The Nucleome Data Bank: web-based resources to simulate and analyze the three-dimensional genome. Nucleic Acids Research, 49(D1), pp.D172-D182.

Parameters

• GROfiles (list of files, required) – List with a single or multiple files in .gro file format. (Default value: None).

• ChromSeq (list of files, optional) – List of files with sequence information for each chromosomal chain. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB). If the chromatin types considered are different from the ones used in the original MiChroM (A1, A2, B1, B2, B3, B4, and NA), the sequence file must be provided, otherwise all the chains will be defined with ‘NA’ type.

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

loadNDB(NDBfiles=None)

Loads a single or multiple .ndb files and gets position and types of the chromosome beads. Details about the NDB file format can be found at the Nucleome Data Bank.

• Contessoto, V.G., Cheng, R.R., Hajitaheri, A., Dodero-Rojas, E., Mello, M.F., Lieberman-Aiden, E., Wolynes, P.G., Di Pierro, M. and Onuchic, J.N., 2021. The Nucleome Data Bank: web-based resources to simulate and analyze the three-dimensional genome. Nucleic Acids Research, 49(D1), pp.D172-D182.

Parameters

NDBfiles (file, required) – Single or multiple files in .ndb file format. (Default value: None).

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

loadPDB(PDBfiles=None, ChromSeq=None)

Loads a single or multiple .pdb files and gets position and types of the chromosome beads. Here we consider the chromosome beads as the carbon-alpha to mimic a protein. This trick helps to use the standard macromolecules visualization software. The type-to-residue conversion follows: {‘ALA’:0, ‘ARG’:1, ‘ASP’:2, ‘GLU’:3,’GLY’:4, ‘LEU’ :5, ‘ASN’ :6}.

Parameters

• PDBfiles (list of files, required) – List with a single or multiple files in .pdb file format. (Default value: None).

• ChromSeq (list of files, optional) – List of files with sequence information for each chromosomal chain. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB). If the chromatin types considered are different from the ones used in the original MiChroM (A1, A2, B1, B2, B3, B4, and NA), the sequence file must be provided, otherwise all the chains will be defined with ‘NA’ type.

Returns

Returns an array of positions.

Return type

\((N, 3)\) numpy.ndarray

loadStructure(data, center=True, massList=None)

Loads the 3D positions of each bead of the chromosome polymer into the OpenMM system.

Parameters

• data (array-like) – The initial positions of the beads. Should be an array of shape (N, 3) or (3, N).

• center (bool, optional) – If True, centers the chromosome’s center of mass at [0, 0, 0]. Defaults to True.

• massList (array-like, optional) – Masses of each chromosome bead in units of μ. If None, all masses are set to 1.0. Defaults to None.

Raises

ValueError – If the input data is not in the correct format or contains NaN values.

printForces()

Prints the energy values for each force applied in the system.

printHeader()

random_ChromSeq(Nbeads)

Creates a random sequence of chromatin types for the chromosome beads.

Parameters

Nbeads (int, required) – Number of beads of the chromosome polymer chain. (Default value = 1000).

Returns

Returns an 1D array of a randomized chromatin type annotation sequence.

Return type

\((N, 1)\) numpy.ndarray

removeFlatBottomHarmonic()

Remove FlatBottomHarmonic force from the system.

removeForce(forceName)

Remove force from the system.

run(nsteps, report=True, interval=10000, totalSteps=None, checkSystem=False, blockSize=1000)

Executes the simulation for a specified number of steps, with optional progress reporting.

This function runs the simulation by advancing it through the defined number of steps (nsteps). If report is set to True, it adds a StateDataReporter to the simulation’s reporters to output progress information to the standard output (stdout) at every interval steps.

Parameters

• nsteps – int Number of steps to run the simulation. Must be a positive integer.

• report – bool, default=True If True, progress will be reported periodically via the StateDataReporter.

• interval – int, default=10000 Number of steps between progress reports. Must be a positive integer.

• totalSteps – int, optional The total number of steps for the entire simulation (used for progress percentage). If None, defaults to nsteps.

• checkSystem – bool, default=False If True, run the simulation in blocks of size blockSize. After each block, the system is checked for NaN energies; if found, it attempts to reset the system with a new velocities assignment.

• blockSize – int, default=1000 Size of each simulation block when checkSystem is True. Must be a positive integer.

Example

`python # Run the simulation for 500,000 steps with progress reporting every 5,000 steps simulation.run(nsteps=500000, report=True, interval=5000) `

saveFolder(folderPath)

Sets the folder path to save data.

Parameters

folderPath (str) – The folder path where simulation data will be saved. If the folder does not exist, it will be created.

saveStructure(fileName=None, mode='gro')

Saves the 3D positions of beads during the simulation in various file formats.

This method exports the simulation’s bead positions to files in formats such as XYZ, PDB, GRO, and NDB. It supports multiple chains and assigns residue names based on bead types. The method ensures that the output directory exists and handles different file formats appropriately.

Parameters

• fileName (str, optional) – The name of the file to save the structure. If None, the fileName is automatically generated using the simulation’s name and current step number with the specified mode as the file extension. (Default: None)

• mode (str, optional) – The file format to save the structure. Supported formats are: - ‘xyz’: XYZ format. - ‘pdb’: Protein Data Bank format. - ‘gro’: GROMACS GRO format. - ‘ndb’: Nucleome Data Bank format. (Default: ‘gro’)

setChains(chains=None)

Sets the configuration of the chains in the system.

This information is used later for adding bonds and angles in the homopolymer potential.

Parameters

chains (list of tuples, optional) – A list of chains in the format [(start, end, isRing)]. isRing is a boolean indicating whether the chromosome chain is circular or not (used to simulate bacterial genomes, for example). The particle range should be semi-open; for example, a chain (0, 3, False) links the particles 0, 1, and 2. If isRing is True, the first and last particles of the chain are linked, forming a ring. Defaults to [(0, None, False)], which links all particles of the system into one chain.

setFibPosition(positions, returnCM=False, factor=1.0)

Distributes the center of mass of chromosomes on the surface of a sphere according to the Fibonacci Sphere algorithm.

Parameters

• positions (\((Nbeads, 3)\) numpy.ndarray, required) – The array of positions of the chromosome chains to be distributed in the sphere surface.

• returnCM (bool, optional) – Whether to return an array with the center of mass of the chromosomes. (Default value: False).

• factor (float, optional) –

  Scale coefficient to be multiplied to the radius of the nucleus, determining the radius of the sphere in which the center of mass of chromosomes will be distributed. The radius of the nucleus is calculated based on the number of beads to generate a volume density of 0.1.

  \(R_{sphere} = factor * R_{nucleus}\)

Returns

Returns an array of positions to be loaded into OpenMM using the function loadStructure.

\((Nchains, 3)\) numpy.ndarray:

Returns an array with the new coordinates of the center of mass of each chain.

Return type

\((Nbeads, 3)\) numpy.ndarray

setPositions(beadsPos, randomize=False, randomOffset=1e-05)

Sets the 3D positions of each bead of the chromosome polymer in the OpenMM system.

Parameters

• beadsPos (numpy.ndarray of shape (N, 3)) – Array of XYZ positions for each bead (locus) in the polymer model.

• randomize (bool, optional) – If True, adds a small random offset to the positions to avoid numerical instability. Defaults to False.

• randomOffset (float, optional) – The magnitude of the random offset to be added if randomize is True. Defaults to 1e-5.

setup(platform='CUDA', gpu='default', integrator='langevin', precision='mixed', deviceIndex='0')

Sets up the simulation environment.

Tries to select the computational platform in the following priority order: the specified platform, then ‘CUDA’, ‘HIP’, ‘OpenCL’, and ‘CPU’. If the preferred platform is not available, it will attempt to use the next available platform in the list and print an informational message.

Parameters

• platform (str, optional) – The preferred computation platform to use. Defaults to ‘CUDA’.

• gpu (str, optional) – GPU device index or ‘default’. Defaults to ‘default’.

• integrator (str or OpenMM Integrator, optional) – The integrator to use for the simulation. Can be a string (e.g., ‘langevin’) or an OpenMM Integrator object. Defaults to ‘langevin’.

• precision (str, optional) – The floating point precision to use (‘mixed’, ‘single’, or ‘double’). Defaults to ‘mixed’.

• deviceIndex (str, optional) – The device index to use if specifying a GPU device. Defaults to ‘0’.

Raises

• ValueError – If an unknown integrator or precision is specified.

• Exception – If no suitable computational platform is available.

OpenMiChroM.Optimization

The Optimization classes perform the energy function parameters training of the chromosomes based on experimental Hi-C data.

class OpenMiChroM.Optimization.AdamTraining(mu=2.0, rc=2.0, eta=0.01, beta1=0.9, beta2=0.999, epsilon=1e-08, it=1, updateNeeded=False, update_storagePath='', method='classic')

Bases: object

The AdamTraining class performs the parameters training for each selected loci pair interaction.

Details about the methodology are decribed in “Zhang, Bin, and Peter G. Wolynes. “Topology, structures, and energy landscapes of human chromosomes.” Proceedings of the National Academy of Sciences 112.19 (2015): 6062-6067.”

The AdamTraining class receive a Hi-C matrix (text file) as input. The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 2.0).

• rc (float, required) – Parameter in the probability of crosslink function. (Default value = 2.0).

• eta (float, required) – Learning rate applied in each step (Default value = 0.01).

• beta1 (float, required) – The hyper-parameter of Adam are initial decay rates used when estimating the first and second moments of the gradient. (Default value = 0.9).

• beta2 (float, required) – The hyper-parameter of Adam are initial decay rates used when estimating the first and second moments of the gradient. (Default value = 0.999).

• it (int, required) – The iteration step

• method (str, required) – ‘classic’: Adam ‘qh’: quassi-hyperbolic momentum Adam

getHiCexp(HiC, centerRemove=False, centrange=[0, 0], norm=True, cutoff_low=0.0, cutoff_high=1.0, KR=False, neighbors=0)

Receives the experimental Hi-C map (Full dense matrix) in a text format and performs the data normalization from Hi-C frequency/counts/reads to probability.

Parameters

• HiC (file, required) – Experimental Hi-C map (Full dense matrix) in a text format.

• centerRemove (bool, optional) – Whether to set the contact probability of the centromeric region to zero. (Default value: False).

• centrange (list, required if centerRemove = True)) – Range of the centromeric region, i.e., centrange=[i,j], where i and j are the initial and final beads in the centromere. (Default value = [0,0]).

• cutoff (float, optional) – Cutoff value for reducing the noise in the original data. Values lower than the cutoff are considered \(0.0\).

getLamb(Lambdas, fixedPoints=None)

Calculates the Lagrange multipliers of each pair of interaction and returns the matrix containing the energy values for the optimization step.

Parameters

• Lambdas (file, required) – The matrix containing the energies values used to make the simulation in that step.

• fixedPoints (list, optional) – List of all pairs (i,j) of interactions that will remain unchanged throughout the optimization procedure.

Returns

Returns an updated matrix of interactions between each pair of bead.

Return type

\((N,N)\) numpy.ndarray

getPars(HiC, centerRemove=False, centrange=[0, 0], cutoff='deprecate', norm=True, cutoff_low=0.0, cutoff_high=1.0, KR=False, neighbors=0)

Receives the experimental Hi-C map (Full dense matrix) in a text format and performs the data normalization from Hi-C frequency/counts/reads to probability.

Parameters

• HiC (file, required) – Experimental Hi-C map (Full dense matrix) in a text format.

• centerRemove (bool, optional) – Whether to set the contact probability of the centromeric region to zero. (Default value: False).

• centrange (list, required if centerRemove = True)) – Range of the centromeric region, i.e., centrange=[i,j], where i and j are the initial and final beads in the centromere. (Default value = [0,0]).

• cutoff (float, optional) – Cutoff value for reducing the noise in the original data. Values lower than the cutoff are considered \(0.0\).

knight_ruiz_balance(matrix, tol=1e-05, max_iter=100)

Perform the Knight-Ruiz matrix balancing.

normalize_matrix(matrix)

Normalize the matrix for simulation optimization. Here the first neighbor should have the probability of contact P=1.0.

probCalc(state)

Calculates the contact probability matrix for a given state.

reset_Pi()

Resets Pi matrix to zeros

class OpenMiChroM.Optimization.CustomMiChroMTraining(ChromSeq='chr_beads.txt', TypesTable=None, mu=3.22, rc=1.78, cutoff=0.0, IClist=None, dinit=3, dend=200)

Bases: object

The CustomMiChroMTraining class performs the parameters training employing MiChroM (Minimal Chromatin Model) energy function.

Details about the methodology are decribed in “Di Pierro, Michele, et al. “Transferable model for chromosome architecture.” Proceedings of the National Academy of Sciences 113.43 (2016): 12168-12173.”

The CustomMiChroMTraining class receive a Hi-C matrix (text file) as input. The parameters \(\mu\) (mi) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

CustomMiChroMTraining optimizes the type-to-type (Types) and the Ideal Chromosome (IC) potential parameters separately.

Parameters

• ChromSeq (file, required) – Chromatin sequence of types file. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB).

• TypesTable (file, required) – A txt/TSV/CSV file containing the upper triangular matrix of the type-to-type interactions. (Default value: None).

• mu (float, required) – Parameter in the probability of crosslink function (Default value = 3.22, for human chromosomes in interphase).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\) (Default value = 1.78, for human chromosomes in interphase).

• cutoff (float, optional) – Cutoff value for reducing the noise in the original data. Values lower than the cutoff are considered \(0.0\).

• IClist (file, required for Ideal Chromosome training) – A one-column text file containing the energy interaction values for loci i and j separated by a genomic distance \(d\). The list should be at least of the size \(dend-dinit\). (Default value: None).

• dinit (int, required) – The first neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential for training. (Default value = 3).

• dend (int, required) – The last neighbor in sequence separation (Genomic Distance) to be considered in the Ideal Chromosome potential for training. (Default value = 200).

calc_phi_exp_IC()

Calculates the contact probability as a function of the genomic distance from the experimental Hi-C for the Ideal Chromosome optimization.

calc_phi_exp_types()

Calculates the average of the contact probability for each chromatin type (compartment annotation) from the experimental Hi-C for the Types optimization.

calc_phi_sim_IC()

Calculates the contact probability as a function of the genomic distance from simulations for the Ideal Chromosome optimization.

calc_phi_sim_types()

Calculates the average of the contact probability for each chromatin type (compartment annotation) from simulation for the Types optimization.

get_HiC_exp(HiC, centerRemove=False, centrange=[0, 0], norm=False, cutoff_low=0.0, cutoff_high=1.0, KR=False, neighbors=0)

Receives the experimental Hi-C map (Full dense matrix) in a text format and performs the data normalization from Hi-C frequency/counts/reads to probability.

Parameters

• HiC (file, required) – Experimental Hi-C map (Full dense matrix) in a text format.

• centerRemove (bool, optional) – Whether to set the contact probability of the centromeric region to zero. (Default value: False).

• centrange (list, required if centerRemove = True)) – Range of the centromeric region, i.e., centrange=[i,j], where i and j are the initial and final beads in the centromere. (Default value = [0,0]).

• cutoff (float, optional) – Cutoff value for reducing the noise in the original data. Values lower than the cutoff are considered \(0.0\).

get_HiC_sim()

Calculates the in silico Hi-C map (Full dense matrix).

get_Pearson()

Calculates the Pearson’s Correlation between the experimental Hi-C used as a reference for the training and the in silico Hi-C obtained from the optimization step.

get_PiPj_sim_IC()

Normalizes the cross term of the Hessian by the number of frames in the simulation for the Ideal Chromosome optimization.

get_PiPj_sim_types()

Normalizes the cross term of the Hessian by the number of frames in the simulation for the Types optimization.

get_chrom_seq(filename)

Reads the chromatin sequence as letters of the types/compartments.

Parameters

filename (file, required) –

Chromatin sequence of types file. The first column should contain the locus index. The second column should have the locus type annotation. A template of the chromatin sequence of types file can be found at the Nucleome Data Bank (NDB).

Returns

Returns an array of the sequence of chromatin types.

Return type

\((N,1)\) numpy.ndarray

get_lambdas_IC(exp_map='file.dense', damp=3e-07, write_error=True)

Calculates the Lagrange multipliers for the Ideal Chromosome optimization and returns a array containing the energy values for the IC optimization step. :param exp_map: The experimental Hi-C map with the .dense file. (Default value: file.dense). :type exp_map: file, required :param damp: The learning parameter for the new lambda. (Default value = \(3*10**-7\)). :type damp: float :param dmax: The maximum distance in the sequence separation (Genomic Distance) to be considered for the convergence of the potential interations. (Default value = 200).

The learning parameter for the new lambda. (Default value = \(3*10**-7\)).

Parameters

write_error (boolean) – Flag to write the tolerance and Pearson’s correlation values. (Default value: True).

get_lambdas_types(exp_map, damp=5e-07, write_error=True)

Calculates the Lagrange multipliers of each type-to-type interaction and returns the matrix containing the energy values for the optimization step.

normalize_matrix(matrix)

Normalize the matrix for simulation optimization. Here the first neighbor should have the probability of contact P=1.0.

prob_calculation_IC(state)

Calculates the contact probability matrix and the cross term of the Hessian for the Ideal Chromosome optimization.

prob_calculation_types(state)

Calculates the contact probability matrix and the cross term of the Hessian for the type-to-type interactions optimization.

class OpenMiChroM.Optimization.FullTraining(expHiC, mu=2.0, rc=2.5, cutoff=0.0, reduce=True, pair_h=2, c_h=0.1, pair_l=4, c_l=0.02)

Bases: object

The FullTraining class performs the parameters training for each selected loci pair interaction.

Details about the methodology are decribed in “Zhang, Bin, and Peter G. Wolynes. “Topology, structures, and energy landscapes of human chromosomes.” Proceedings of the National Academy of Sciences 112.19 (2015): 6062-6067.”

The FullTraining class receive a Hi-C matrix (text file) as input. The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 2.0).

• rc (float, required) – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 2.5).

• cutoff (float, optional) – Cutoff value for reducing the noise in the original data. Values lower than the cutoff are considered \(0.0\).

• reduce (bool, optional) – Whether to reduce the number of interactions to be considered in the inversion. If False, it will consider every possible interaction \((N*(N-1)/2)\). If True, it is necessary to give values for the lower and higher cutoffs. (Default value: True).

• pair_h (int, required if reduce = True) – Loci selection to apply the high-resolution cutoff. If pair_h = 2, the interaction in the high-resolution index grid \(2 : 2 : N × 2:2:N\) are subject to a cutoff value c_h, where N is the total number of monomers interactions (Default value = 2).

• c_h (float, required if reduce = True)) – The the high-resolution cutoff. (Default value = 0.1).

• pair_l (int, required if reduce = True)) – Loci selection to apply the high-resolution cutoff. If pair_l = 4, the interaction in the low-resolution index grid \(1:4:N×1:4:N\) are subject to a cutoff value c_l, where N is the total number of monomers interactions (Default value = 4).

• c_l (float, required if reduce = True)) – The the low-resolution cutoff. (Default value = 0.02).

appCutoff(pair_h, c_h, pair_l, c_l)

Applies the cutoff for low- and high-resolution values.

createInitialLambda(sequenceFile, outputPath='.', initialGuess=0.0, baseLine=- 0.2)

getLambdas()

Calculates the Lagrange multipliers of each selected interaction and returns the matrix containing the energy values for the optimization step.

getPearson()

Calculates the Pearson’s Correlation between the experimental Hi-C used as a reference for the training and the in silico Hi-C obtained from the optimization step.

get_indices(hic)

Receives non-zero interaction indices, i.e., the loci pair i and j which interaction will be optimized.

probCalc(state)

Calculates the contact probability matrix and the cross term of the Hessian.

saveLambdas(sequenceFile, data, outputPath, name)

OpenMiChroM.CndbTools

The cndbTools class perform analysis from cndb or ndb - (Nucleome Data Bank) file format for storing an ensemble of chromosomal 3D structures. Details about the NDB/CNDB file format can be found at the Nucleome Data Bank.

class OpenMiChroM.CndbTools.cndbTools

Bases: object

compute_Chirality(xyz, neig_beads=4)

Calculates the Chirality parameter \(\Psi\). Details are decribed in “Zhang, B. and Wolynes, P.G., 2016. Shape transitions and chiral symmetry breaking in the energy landscape of the mitotic chromosome. Physical review letters, 116(24), p.248101.”

Parameters

• xyz (\((frames, beadSelection, XYZ)\) numpy.ndarray, required) – Array of the 3D position of the selected beads for different frames extracted by using the xyz() function.

• neig_beads (int, required) – Number of neighbor beads to consider in the calculation (Default value = 4).

Returns

Returns the Chirality parameter \(\Psi\) for each bead.

Return type

numpy.ndarray

compute_FFT_from_Oij(Oijy, lowcut=1, highcut=500, order=5)

compute_GyrTensorEigs(xyz)

Calculates the eigenvalues of the Gyration tensor: For a cloud of N points with positions: {[xi,yi,zi]},gyr tensor is a symmetric matrix defined as,

gyr= (1/N) * [[sum_i(xi-xcm)(xi-xcm) sum_i(xi-xcm)(yi-ycm) sum_i(xi-xcm)(zi-zcm)],

[sum_i(yi-ycm)(xi-xcm) sum_i(yi-ycm)(yi-ycm) sum_i(yi-ycm)(zi-zcm)], [sum_i(zi-zcm)(xi-xcm) sum_i(zi-zcm)(yi-ycm) sum_i(zi-zcm)(zi-zcm)]]

the three non-negative eigenvalues of gyr are used to define shape parameters like radius of gyration, asphericity, etc

Parameters

xyz (:math:`(frames, beadSelection, XYZ) – TxNx3), required): Array of the 3D position of the selected beads for different frames extracted by using the xyz() function.

Returns

Tx3):

Returns the sorted eigenvalues of the Gyration Tensor.

Return type

numpy.ndarray (dim

compute_MSD(xyz)

Calculates the Mean-Squared Displacement using Fast-Fourier Transform. Uses Weiner-Kinchin theorem to compute the autocorrelation, and a recursion realtion from the following reference: see Sec. 4.2 in Calandrini V, et al. (2011) EDP Sciences (https://doi.org.10.1051/sfn/201112010). Also see this stackoverflow post: https://stackoverflow.com/questions/34222272/computing-mean-square-displacement-using-python-and-fft

Parameters

xyz (:math:`(frames, beadSelection, XYZ) – TxNx3), required): Array of the 3D position of the selected beads for different frames extracted by using the xyz() function.

Returns

NxT):

Returns the MSD of each particle over the trajectory.

Return type

numpy.ndarray (dim

compute_Orientation_OP(xyz, chrom_start=0, chrom_end=1000, vec_length=4)

compute_RDP(xyz, beadSelection=None, radius=20.0, bins=200)

Calculates the RDP - Radial Distribution Probability. Details can be found in the following publications:

• Oliveira Jr., A.B., Contessoto, V.G., Mello, M.F. and Onuchic, J.N., 2021. A scalable computational approach for simulating complexes of multiple chromosomes. Journal of Molecular Biology, 433(6), p.166700.

• Di Pierro, M., Zhang, B., Aiden, E.L., Wolynes, P.G. and Onuchic, J.N., 2016. Transferable model for chromosome architecture. Proceedings of the National Academy of Sciences, 113(43), pp.12168-12173.

Parameters

• xyz (\((frames, XYZ)\) numpy.ndarray, required) – Array of the 3D position of the frames extracted by using the xyz() function.

• beadSelection (\((beadSelection)\) numpy.ndarray) – The index of the beads to be sliced from xyz that you want to compute RDP. Usualy, you can use the internal selection using the dictChromSeq[‘types’] with ‘types’ been the selection that you want.

• radius (float, required) – Radius of the sphere in units of \(\sigma\) to be considered in the calculations. The radius value should be modified depending on your simulated chromosome length. (Default value = 20.0).

• bins (int, required) – Number of slices to be considered as spherical shells. (Default value = 200).

Returns

Returns the radius of each spherical shell in units of \(\sigma\). \((N, 1)\) numpy.ndarray:

Returns the RDP - Radial Distribution Probability for each spherical shell.

Return type

\((N, 1)\) numpy.ndarray

compute_RG(xyz)

Calculates the Radius of Gyration.

Parameters

xyz (:math:`(frames, beadSelection, XYZ) – TxNx3), required): Array of the 3D position of the selected beads for different frames extracted by using the xyz() function.

Returns

Tx1):

Returns the Radius of Gyration in units of \(\sigma\).

Return type

numpy.ndarray (dim

compute_RadNumDens(xyz, dr=1.0, ref='centroid', center=None)

Calculates the radial number density of monomers; which when integrated over the volume (with the appropriate kernel: 4*pi*r^2) gives the total number of monomers.

Parameters

• xyz (:math:`(frames, beadSelection, XYZ) – TxNx3), required): Array of the 3D position of the selected beads for different frames extracted by using the xyz() function.

• dr (float, required) – mesh size of radius for calculating the radial distribution. can be arbitrarily small, but leads to empty bins for small values. bins are computed from the maximum values of radius and dr.

• ref (string) –

  defines reference for centering the disribution. It can take three values:

  ’origin’: radial distance is calculated from the center

  ’centroid’ (default value): radial distributioin is computed from the centroid of the cloud of points at each time step

  ’custom’: user defined center of reference. ‘center’ is required to be specified when ‘custom’ reference is chosen

• center (list of float, len 3) – defines the reference point in custom reference. required when ref=’custom’

Returns

the number density

bins:class:numpy.ndarray:

bins corresponding to the number density

Return type

num_density:class:numpy.ndarray

load(fileName)

Receives the path to cndb or ndb file to perform analysis.

Parameters

fileName (file, required) – Path to cndb or ndb file. If an ndb file is given, it is converted to a cndb file and saved in the same directory.

ndb2cndb(fileName)

Converts an ndb file format to cndb.

Parameters

filename (path, required) – Path to the ndb file to be converted to cndb.

traj2HiC(xyz, mu=3.22, rc=1.78)

Calculates the in silico Hi-C maps (contact probability matrix) using a chromatin dyamics trajectory.

The parameters \(\mu\) (mu) and rc are part of the probability of crosslink function \(f(r_{i,j}) = \frac{1}{2}\left( 1 + tanh\left[\mu(r_c - r_{i,j}\right] \right)\), where \(r_{i,j}\) is the spatial distance between loci (beads) i and j.

Parameters

• mu (float, required) – Parameter in the probability of crosslink function. (Default value = 3.22).

• rc – Parameter in the probability of crosslink function, \(f(rc) = 0.5\). (Default value = 1.78).

xyz(frames=None, beadSelection=None, XYZ=[0, 1, 2])

Get the selected beads’ 3D position from a cndb or ndb for multiple frames.

Parameters

• frames (list, required) – Define which frames to extract the position of the selected beads. The list should include all frames to iterate over, or a iterator that goes over the desired frames. (Default value: None, all frames)

• beadSelection (list of ints, required) – List of beads to extract the 3D position for each frame. The list is defined by beadSelection=[0,1,2,…,N-1]. (Default value: None, all beads)

• XYZ (list, required) – List of the axis in the Cartesian coordinate system that the position of the bead will get extracted for each frame. The list is defined by XYZ=[0,1,2]. where 0, 1 and 2 are the axis X, Y and Z, respectively. (Default value: XYZ=[0,1,2])

Returns

Returns an array of the 3D position of the selected beads for different frames.

Return type

(\(N_{frames}\), \(N_{beads}\), 3) numpy.ndarray

OpenMiChroM.Integrators

The ActiveBrownianIntegrator class is a custom Brownian integrator. Just like all Brownian integrators there are no velocities, there are only forces and displacements. Here we use the velocities as a proxy to keep track of the active force for each monomer. Hence velocity v of a monomer represents the active force f = gamma * v.

class OpenMiChroM.Integrators.ActiveBrownianIntegrator(*args: Any, **kwargs: Any)

Bases: openmm.CustomIntegrator

The ActiveBrownianIntegrator class is a custom Brownian integrator. Just like all Brownian integrators there are no velocities, there are only forces and displacements. Here we use the velocities as a proxy to keep track of the active force for each monomer. Hence velocity v of a monomer represents the active force f = gamma * v.