© 2021 Rebecca Rousseau. This work is licensed under a Creative Commons Attribution License CC-BY 4.0. All code contained herein is licensed under an MIT license.
Objective: Last class, we used a random number generator to plot stochastic trajectories of diffusive particles and get a sense of the evolution of the probability distribution for the random walk position over time. Today, we will arrive at this distribution more directly by using numerical integration of a master equation for diffusion. We will then apply this approach to simulate the FRAP experiment and analyze the evolution of the mRNA distribution for a constitutive promoter.
Assuming one-dimensional motion along the x-axis, we want to determine the probability of a given diffusive particle being at position $x$ at a given timepoint. This diffusive particle can take at most one step in time $\Delta t$ of size $a$ in either direction along the x axis. There are then three possible events that could occur in time $\Delta t$ to arrive at position $x$: (1) The given particle was at position x at time t and remained there, (2) the particle jumped from x-a to x with a jump probability $k\Delta t$, or (3) the particle jumped from x+a to x with jump probability $k\Delta t$. Therefore,
$$ P(x,t+\Delta t) = (1-2k\Delta t)P(x,t) + (k\Delta t)P(x-a,t) + (k\Delta t)P(x+a, t)\\ = P(x,t)+ (k\Delta t)P(x-a,t) + (k\Delta t)P(x+a, t)-2k\Delta t P(x,t)$$The master equation at the boundaries of the random walk are slightly modified from the above as
$$ P(0, t+\Delta t) = P(0,t) + (k\Delta t)P(a, t)-k\Delta t P(0,t)\\ P(N, t+\Delta t) = P(N,t) + (k\Delta t)P(N-a, t)-k\Delta t P(N,t).$$For our first numerical integration, we will model how an initial point source of particles diffuses over time in 1D. To do this, we will set up an array of boxes through which our particles will diffuse, and specify the box in the center to initially have probability 1 and all other boxes to have probability 0 (i.e., all particles are initially found in the central box).
# Import modules
import matplotlib.pyplot as plt
import numpy as np
import panel as pn
pn.extension()
import scipy.stats as stats
import seaborn as sns
sns.set()
%config InlineBackend.figure_format = 'retina'
Now, let's set up our array of boxes and our initial conditions.
# number of boxes
n_boxes = 40
# number of steps
n_steps = 100
# time increment in seconds (make sure time increment < hopping timescale for accuracy)
dt = 0.1
# hopping rate constant (k) in sec^-1
k = 1
# Initialize the probability distribution
prob = np.zeros([n_boxes, n_steps]) # each column represents the distribution across boxes at the timeslice
# Find the middle box
n_center = int(n_boxes / 2)
# Set middle box to have probability 1 at the initial time step
prob[n_center, 0] = 1
# Loop through time points and generate the probability distribution across all positions
# Loop through time points
for t in np.arange(1, n_steps):
# Loop through positions
# FIRST: loop probabilities at nonzero time for NON-boundary states
for n in np.arange(1, n_boxes-1):
prob[n, t] = prob[n, t-1] + k * dt * prob[n-1, t-1] + k * dt * prob[n+1, t-1] - 2 * k * dt * prob[n, t-1]
# SECOND: handle the boundaries at time step t
prob[0, t] = prob[0, t-1] + k * dt * prob[1, t-1] - k * dt * prob[0, t-1]
prob[n_boxes - 1, t] = prob[n_boxes - 1, t-1] + k * dt * prob[n_boxes - 2, t-1] - k * dt * prob[n_boxes - 1, t-1]
# Interactive plot of the probability distribution over time
frame_player = pn.widgets.Player(
name='movie player',
start=0,
end=n_steps,
step=1,
interval=50,
)
@pn.depends(frame_player.param.value) # make a new plot when the parameters in frame_player change
def movie_player(frame):
# Initialize plot
fig, ax = plt.subplots(1,1) # Making one figure in our interactive plot
# Generate bar plot of distribution
ax.bar(np.arange(n_boxes), prob[:, frame])
ax.set_xlabel('position')
# Set height for all plots
ax.set_ylim([0, 1])
plt.close(fig)
return fig
pn.Row(movie_player, frame_player)
We see here, as in the coin-flip exercise, that the probability distribution broadens over time while remaining centered at the origin. In this simulation, however, we derived the probability distribution explicitly rather than stochastically (i.e., using coin flips).
We now apply this approach of numerically integrating the chemical master equation to the problem of fluorescence recovery after photobleaching, aka FRAP.
In FRAP experiments, the fluorescently labeled molecules initially have a uniform distribution. Applying high intensity light to a specific region causes the fluorescent molecules there to be photobleached so they no longer fluoresce, but over time, fluorescent molecules that were outside the photobleached region will diffuse into it, ultimately recovering the uniform fluorescence distribution. By estimating the timescale of fluorescence recovery, experimentalists can then determine the diffusion constant of fluorescently labeled molecules.
We will now simulate the 1D version of the FRAP experiment, initializing the probability distribution as zero in the middle region and uniform outside it to indicate that the central region starts out photobleached. We will then observe the process of recovery as the gap gets "filled" with probability from diffusion into the photobleached region.
# Number of steps
n_steps = 1000
# Initialize probability distribution
p = np.zeros([n_boxes, n_steps])
# Assign nonzero values to the first time step
p[:,0] = 1 # Uniform, will normalize after photobleaching effect introduced
# Number of boxes to photobleach
n_bleach = 10
# Find start and end positions for the photobleached region
start = n_center - int(n_bleach/2)
end = n_center + int(n_bleach/2)
# "Do" the photobleach by setting the probability to 0 in that region
p[start:end, 0] = 0
# Renormalize so probability sums to 1
p[:,0] = p[:,0]/np.sum(p[:,0])
We can now run the same code as before but for these new initial conditions.
# Loop through time points and generate the probability distribution across all positions
# Loop through time points
for t in np.arange(1, n_steps):
# Loop through positions
# FIRST: loop probabilities at nonzero time for NON-boundary states
for n in np.arange(1, n_boxes-1):
p[n, t] = p[n, t-1] + k * dt * p[n-1, t-1] + k * dt * p[n+1, t-1] - 2 * k * dt * p[n, t-1]
# SECOND: handle the boundaries at time step t
p[0, t] = p[0, t-1] + k * dt * p[1, t-1] - k * dt * p[0, t-1]
p[n_boxes - 1, t] = p[n_boxes - 1, t-1] + k * dt * p[n_boxes - 2, t-1] - k * dt * p[n_boxes - 1, t-1]
# Interactive plot of the probability distribution over time
frame_player = pn.widgets.Player(
name='movie player',
start=0,
end=n_steps,
step=1,
interval=50,
)
@pn.depends(frame_player.param.value)
def movie_player(frame):
# Initialize plot
fig, ax = plt.subplots(1,1)
# Generate bar plot of distribution
ax.bar(np.arange(n_boxes), p[:, frame])
ax.set_xlabel('position')
# Set height for all plots
ax.set_ylim([0, 0.05])
plt.close(fig)
return fig
pn.Row(movie_player, frame_player)
So far, we have considered unbiased systems of random motion where the hopping rate between neighboring states is the same regardless of direction and mass is conserved. When considering the production of mRNAs from gene transcription, however, forward and reverse rates are not equal. Specifically, we must introduce the notion of production and degradation rates, with production generally favored. Letting $r$ be the production rate and $\gamma$ be the degradation rate, the adapted master equation is
$$ P(n, t+\Delta t) = P(n, t) + r\Delta t P(n-1, t) + \gamma(n+1)\Delta t P(n+1, t) - r\Delta t P(n, t) - \gamma n\Delta t P(n, t),$$where the probability of degradation depends on rate $\gamma$ and the number of mRNA transcripts present (i.e., available to degrade).
# Define the parameters for mRNA production and degradation
# production rate
r = 1
# degradation rate
gamma = 1 / 10
# Define number of time steps for integration
n_steps = 500
# Define dt
dt = 0.1
# Define the number of mRNA copies to keep track of
n_mRNA = 30
# Initialize probability distribution
P = np.zeros([n_mRNA, n_steps])
# Set initial condition, starting at the beginning of transcription event
P[0, 0] = 1
# Loop through time points
for t in range(1, n_steps):
# Loop through mRNA copy number
for n in range(1, n_mRNA - 1):
P[n, t] = P[n, t-1] + r * dt * P[n-1, t-1] + gamma * (n+1) * dt * P[n+1, t-1] - r * dt * P[n, t-1] - \
gamma * n * dt * P[n, t-1]
# Boundary states
P[0, t] = P[0, t-1] + gamma * dt * P[1, t-1] - r * dt * P[0, t-1]
P[-1, t] = P[-1, t-1] + r * dt * P[-2, t-1] - r * dt * P[-1, t-1] - gamma * n_mRNA * dt * P[-1, t-1]
frame_player = pn.widgets.Player(
name='movie player',
start=0,
end=n_steps,
step=1,
interval=50,
)
@pn.depends(frame_player.param.value)
def movie_player(frame):
# Initialize plot
fig, ax = plt.subplots(1, 1)
# Generate bar plot of distribution
ax.bar(range(n_mRNA), P[:, frame])
# Generate Poisson distribution
poisson_distr = stats.poisson.pmf(range(n_mRNA), r / gamma) # mean mu = r/gamma
# Plot steady state distribution
ax.plot(range(n_mRNA), poisson_distr)
ax.set_xlabel('mRNA / cell')
# Set y axis limits
ax.set_ylim([0, 1])
plt.close(fig)
return fig
pn.Row(
movie_player, frame_player
)
As we approach the steady state, the mean tends to scale with the variance of the distribution, fitting a Poisson curve.