© 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.

Simulating random walk from the chemical master equation

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).

Now, let's set up our array of boxes and our initial conditions.

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.

Application: 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.

We can now run the same code as before but for these new initial conditions.

Application: mRNA distribution over time for constitutive gene

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).

As we approach the steady state, the mean tends to scale with the variance of the distribution, fitting a Poisson curve.