{ "cells": [ { "cell_type": "markdown", "id": "2b296fce9a4c6c8b", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "# Tutorial 5: Two-Mark Chromatin with Competitive Binding\n", "\n", "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. " ] }, { "cell_type": "markdown", "id": "a9e7273045e1037f", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Import Modules" ] }, { "cell_type": "code", "execution_count": 1, "id": "4bececc8d946b41d", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:09.652268Z", "start_time": "2024-07-03T21:45:08.831113Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Built-in modules\n", "import os\n", "import sys\n", "from inspect import getmembers, isfunction\n", "\n", "# Third-party modules\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "\n", "# Custom modules\n", "from chromo.binders import get_by_name, make_binder_collection\n", "from chromo.polymers import Chromatin\n", "from chromo.fields import UniformDensityField\n", "import chromo.mc as mc\n", "import chromo.mc.mc_controller as ctrl\n", "from chromo.util.reproducibility import get_unique_subfolder_name\n", "import chromo.util.mu_schedules as ms" ] }, { "cell_type": "markdown", "id": "45e0236b607bbe55", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Specify Binders" ] }, { "cell_type": "code", "execution_count": 2, "id": "98a4a592d7643e5e", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:11.753345Z", "start_time": "2024-07-03T21:45:11.751467Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Instantiate HP1 and PRC1 reader proteins\n", "hp1 = get_by_name('HP1')\n", "prc1 = get_by_name('PRC1')" ] }, { "cell_type": "markdown", "id": "289fa3b6d89ca428", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "In this demonstration, we will increase the chemical potentials of the reader proteins to force competition for binding sites." ] }, { "cell_type": "code", "execution_count": 3, "id": "3f7232ce6ee1e13b", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:12.724172Z", "start_time": "2024-07-03T21:45:12.722269Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Set the HP1-PRC1 cross interaction to 1.0 kT\n", "hp1.cross_talk_interaction_energy[\"PRC1\"] = 1.0\n", "\n", "# Set the HP1 and PRC1 chemical potentials to 1.0 kT\n", "hp1.chemical_potential = 1.0\n", "prc1.chemical_potential = 1.0" ] }, { "cell_type": "code", "execution_count": 4, "id": "56c348cd85090d35", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:13.258043Z", "start_time": "2024-07-03T21:45:13.248770Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Create a binder collection with the HP1 and PRC1 reader proteins\n", "binder_collection = make_binder_collection([hp1, prc1])" ] }, { "cell_type": "markdown", "id": "ff74fb440f99be89", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Specify the Confinement" ] }, { "cell_type": "code", "execution_count": 5, "id": "3221ef96e1a4e0c6", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:14.144582Z", "start_time": "2024-07-03T21:45:14.142848Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "confine_type = \"Spherical\"\n", "confine_radius = 200.0" ] }, { "cell_type": "markdown", "id": "3b3743be8c6064ed", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Instantiate Epigenetic Mark Pattern" ] }, { "cell_type": "code", "execution_count": 6, "id": "f3c443113199c929", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:15.061185Z", "start_time": "2024-07-03T21:45:15.058926Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Shape of mark pattern: (1000, 2)\n" ] } ], "source": [ "num_beads = 1000\n", "mark_pattern = np.zeros((num_beads, 2), dtype=int)\n", "mark_names = np.array([\"H3K9me3\", \"H3K27me3\"])\n", "mark_pattern[:500, 0] = 2\n", "mark_pattern[500:, 1] = 2\n", "print(\"Shape of mark pattern:\", mark_pattern.shape)" ] }, { "cell_type": "markdown", "id": "b0394d4ca8eda3be", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Specify Initial Reader Protein Binding States" ] }, { "cell_type": "code", "execution_count": 7, "id": "d8f4eeb964d1b597", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:16.190355Z", "start_time": "2024-07-03T21:45:16.188308Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "binder_names = np.array([hp1.name, prc1.name])\n", "states = np.zeros((num_beads, 2), dtype=int)" ] }, { "cell_type": "markdown", "id": "d39844c4211409a5", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Instantiate the Polymer\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 8, "id": "bfd8802c9ded69f6", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:17.701873Z", "start_time": "2024-07-03T21:45:17.625099Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Specify the name, number of beads, and bead spacing of the chromatin fiber\n", "name = \"Chr\"\n", "bead_spacing = np.ones(num_beads - 1) * 16.5\n", "\n", "# Instantiate the chromatin fiber\n", "poly = Chromatin.confined_gaussian_walk(\n", " name,\n", " num_beads,\n", " bead_spacing,\n", " confine_type=confine_type,\n", " confine_length=confine_radius,\n", " states=states,\n", " binder_names=binder_names,\n", " chemical_mods=mark_pattern,\n", " chemical_mod_names=mark_names,\n", " max_binders=2\n", ")" ] }, { "cell_type": "markdown", "id": "ef7a29ca994c6fb1", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Specify the Uniform Density Field" ] }, { "cell_type": "code", "execution_count": 9, "id": "a190adaba56bf3f9", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:18.729490Z", "start_time": "2024-07-03T21:45:18.717931Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ " # Specify the dimensions of the field\n", "n_accessible = int(np.round((63 * confine_radius) / 900))\n", "n_buffer = 2\n", "n_bins_x = n_accessible + n_buffer\n", "x_width = 2 * confine_radius * (1 + n_buffer/n_accessible)\n", "n_bins_y = n_bins_x\n", "y_width = x_width\n", "n_bins_z = n_bins_x\n", "z_width = x_width\n", "\n", "# Initialize the uniform density field\n", "udf = UniformDensityField(\n", " [poly],\n", " binder_collection,\n", " x_width,\n", " n_bins_x,\n", " y_width,\n", " n_bins_y,\n", " z_width,\n", " n_bins_z,\n", " confine_type=confine_type,\n", " confine_length=confine_radius,\n", " chi=1,\n", " assume_fully_accessible=1,\n", " fast_field=0\n", ")" ] }, { "cell_type": "markdown", "id": "758cdf2a6f6a8233", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Specify the Simulation Parameters\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 10, "id": "48603189e99e4c4e", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:21.502282Z", "start_time": "2024-07-03T21:45:21.499354Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "amp_bead_bounds, amp_move_bounds = mc.get_amplitude_bounds(\n", " polymers = [poly]\n", ")" ] }, { "cell_type": "code", "execution_count": 11, "id": "f6c5c93e5153b197", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:21.955525Z", "start_time": "2024-07-03T21:45:21.952690Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Specify the number of snapshots and the number of MC steps to attempt per snapshot\n", "num_snapshots = 200\n", "mc_steps_per_snapshot = 200 # Reduce the number of MC steps per snapshot to speed up the simulation\n", "# TODO: If you want to run a more rigorous simulation, increase the number of MC steps per snapshot" ] }, { "cell_type": "code", "execution_count": 12, "id": "4c64dc423fe0a93", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T21:45:22.414512Z", "start_time": "2024-07-03T21:45:22.412307Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Create a list of simulated annealing schedules, which are defined in another file\n", "schedules = [func[0] for func in getmembers(ms, isfunction)]\n", "\n", "# Select a schedule for slowly adding HP1 into the system\n", "select_schedule = \"linear_step_for_negative_cp\"\n", "mu_schedules = [\n", " ms.Schedule(getattr(ms, func_name)) for func_name in schedules\n", "]\n", "mu_schedules = [sch for sch in mu_schedules if sch.name == select_schedule]" ] }, { "cell_type": "markdown", "id": "e719d3a194dd18c7", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Run the Simulation" ] }, { "cell_type": "code", "execution_count": null, "id": "97b286e8c133d809", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T22:22:58.623260Z", "start_time": "2024-07-03T21:45:23.414385Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "# Specify the output directory\n", "out_dir = \"output_demo\"\n", "\n", "# Run the simulation\n", "polymers = mc.polymer_in_field(\n", " [poly],\n", " binder_collection,\n", " udf,\n", " mc_steps_per_snapshot,\n", " num_snapshots,\n", " amp_bead_bounds,\n", " amp_move_bounds,\n", " output_dir=out_dir,\n", " mu_schedule=mu_schedules[0],\n", ")" ] }, { "cell_type": "markdown", "id": "d6e94d42c83c75cb", "metadata": { "collapsed": false, "jupyter": { "outputs_hidden": false } }, "source": [ "|\n", "\n", "#### Print the Resulting Binding Patterns\n", "\n", "In this demonstration, we will only highlight the differences between the HP1 and PRC1 binding patterns to demonstrate that only two reader proteins can bind each nucleosome at a time." ] }, { "cell_type": "code", "execution_count": 14, "id": "1c3900c9cc2f6d0a", "metadata": { "ExecuteTime": { "end_time": "2024-07-03T22:22:58.628528Z", "start_time": "2024-07-03T22:22:58.624584Z" }, "collapsed": false, "jupyter": { "outputs_hidden": false } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "HP1 binding pattern: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 0 2 2 2 2 2 2 0 2 0 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 0 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 2 2 2 2 2 2 2 2 2\n", " 2 1 2 2 2 2 2 1 2 2 2 2 2 2 2 0 0 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 1 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 1 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 1 0 2 2 2 2 0 2 2 2 2 2 2 0 2 2 2 2 2 2\n", " 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 0 2 2 2 2 0 2 2 2 2 2 2 2 0\n", " 2 2 2 2 2 2 2 2 2 2 2 1 2 0 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 1 0 0 0 0 0 2 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 2 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 2 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0\n", " 0 0 2 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 2\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0\n", " 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 0 2 0 0 0 0 0\n", " 0]\n", "\n", "PRC1 binding pattern: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 2 2 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 2 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 2 0 0 0 0 2 0 0 0 0 0 0 0 2\n", " 0 0 0 0 0 0 0 0 0 0 0 1 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n", " 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 1 2 2 2 1 2 0 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 0 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n", " 2 2 0 2 1 2 2 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2\n", " 2 2 0 2 2 2 2 0 2 2 2 2 2 0 2 2 2 2 2 2 2 2 1 0 2 2 2 2 2 2 2 2 2 2 2 2 0\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2\n", " 0 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 2 2 2 2 2 2 2 2 2 2\n", " 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 1 0 2 2 2 2 2 2 2 2 0 2 2 2 2 2\n", " 2]\n" ] } ], "source": [ "HP1_pattern = np.asarray(polymers[0].states[:, 0]).flatten()\n", "PRC1_pattern = np.asarray(polymers[0].states[:, 1]).flatten()\n", "\n", "print(\"HP1 binding pattern:\", HP1_pattern)\n", "print()\n", "print(\"PRC1 binding pattern:\", PRC1_pattern)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.7" } }, "nbformat": 4, "nbformat_minor": 5 }