chromo.mc package

Modules in this package operate and adapt the Monte Carlo simulation.

chromo.mc.mc_controller module

Control MC simulation parameters.

class chromo.mc.mc_controller.Controller(mc_adapter: MCAdapter, bead_amp_bounds: Tuple[int, int], move_amp_bounds: Tuple[float, float])

Bases: ABC

Class representation of controller for MC simulation parameters.

abstract update_amplitudes()

Update bead selection and move amplitudes.

abstract update_bead_amplitude()

Update the bead selection amplitude.

Use acceptance rate self.move.acceptance_tracker.acceptance_rate to update the move amplitude.

abstract update_move_amplitude()

Update the move amplitude.

Use acceptance rate self.move.acceptance_tracker.acceptance_rate to update the move amplitude.

class chromo.mc.mc_controller.NoControl(mc_adapter: MCAdapter, bead_amp_bounds: Tuple[int, int], move_amp_bounds: Tuple[float, float])

Bases: Controller

Class representation of a trivial, non-acting MC controller.

update_amplitudes()

Trivially maintain current move and bead amplitudes.

update_bead_amplitude()

Trivially maintain current bead selection amplitude.

update_move_amplitude()

Trivially maintain current move amplitude.

class chromo.mc.mc_controller.SimpleControl(mc_adapter: MCAdapter, bead_amp_bounds: Tuple[int, int], move_amp_bounds: Tuple[float, float])

Bases: Controller

Apply factor changes to move or bead amplitude based on acceptance rate.

A factor change is applied to the move amplitude based on the acceptance rate – if the acceptance rate is less than a setpoint (default 0.5), move amplitude is increased; otherwise, the move amplitude is decreased. Once the move amplitude reaches an upper or lower threshold, rather than changing the move amplitude, the bead amplitude is increased or decreased by a certain length, and the move amplitude is reset to its lowest value.

update_amplitudes(setpoint_acceptance: Optional[float] = 0.5, move_adjust_factor: Optional[float] = 0.95, num_delta_beads: Optional[int] = 1)

Update move and bead amplitudes using a simple controller.

Parameters:

• setpoint_acceptance (Optional[float]) – Optional acceptance rate setpoint (default = 0.5)

• move_adjust_factor (Optional[float]) – Factor by which to multiply or divide move amplitude in response to current acceptance rate, between 0 and 1 (default = 0.95)

• num_delta_beads (Optional[int]) – Number of beads by which to adjust bead amplitude (default = 1)

update_bead_amplitude(increase: bool, num_delta_beads: Optional[int] = 1)

Update the bead amplitude based on the acceptance rate.

Parameters:

• increase (bool) – True or false indicator for whether to increase bead amplitude

• num_delta_bead (Optional[int]) – Number of beads by which to increase bead amplitude (default = 1)

Returns:

prop_bead_amp – Proposed bead amplitude after adjustment

Return type:

int

update_move_amplitude(setpoint_acceptance: Optional[float] = 0.5, move_adjust_factor: Optional[float] = 0.95, num_delta_beads: Optional[int] = 1)

Update move amplitude based on acceptance rate.

Parameters:

• setpoint_acceptance (Optional[float]) – Optional acceptance rate setpoint (default = 0.5)

• move_adjust_factor (Optional[float]) – Factor by which to multiply or divide move amplitude in response to current acceptance rate, between 0 and 1 (default = 0.95)

• num_delta_beads (Optional[int]) – Number of beads by which to adjust bead amplitude (default = 1)

chromo.mc.mc_controller.all_moves(log_dir: str, bead_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], move_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], controller: ~typing.Optional[~chromo.mc.mc_controller.Controller] = <class 'chromo.mc.mc_controller.NoControl'>) → List[Controller]

Generate a list of controllers for all adaptable MC moves (single poly).

Parameters:

• log_dir (str) – Path to the directory in which to save log files

• bead_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of bead selection amplitude bounds where keys specify the name of the MC move and the values specify the selection amplitude bounds for the move in the form (lower bound, upper bound)

• move_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of maximum move amplitude bounds where keys specify the name of the MC move and the values specify the move amplitude bounds for the move in the form (lower bound, upper bound)

• controller (Optional[Controller]) – Bead and move amplitude controller (default = NoControl)

Return type:

List of controllers for all adaptable MC moves.

chromo.mc.mc_controller.all_moves_except_binding_state(log_dir: str, bead_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], move_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], controller: ~typing.Optional[~chromo.mc.mc_controller.Controller] = <class 'chromo.mc.mc_controller.NoControl'>) → List[Controller]

