Random walk, brownian motion, and diffusion

This notebook illustrates some interesting applications of random walk models in biology

import numpy as np
from scipy import stats
from matplotlib import pyplot as plt

I. Polymer models of noncondensed chromosomes

For this toy model of chromosomes (the Rouse model), import a pre-existing python package "rouse"

Link: https://github.com/OpenTrajectoryAnalysis/rouse

The basic logic of the Rouse model is very simple. Each monomer wiggles in space due to thermal fluctuation ("random walk" dynamics). In addition to this, each pair of bonded monomers attract each other like a spring with a vanishingly small relaxation length ("non-random" dynamics)

import rouse

Two parameters to initialize the model: $N$ is the number of monomers, and $d$ is space dimension (default=3)

model = rouse.Model(N=10, d=2)
print(model)

rouse.Model(N=10, D=1.0, k=1.0, d=2)

# Define a function to plot a few polymer conformations
# we'll use this repeatedly below
def plot_conformations(conformations):
    n = len(conformations)
    fig, axs = plt.subplots(1, n, figsize=[3*n, 3])

    for i, ax in enumerate(axs):
        conf = conformations[i]
        ax.plot(conf[:, 0], conf[:, 1],
                marker='o',
                color=f"C{i}",
               )
        ax.axis('square')
        ax.axis('off')

    plt.show()

Here we show some equilibrium conformations of the polymer chain under the Rouse model

model = rouse.Model(N=10, d=3)
conformations = [model.conf_ss() for _ in range(5)]
plot_conformations(conformations)

model = rouse.Model(N=100, d=3)
conformations = [model.conf_ss() for _ in range(5)]
plot_conformations(conformations)

model = rouse.Model(N=1000, d=3)
conformations = [model.conf_ss() for _ in range(5)]
plot_conformations(conformations)

The equilibrium conformation of the Rouse model is a collapsed polymer. Let's inspect the collapse process when we evolve a stretched chain freely

N_monomer = 100
Δt = 200
n_steps=4

model = rouse.Model(N=N_monomer)


init = np.array([np.arange(N_monomer), np.arange(N_monomer)]).T # diagonally stretched linear conformation


t = np.arange(0,Δt*(n_steps+1),Δt)

for _ in range(1):
    
    # Get time trace of conformations
    conformations = [init]
    for dt in np.diff(t):
        current = conformations[-1]
        evolved = model.evolve(current, dt=dt)
        conformations.append(evolved)

    # Plot
    n = len(t) # number of subplots
    fig, axs = plt.subplots(1, n,
                            figsize=[n*3, 3],
                            sharex=True,
                            sharey=True,
                           )

    for i, (ax, conf, t_cur) in enumerate(zip(axs, conformations, t)):
        ax.plot(conf[:, 0], conf[:, 1],
                marker='o',
                color=f"C{i}",
               )
        ax.axis('square')
        ax.axis('off')
        ax.set_title(f"t = {t_cur}")

    plt.show()

Despite the collapse, neighboring monomers are more likely to "contact each other". That's because they are more tightly linked and so their spatial positions are more correlated. In modern chromosomal research, this information is captured by the so-called "HiC" analysis). A standard output is the following contact probability map between each pair of loci on a DNA segment.