Math Insight

Let's imagine the molecule jumps left or right a distance of 1 mm at the rate of $\lambda = 0.01$ per second. (We are intentionally using ridiculously large jumps at a ridiculously slow rate for this initial approximation.)

When a jump occurs, $X(t)$ is equally likely to jump left or jump right by 1 mm. If we measure distance in millimeters, this process is equivalent to a birth-death process with birth and death rate both equal to 0.005.

In a small time window of length $\Delta t$, the probability that $X(t)$ increases from $x$ to $x+1$ is
$P(X(t + \Delta t)=x+1 \,|\, X(t) = x) =$

and the probability that it decreases from $x$ to $x-1$ is
$P(X(t + \Delta t)=x-1 \,|\, X(t) = x) = $
.

Since we only allow jumping up or down by one, the probability that it stays the same is
$P(X(t + \Delta t)=x \,|\, X(t) = x) = $
, and the probability of jumping to any other location $y \ne x-1, x,$ or $x+1$ is
$P(X(t + \Delta t)=y \,|\, X(t) = x) = $
.

Let's imagine that the molecule cannot move to a negative position. (Perhaps $X=0$ represents a cell membrane that the molecule cannot cross.) Instead, when $X(t)=0$ at the time of a jump, the molecule always jumps to the right. Therefore, the probability of a rightward jump in a time interval $\Delta t$ is
$P(X(t + \Delta t)=1 \,|\, X(t) = 0) =$

and the probability of no jump is
$P(X(t + \Delta t)=0\,|\, X(t) = 0) = $
.

Again, the probability of any other jump is zero, so for $y > 1$
$P(X(t + \Delta t)=y \,|\, X(t) = 0) = $
.

The resulting random left and right motion of the molecule is a random walk that is approximates diffusion. Using the R script diffusion_from_birth_death.R, you can simulate this diffusion process. First set n_samples to 1 and X_max to 10 so that the diffusion will keep going until the molecule moves 10 units to the right.

The script plots the resulting position versus time, except that it flips the axes so that the position $X$ is shown on the $x$-axis and the time $t$ is shown on the $y$-axis. Notice how, if the molecule hits the left “wall” at $X=0$, it must move rightward on the next jump. If the molecule hits the right side at $X=10$, the random walk ends and the script reports the time (in seconds) it took the molecule to diffuse from $X=0$ to $X=10$ mm.

Repeat the simulation to observe the large variation in these random walks. You can increase n_samples to 5 or 10 to see multiple walks plotted together.

Increase n_samples to 100, and the plot becomes too messy to convey much useful information, though one can visualize all the times when the walk hits 10 at the right. You can set plot_paths=FALSE to turn off plotting the random walks and instead set plot_histogram=TRUE to turn on a histogram of the times at which the molecule hits $X=10$.

The code also displays the average time for the molecule to diffuse to 10. If you increase n_samples as large as 10000, this sample average will be close to the actual average (or expected value) of this crossing time. The expected crossing time is actually a round number. You won't get it exactly from your samples, but you should be able to guess how large it is. The expected time for the molecule to diffuse from $X=0$ until $X=10$ is
seconds.

Now, let's make the random walk really simple and set X_max to 1. Before you run the code, can you guess what the average time should be for the molecule to take one jump from $X=0$ to $X=1$, given that the jump to the right occurs at a rate of $\lambda=0.01$? The expected time for the molecule to jump (not really diffuse, in this case) from $X=0$ to $X=1$ is
seconds. Verify that you get close to this average time if you make n_sample large.

The change we'll make is set X_max to 2, so that the molecule has to make two jumps to reach the right boundary. Before you do any simulations, do you have any idea what the average time would be for the molecule to reach the end? What's your guess?
seconds.

Now simulate the random walk for a large number of samples to get a good idea of the time to reach 2 (it's a round number). The expected time to reach 2 is
seconds.

When you doubled the distance, the expected time for the molecule to travel that distance was multiplied by
. Let's see if that trend continues. Set X_max to 4 and estimate the expected time for the molecule to reach 4 mm. It is
seconds. When you doubled X_max from 2 to 4, the expected travel time was multiplied by
. This means, when you quadrupled X_max from 1 to 4, the expected travel time was multiplied by
.

When X_max=7, what is the expected travel time?
seconds. Compared to the travel time when X_max was 1, the expected travel time was multiplied by
.

In general, how long will it take the molecule to travel a distance of $x$ mm?
seconds.

The fact that the travel time depends quadratically on the distance $x$ is a fundamental property of diffusive processes. If you remember just one thing about diffusion, remember this fact: if $T$ is the travel time and $x$ is the distance, then $T \propto x^2$ (where $\propto$ means "is proportional to").

You can experiment with different distances X_max to convince yourself that this relationship is true.

The fact that travel time increases quadratically with distance means the travel time across long distance will be much shorter slightly longer slightly shorter much longer
than travel time across short distances. Hence, using diffusion for the transport of molecules across small distances could be reasonable but replying on diffusion for longer distance transport would be very inefficient would be efficient would be reasonable would be slightly inefficient
.

Above, the waiting time until the next event was an
random variable with an expected value of
seconds. When an event occurred, the position of the molecule jumped to the right with probability
and jumped to the left with probability
(except for the case when the molecule's position was $X=0$, in which case the molecule jumped to the right with probability
).

Now, if you change X_max to 1, you should find that there is no variability and that the time for the molecule to reach X_max is exactly
seconds. However, if you increase X_max to 10, you won't notice a lot of difference (other than that your code runs faster as it doesn't have to draw exponential random variables).