Generate list of controllers for physical MC moves (single poly).

Parameters:

• log_dir (str) – Path to the directory in which to save log files

• bead_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of bead selection amplitude bounds where keys specify the name of the MC move and the values specify the selection amplitude bounds for the move in the form (lower bound, upper bound)

• move_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of maximum move amplitude bounds where keys specify the name of the MC move and the values specify the move amplitude bounds for the move in the form (lower bound, upper bound)

• controller (Optional[Controller]) – Bead and move amplitude controller (default = NoControl)

Return type:

List of controllers for all adaptable MC moves.

chromo.mc.mc_controller.only_binding_move(log_dir: str, bead_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], move_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], controller: ~typing.Optional[~chromo.mc.mc_controller.Controller] = <class 'chromo.mc.mc_controller.NoControl'>) → List[Controller]

Generate list containing controller for the binding move (single poly).

Parameters:

• log_dir (str) – Path to the directory in which to save log files

• bead_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of bead selection amplitude bounds where keys specify the name of the MC move and the values specify the selection amplitude bounds for the move in the form (lower bound, upper bound)

• move_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of maximum move amplitude bounds where keys specify the name of the MC move and the values specify the move amplitude bounds for the move in the form (lower bound, upper bound)

• controller (Optional[Controller]) – Bead and move amplitude controller (default = NoControl)

Return type:

List of controllers for all adaptable MC moves.

chromo.mc.mc_controller.specific_move(move_fxn, log_dir: str, bead_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], move_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], controller: ~typing.Optional[~chromo.mc.mc_controller.Controller] = <class 'chromo.mc.mc_controller.NoControl'>) → List[Controller]

Generate a list of controllers for a specific MC move.

Parameters:

• move_fxn (Callable) – Monte Carlo move function (from move_funcs.pyx)

• log_dir (str) – Path to the directory in which to save log files

• bead_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of bead selection amplitude bounds where keys specify the name of the MC move and the values specify the selection amplitude bounds for the move in the form (lower bound, upper bound)

• move_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of maximum move amplitude bounds where keys specify the name of the MC move and the values specify the move amplitude bounds for the move in the form (lower bound, upper bound)

• controller (Optional[Controller]) – Bead and move amplitude controller (default = NoControl)

Return type:

List of controllers for all adaptable MC moves.

chromo.mc.mc_controller.specific_moves(move_fxns: ~typing.List[~typing.Callable], log_dir: str, bead_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], move_amp_bounds: ~typing.Dict[str, ~typing.Tuple[float, float]], controller: ~typing.Optional[~chromo.mc.mc_controller.Controller] = <class 'chromo.mc.mc_controller.NoControl'>) → List[Controller]

Generate a list of controllers for all adaptable physical MC moves.

Parameters:

• move_fxns (List[Callable]) – List of Monte Carlo move functions (from move_funcs.pyx)

• log_dir (str) – Path to the directory in which to save log files

• bead_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of bead selection amplitude bounds where keys specify the name of the MC move and the values specify the selection amplitude bounds for the move in the form (lower bound, upper bound)

• move_amp_bounds (Dict[str, Tuple[float, float]]) – Dictionary of maximum move amplitude bounds where keys specify the name of the MC move and the values specify the move amplitude bounds for the move in the form (lower bound, upper bound)

• controller (Optional[Controller]) – Bead and move amplitude controller (default = NoControl)

Return type:

List of controllers for all adaptable MC moves.

chromo.mc.mc_sim module

Routines for performing Monte Carlo simulations.

chromo.mc.mc_sim.mc_sim()

Perform Monte Carlo simulation.

Pseudocode

Repeat for each Monte Carlo step:

Repeat for each adaptable move:

If active, apply move to each polymer

Notes

To report process time for each interval of MC steps, add the following:

if (k+1) % 50000 == 0:

print(“MC Step “ + str(k+1) + “ of “ + str(num_mc_steps)) print(

“Time for previous 50000 MC Steps (in seconds): “, round(

process_time()-t1_start, 2

)

) t1_start = process_time()

TODO: To make the simulator more reproducible, consider setting a random seed for the CYTHON code. Note, this uses srand (in libc.stdlib). However, be very careful not to accidentally set the random seed in the wrong place so as not to disrupt the randomness of the simulation. This is a potential source of error.

param polymers:

Polymers affected by Monte Carlo simulation

type polymers:

List[Polymer]

param readerproteins:

Specification of reader proteins on polymer

type readerproteins:

List[ReaderProteins]

param num_mc_steps:

