Math Insight

Introduction to diffusion

Math 2241, Spring 2023
ID #:
Due date:
Table/group #:
Group members:
  1. A molecule is moving randomly left and right. Let $X(t)$ be its position at time $t$. For now, imagine that $X(t)$ jumps left and right in a manner similar to a birth-death process.
    1. 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) = $

    2. 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) = $

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

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

    5. 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?

      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

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

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

      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.

    8. The fact that travel time increases quadratically with distance means the travel time across long distance will be
      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

  2. We introduced the diffusion process by drawing an analogy to a birth-death process. It turns out that we can obtain similar behavior even if we simplify the process.
    1. 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

    2. This process had two sources of randomness, the random waiting time and the random jump direction. We can simplify the process by eliminating the first source of randomness and just making the waiting time be fixed to the value of 100 seconds.

      To make this change, eliminate the variable lambda in the R script and replace it with the variable dt=100. Then, find the line where the time is increased by an exponential random variable,

      t = t + rexp(1,rate=lambda)
      and replace it with the simpler statement that the time increases by the value dt.
      t = t + dt

      If you like, you can rename the script to something ilke discrete_diffusion.R, as it is no longer a diffusion based on a birth-death process.

    3. 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=$

  3. (Optional bonus question) So far, we've been making the molecule take ridiculously large jumps of 1 mm at ridiculously large time intervals of 100 seconds. In reality, a molecule that is diffusing will move only a tiny amount at a time before changing directions (say by colliding with another molecule).
    1. 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.)

    2. Now, we can start changing the jump size $\Delta x$ to smaller values in an attempt to make the diffusion more representative of a molecule's behavior (though, of course, we are still allowing the molecule to move in only one dimension, left and right along the interval from 0 to 10 mm).

      As we decrease the jump size $\Delta x$, we should simultaneously make the time in between jumps $\Delta t$ be
      if we want the time to diffuse to $X_{max}$ roughly the same value.

      To begin, cut $\Delta x$ in half to 0.5 mm. One might expect that a reasonable option would be to simultaneously cut the time interval $\Delta t$ in half to
      . When you do that, what happens to the time to reach $X_{max}=10$ mm? It approximately

      Continue shrinking the jump size $\Delta x$ to 0.25 mm and then to 0.1 mm, at the same time shrinking time between jumps $\Delta t$ proportionally to 25 seconds and then to 10 seconds. What happens to the time to reach $X_{max}=10$ mm?
      In other words, as the jump size decreases, the molecule is moving

      When we decrease the jump size $\Delta x$ and the time between jumps $\Delta t$ by the same factor, what happens to the ratio $\frac{\Delta x}{\Delta t}$? The ratio
      . Since we want to represent the diffusion of the molecule in a manner that doesn't depend on the actual jump size and we want to work with continuous space (rather than just intervals of width $\Delta x$), we are interested in moving to the limit where $\Delta x$ approaches zero. If we let $\Delta x$ approach zero and, at the same time, keep the ratio $\frac{\Delta x}{\Delta t}$ constant, what do you think would happen to the time for the molecule to diffuse to $X_{max}$? The time would approach
      . What happens to the speed of molecules diffusion? Its speed approaches

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

    4. 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 = _
      then you can set $\Delta t$ to the correct value depending on the value of $\Delta x$ with the line
      Delta_t = Delta_x^2/(2*D)
    5. 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).

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