Just like the previous case, the expected time for the molecule to diffuse to a distance $x$ will be $T=$
seconds.

To explore the effects of decreasing the jump size, download the R script discrete_diffusion_across_interval.R. This file is similar to the script you have been using so far, except with some minor changes to make it better handle jump sizes that aren't integers (and the round off error that can occur when computers compute with numbers that include decimals). We'll denote the jump size by $\Delta x$ (which is written as Delta_x in the script) and the time between jumps by $\Delta t$ (which is written as Delta_t in the script).

Set $X_{max}$ (written as X_max in the script) to 10, the jump size $\Delta x$ to 1 and the time $\Delta t$ between jumps to 100. This parameter set should be identical to the first set of parameters used above. The expected time to reach 10 mm is the same round number of
seconds. (Remember, you can increase the variable n_samples to increase the number of samples paths of the molecule to average over.)

Since we don't want the molecule to stop moving as we let the jump size approach zero, reducing $\Delta t$ to keep it proportional to $\Delta x$ is the wrong idea. To make the molecule move faster, we have to decrease $\Delta t$ faster than $\Delta x$ (so the rate of the jumps increases faster than the jump size).

The clue for what to do with $\Delta t$ can be seen in our earlier experiments with increasing $X_{max}$ while keeping $\Delta x$ fixed at one. Cutting $\Delta x$ in half is equivalent increasing $X_{max}$ by a factor of
, in that both manipulations triple cut in half keep the same double
the number of rightward jumps needed to reach $X_{max}$. Above, when we doubled $X_{max}$, the time for the molecule to diffuse to $X_{max}$ increased by a factor of
. Therefore, if we keep $X_{max}$ and $\Delta t$ the same, but cut $\Delta x$ in half, the time to diffuse to $X_{max}$ should also increase by a factor of
. If we want to keep the diffusion time the same, we should decrease the time between jumps $\Delta t$ by a factor of
.

If we decrease $\Delta x$ by factor of two (cut it in half) and decrease $\Delta t$ by a factor of four, what happens to the ratio $\frac{(\Delta x)^2}{\Delta t}$? The ratio decreases stays the same increases
. If we keep the ratio $\frac{(\Delta x)^2}{\Delta t}$ constant as we decrease $\Delta x$ toward zero, what should happen to the time to diffuse to $X_{max}$? It approximately stays the same decreases increases
.

The quantity $$D = \frac{ (\Delta x)^2}{2\Delta t}$$ is called the diffusion coefficient. (Don't worry about the extra factor of 1/2. This version of $D$ is actually a more natural choice than the ratio $(\Delta x)^2/\Delta t$ for reasons we won't discuss.) It captures how fast the molecule is diffusing, independent of our choice of $\Delta x$ or $\Delta t$.

Experiment with the simulation to convince yourself that $D$ is the quantity that determines how long it takes the molecule to diffuse to $X_{max}$. Let's use the diffusion coefficient determined from our original simulations, when $\Delta x=1$ mm and $\Delta t = 100$ seconds. With this choice $D = $
$\frac{mm^2}{s}$. (The strange units of square millimeters per second just turns out to be what we need for a diffusion coefficient.)

If you create a variable for $D$ with the line

D = ＿

Delta_t = Delta_x^2/(2*D)

Set $X_{max} $ to 1 mm and use this value of $D$. Estimate how long it takes the molecule to diffuse to $X_{max}$ using different jump sizes of $\Delta x$ between 1 and 0.01. (For smaller values of $\Delta x$, you may need to decrease n_samples so that the code doesn't take too long.) The mean transport time should be a round number:
seconds. It shouldn't depend on the value of $\Delta x$ (although unless you run it for a huge number of samples, the reported mean value will vary a lot between runs).

Estimate how long it would take the molecule to diffuse 0.01 mm.
seconds.