Number of Monte Carlo steps to take between save points

type num_mc_steps:

int

param mc_move_controllers:

List of controllers for active MC moves in simulation

type mc_move_controllers:

List[Controller]

param field:

Field affecting polymer in Monte Carlo simulation

type field:

FieldBase (or subclass)

param mu_adjust_factor:

Adjustment factor applied to the chemical potential in response to simulated annealing

type mu_adjust_factor:

double

param random_seed:

Randoms seed with which to initialize simulation

type random_seed:

int

chromo.mc.mc_sim.mc_step()

Compute energy change and determine move acceptance.

Notes

Get the proposed state of the polymer. Calculate the total (polymer + field) energy change associated with the move. Accept or reject the move based on the Metropolis Criterion.

After evaluating change in energy (from the polymer and the field), the try-except statement checks for RuntimeWarning if the change in energy gets too large.

Parameters:

• adaptible_move (MCAdapter) – Move applied at particular Monte Carlo step

• poly (PolymerBase) – Polymer affected by move at particular Monte Carlo step

• readerproteins (List[ReaderProteins]) – Reader proteins affecting polymer configuration

• field (FieldBase (or subclass)) – Field affecting polymer in Monte Carlo step

• active_field (bool) – Indicator of whether or not the field is active for the polymer

• update_distances (bint) – Update pairwise distances between beads – only relevant if the polymer is an instance of DetailedChromatinWithSterics, which tracks pairwise distances between beads

chromo.mc.move_funcs module

Functions for performing Monte Carlo transformations.

Notes

TODO: No need to make a numpy inds array when the inds are continuous!

chromo.mc.move_funcs.change_binding_state()

Flip the binding state of a reader protein.

Notes

Before applying the move, check that the polymer has any reader proteins at all.

Begin the move by identifying the number of binding sites that each bead has for the particular reader protein. Then randomly select a bead in the chain. Select a second bead from a two sided decaying exponential distribution around the index of the first bead. Replace the binding state of the selected beads.

In our model, we do not track the binding state of individual tails. We care only about how many tails are bound. Therefore, for each move, we will generate a random order of bound and unbound tails and will flip the state of the first M tails in that order.

This function assumes that polymers are formed from a single type of bead object. if different types of bead objects exist, then the first two if-statements that reference bead[0] will need to be generalized.

Parameters:

• polymer (poly.PolymerBase) – Polymer object on which the move is applied

• amp_move (double) – Magnitude of the MC move; this attribute is included simply for consistency with the other MC move functions and does not affect the binding state move

• amp_bead (int) – Maximum range of beads to which a binding state swap will take palce

Returns:

Indices of M beads affected by the MC move

Return type:

array_like (M,) of long

chromo.mc.move_funcs.crank_shaft()

Rotate section of polymer around axis connecting two end beads.

Notes

The crank shaft move affects a continuous range of beads; the internal configuration is unaffected by the crank shaft move. Therefore, when evaluating polymer energy change associated with the move, it is sufficient to look only at the ends of the affected segment.

Begin by randomly selecting a starting and ending index for the crank-shaft move based on the bead amplitude. Then generate a random rotation angle for the move based on the move amplitude. Obtain the axis of rotation from the change in position between the starting and ending beads. Finally, generate the rotation matrix and obtain trial positions and tangents for evaluation.

TODO: If greater speed boost is needed, return the length of inds.

Parameters:

• polymer (poly.PolymerBase) – Polymer object

• amp_move (double) – Maximum amplitude (rotation angle) of the crank-shaft move

• amp_bead (int) – Maximum amplitude (number) of beads affected by the crank-shaft move

Returns:

Indices of M beads affected by the MC move

Return type:

array_like (M,) of long

chromo.mc.move_funcs.end_pivot()

Randomly rotate segment from end of polymer about random axis.

Notes

Begin by selecting a random rotation angle based on the move amplitude. The isolate a random subset of beads on either the LHS or RHS of the polymer; the size of this subset is affected by the bead amplitude. Update the transformation matrix to reflect the end-pivot of this polymer segment about a random rotation axis. Apply the transformation matrix to the polymer position and orientations to generate trial positions and orientations dictated by the move.

TODO: Remove the modulo operation when determining whether or not to rotate the LHS of the polymer. Set rotate_lhs = rand(). Change the rotate_lhs condition to if rotate_lhs < RAND_MAX/2.

Parameters:

• polymer (poly.PolymerBase) – Polymer object

• amp_move (double) – Maximum amplitude (rotation angle) of the end-pivot move

