Tutorial 5: Two-Mark Chromatin with Competitive Binding

In tutorial_4, we introduced a model of chromatin with two epigenetic marks and their reader proteins. We assume that each reader protein binds its respective epigenetic mark without directly interfering with the binding of the other reader. In this tutorial, we will add a constraint on the reader protein binding. We will assume that a maximum of two reader proteins can bind each nucleosome at a time. This demonstration almost directly follows the previous tutorial, with the addition of the constraint on the reader protein binding. The descriptions included in this notebook will only highlight differences associated with the binding constraint.


Import Modules

[1]:

# Built-in modules
import os
import sys
from inspect import getmembers, isfunction

# Third-party modules
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Custom modules
from chromo.binders import get_by_name, make_binder_collection
from chromo.polymers import Chromatin
from chromo.fields import UniformDensityField
import chromo.mc as mc
import chromo.mc.mc_controller as ctrl
from chromo.util.reproducibility import get_unique_subfolder_name
import chromo.util.mu_schedules as ms

Specify Binders

[2]:

# Instantiate HP1 and PRC1 reader proteins
hp1 = get_by_name('HP1')
prc1 = get_by_name('PRC1')

In this demonstration, we will increase the chemical potentials of the reader proteins to force competition for binding sites.

[3]:

# Set the HP1-PRC1 cross interaction to 1.0 kT
hp1.cross_talk_interaction_energy["PRC1"] = 1.0

# Set the HP1 and PRC1 chemical potentials to 1.0 kT
hp1.chemical_potential = 1.0
prc1.chemical_potential = 1.0

[4]:

# Create a binder collection with the HP1 and PRC1 reader proteins
binder_collection = make_binder_collection([hp1, prc1])

Specify the Confinement

[5]:

confine_type = "Spherical"
confine_radius = 200.0

Instantiate Epigenetic Mark Pattern

[6]:

num_beads = 1000
mark_pattern = np.zeros((num_beads, 2), dtype=int)
mark_names = np.array(["H3K9me3", "H3K27me3"])
mark_pattern[:500, 0] = 2
mark_pattern[500:, 1] = 2
print("Shape of mark pattern:", mark_pattern.shape)

Shape of mark pattern: (1000, 2)

Specify Initial Reader Protein Binding States

[7]:

binder_names = np.array([hp1.name, prc1.name])
states = np.zeros((num_beads, 2), dtype=int)

Instantiate the Polymer

We specify the maximum number of binders on each nucleosome using the max_binders named argument when we instantiate our polymer. By default, max_binders=-1 to indicate that there is no direct competition for binding sites. In this case, we will set max_binders=2 to indicate that a maximum of two reader proteins can bind each nucleosome at a time.

[8]:

# Specify the name, number of beads, and bead spacing of the chromatin fiber
name = "Chr"
bead_spacing = np.ones(num_beads - 1) * 16.5

# Instantiate the chromatin fiber
poly = Chromatin.confined_gaussian_walk(
    name,
    num_beads,
    bead_spacing,
    confine_type=confine_type,
    confine_length=confine_radius,
    states=states,
    binder_names=binder_names,
    chemical_mods=mark_pattern,
    chemical_mod_names=mark_names,
    max_binders=2
)

Specify the Uniform Density Field

[9]:

 # Specify the dimensions of the field
n_accessible = int(np.round((63 * confine_radius) / 900))
n_buffer = 2
n_bins_x = n_accessible + n_buffer
x_width = 2 * confine_radius * (1 + n_buffer/n_accessible)
n_bins_y = n_bins_x
y_width = x_width
n_bins_z = n_bins_x
z_width = x_width

# Initialize the uniform density field
udf = UniformDensityField(
    [poly],
    binder_collection,
    x_width,
    n_bins_x,
    y_width,
    n_bins_y,
    z_width,
    n_bins_z,
    confine_type=confine_type,
    confine_length=confine_radius,
    chi=1,
    assume_fully_accessible=1,
    fast_field=0
)

Specify the Simulation Parameters

For the purposes of this demonstration, I will reduce the number of mc steps per snapshot. This will allow us to run the simulation faster. For more rigorous simulations, increase the number of mc steps per snapshot to ensure that the system reaches equilibrium.

[10]:

amp_bead_bounds, amp_move_bounds = mc.get_amplitude_bounds(
    polymers = [poly]
)

[11]:

# Specify the number of snapshots and the number of MC steps to attempt per snapshot
num_snapshots = 200
mc_steps_per_snapshot = 200  # Reduce the number of MC steps per snapshot to speed up the simulation
# TODO: If you want to run a more rigorous simulation, increase the number of MC steps per snapshot

[12]:

# Create a list of simulated annealing schedules, which are defined in another file
schedules = [func[0] for func in getmembers(ms, isfunction)]

# Select a schedule for slowly adding HP1 into the system
select_schedule = "linear_step_for_negative_cp"
mu_schedules = [
    ms.Schedule(getattr(ms, func_name)) for func_name in schedules
]
mu_schedules = [sch for sch in mu_schedules if sch.name == select_schedule]

Run the Simulation

[ ]:

# Specify the output directory
out_dir = "output_demo"

# Run the simulation
polymers = mc.polymer_in_field(
    [poly],
    binder_collection,
    udf,
    mc_steps_per_snapshot,
    num_snapshots,
    amp_bead_bounds,
    amp_move_bounds,
    output_dir=out_dir,
    mu_schedule=mu_schedules[0],
)