• amp_bead (int) – Maximum amplitude (number) of beads affected by the end-pivot move

Returns:

Indices of M beads affected by the MC move

Return type:

array_like (M,) of long

chromo.mc.move_funcs.full_chain_rotation()

Rotate an entire polymer about an arbitrary axis.

Notes

This move does not change a polymer’s internal configurational (elastic) energy and is only relevant to simulations involving more than one polymer.

The rotation takes place about a random axis, sampled uniformally from the unit sphere. The fulcrum of the rotation is a random bead of the polymer.

Parameters:

• polymer (poly.PolymerBase) – Polymer object

• amp_move (double) – Maximum amplitude (rotation angle) of the rotation move

• amp_bead (long) – Maximum amplitude of the window of beads selected for the move; this attribute is included simply to match the formatting of the other move functions and does not affect the full chain rotation

Returns:

Indices of beads affected by the MC move (will always be entire polymer)

Return type:

array_like (N,) of long

chromo.mc.move_funcs.full_chain_translation()

Translate an entire polymer in an arbitrary direction.

Notes

This move does not change a polymer’s internal configurational (elastic) enregy and is only relevant to simulations involving more than one polymer.

The translation occurs in a random direction, sampled uniformally from a unit sphere.

Parameters:

• polymer (poly.PolymerBase) – Polymer object

• amp_move (double) – Maximum amplitude (translation distance) of the slide move

• amp_bead (long) – Maximum amplitude of the window of beads selected for the move; this attribute is included simply to match the formatting of the other move functions and does not affect the full chain translation

Returns:

Indices of beads affected by the MC move (will always be entire polymer)

Return type:

array_like (N,) of long

chromo.mc.move_funcs.slide()

Randomly translate a segment of the polymer.

Notes

Randomly set the translation distance based on the slide move amplitude. Then randomly pick a direction by sampling uniformally from a unit sphere. Split the translation distance into x, y, z components using the direction. Select a random segment of beads to move based on the bead amplititude, and identify the affected indices. Update and apply the transformation matrix for the translation to the affected beads.

There is a difference in how we define the move amplitude between this code and the original FORTRAN codebase. In the original codebase, the move amplitude specifies the maximum translation in each dimension, while in this code, the move amplitude specifies the maximum magnitude of translation.

Parameters:

• polymer (poly.PolymerBase) – Polymer object

• amp_move (double) – Maximum amplitude (translation distance) of the slide move

• amp_bead (int) – Maximum amplitude (number) of beads affected by the slide move

Returns:

Indices of M beads affected by the MC move

Return type:

array_like (M,) of long

chromo.mc.move_funcs.tangent_rotation()

Random bead rotation for a random selection of beads.

Notes

Select a random rotation angle. Then select some random number of beads to rotate, requiring that at least one bead be selected if the move is on. Call rotate_selected_beads to generate a random axis of rotation for each bead, update the transformation matrix, and apply the rotation to each selected bead.

Parameters:

• polymer (poly.PolymerBase) – Polymer object affected by tangent rotation

• amp_move (double) – Range of angles allowed for single tangent rotation move

• amp_bead (int) – Number of beads to randomly rotate

Returns:

Indices of M beads affected by the MC move

Return type:

array_like (M,) of long

chromo.mc.move_funcs.transform_r_t3_t2()

Apply transformation operations to r_trial, t3_trial, t2_trial.

Parameters:

• polymer (poly.PolymerBase) – Polymer affected by the transformation (which is also storing the transformation matrix).

• inds (long[:]) – Bead indices affected by the MC move

• n_inds (long) – Number of indices affected by the move

chromo.mc.moves module

Utilities for proposing Monte Carlo moves.

This module includes an MC adapter class which will drive move adaption and functions applying transformations made by each move.

class chromo.mc.moves.Bounds

Bases: object

Class representation of move or bead amplitude bounds.

to_csv()

Save Polymer object to CSV file as DataFrame.

to_dataframe()

Express the Bounds using a dataframe.

to_file()

Synonym for to_csv to conform to make_reproducible spec.

class chromo.mc.moves.MCAdapter

Bases: object

Track success rate and adjust parameters for Monte Carlo moves.

In order to ensure that a Monte Carlo simulation equilibrates as quickly as possible, we provide an automated system for tuning the “aggressiveness” of each Monte Carlo move in order to optimize the tradeoff between

1. The move being large enough that it does not take too many moves to equilibrate, and

2. The move not being so large that it will never be accepted, or the computational cost becomes prohibitive (especially in the case when the computational cost scales super-linearly).

All Monte-Carlo moves we currently use for our polymers have two “amplitudes”. The “bead” amplitude controls how many beads are typically affected by the move, and the “move” amplitude controls how “far” those beads are moved. Typically, actually bead counts and move sizes are selected randomly, so these amplitudes will correspond to averages.

Success rate will be tracked by a AcceptanceTracker object specified in the chromo/util/mc_stat.py module. The performance tracker will indicate an overall acceptance rate based on the total numbers of moves attempted and accepted. A running acceptance rate will also be maintained, which decays the weight of historic values using a decay factor.

In the future, we want to do this optimization by looking at actual metrics based on the simulation output, but for now only implementation of an MCAdapter (MCAdaptSuccessRate) handles this optimization by simply trying to maintain a fixed target success rate.

accept()

Update polymer with new state and update proposal stats.

Update all elements of poly for which proposed state is not None. Log the move acceptance/rejection in the move acceptance tracker.

replace_none is outdated with cython implementation of move_funcs

Parameters:

• poly (PolymerBase) – Polymer affected by the MC move

• dE (float) – Change in energy associated with the move

• inds (long[:]) – Indices of beads affected by the MC move

• n_inds (long) – Number of beads affected by a move

• log_move (bint) – Indicator for whether (1) or not (0) to log the move to a list later outputted to a CSV file

• log_update (bint) – Indicator for whether (1) or not (2) to record the updated acceptance rate after the MC move

• update_distances (bint) – Update pairwise distances between beads – only relevant if the polymer is an instance of DetailedChromatinWithSterics, which tracks pairwise distances between beads.

reject()

Reject a proposed Monte Carlo move.

Log the rejected move in the acceptance tracker.

Parameters:

• poly (PolymerBase) – Polymer affected by the MC move

• dE (float) – Change in energy associated with the move

• inds (long[:]) – Indices of the beads that were affected by the move

• n_inds (long) – Number of beads that were affected by the move

• log_move (bint) – Indicator for whether (1) or not (0) to log the move to a list later outputted to a CSV file

• log_update (bint) – Indicator for whether (1) or not (2) to record the updated acceptance rate after the MC move

• update_distances (bint) – Update pairwise distances between beads – only relevant if the polymer is an instance of DetailedChromatinWithSterics, which tracks pairwise distances between beads.

chromo.mc.moves.all_moves()

Generate a list of all adaptable MC move objects.

NOTE: Use all_moves function in chromo.mc.mc_controller module to create list of controllers for all moves. This function only creates list of all moves, which may not be compatible with mc_sim function in chromo.mc.mc_sim module.

Parameters:

log_dir (str) – Path to the directory in which to save log files

Return type:

List of all adaptable MC move objects

chromo.mc.__init__

Monte Carlo simulations of a discrete wormlike chain.

This module runs Monte Carlo simulations for the reconfiguration of multiple discrete wormlike chains when placed in a user-specified field and labeled by multiple epigenetic marks.

chromo.mc.continue_polymer_in_field_simulation(polymer_class, binders: List[ReaderProtein], field: FieldBase, output_dir: str, num_save_mc: int, num_saves: int, mc_move_controllers: Optional[List[Controller]] = None, random_seed: Optional[int] = 0)

Continue a simulation of a polymer in a field.

Parameters:

• polymer_class (PolymerBase) – Class of polymers for which to continue the simulation

• binders (List[ReaderProteins]) – List of reader proteins active on polymer

• field (FieldBase) – The discretization of space in which to simulate the polymers

• output_dir (str) – Output directory from which to continue the simulation

• num_save_mc (int) – Number of steps per snapshot in continued simulation

• num_saves (int) – Number of additional save points to collect

• mc_move_controllers (Optional[List[Controller]]) – Controllers for monte carlo moves desired; default of None activates SimpleControl for all MC moves

• random_seed (Optional[SEED]) – Random seed for replication of simulation (default = 0)

chromo.mc.get_amplitude_bounds(polymers: List[PolymerBase]) → Tuple[Dict[str, Tuple[int, int]], Dict[str, Tuple[float, float]]]

Get lower and upper bound for bead selection and move amplitudes.

Parameters:

polymers (List[PolymerBase]) – List of polymers involved in simulation

Returns:

• Dict[str, Tuple[int, int]] – Dictionary of bead selection bounds for each move type, where keys are the names of the move types and values are tuples in the form (lower bound, upper bound)

• Dict[str, Tuple[int, int]] – Dictionary of move amplitude bounds for each move type, where keys are the names of the move types and values are tuples in the form (lower bound, upper bound)