In 1-D, the concentration C of a particular chemical species can be expressed as a function of 2 variables, the linear position and the time:
\(C(x, t)\)
The method BioSim1D.diffuse_step(time_step) performs a single diffusion step, for a single chemical species.
The position is discretized over n_bins points
The rate of change with time of the concentration function is given by:
\({\partial \over \partial t} C(x, t) = D {\partial^2 \over \partial x^2} C(x,t)\)
provided that the Diffusion coefficient D is a constant. That's often assumed, but not always the case; for example, an argument can be made that it has some dependency on the concentration (see for example this paper.) In a later phase, we'll explore modeling variations in the diffusion coefficient, but for now we'll assume it to be constant.
Plain-English intuition:
If the concentration is the same everywhere in space, then the right-hand side becomes zero... and there's no change in time: the system has equilibrated.
If the concentration has a constant gradient throughout space (imagine the slope of a sand dune), then the first derivative (rate of change) with respect to position is constant... and therefore the second derivative (change of the change) is zero. Again, no change in time. Imagine the sand dune in equilibrium: at points on the dune, sand falls in from above - and falls out to below at the same rate. (At the boundaries, special intervention must be applied - such as constant infusion at the high end and constant drain at the low end.)
So, all said and done, the equation says that rate of change with time depends on how changeable the gradient throughout space is.
Want to follow the math, but not feeling in your element?
Check out these
brilliant short videos, with clear explanations and immensely helpful animations: a YouTube channel
with 4 million subscribers, well deserved!We'll use the following 2 boundary conditions:
1. We specify the initial concentrations at each position, at time 0 (essentially, an "initial condition"):
\(C(x_i, 0) = C_0(x_i)\)
2. No diffusion (flux) at the end bins, into the outside. The container is considered "sealed." A class variable BioSim1D.sealed is by default True, with the expectation that later similations will explore scenarios where the container is not sealed. In the hypothetical space outside the boundary of the system, the concentrations remain fixed to the value of the edge bin.
If xL and xR are, respectively, the leftmost and righmost coordinates of our system, we'll require:
\(C(x, t) = C(x_L, t)\) for any x < xL and any t
and
\(C(x, t) = C(x_R, t)\) for any x > xR and any t
We'll discretize C(x, t) as \(C(x_i, t_\mu)\) , where xi and tμ are discrete values for position and time, respectively.
The current simulation, as done by the class BioSim1D, makes use of a set of "bins" (discrete x-positions), numbered in the range from 0 to (n_bins-1), inclusive.
We use the term "bin" rather than "cell" to avoid confusion with modeling of biological cells.
And we're using a Greek letter as an index to discrete time, to avoid confusion with the spacial index. As a mnemonic, tμ , "t sub mu" spells "time".
The following is the basis for the computations in the method BioSim1D.diffuse_step
The most straightforward approach is to use multivariate finite differences to approximate the partial derivatives.
In the 1-D diffusion equation, we'll use the forward difference to approximate the time derivative, and the central difference to approximate the 2nd spacial derivative (this approach is known as the "3+1 stencil" or "Explicit Forward-Time Centered Space"):
\({\partial \over \partial t} C(x, t) \approx { C(x \ ,\ t + \Delta t) - C(x \ , \ t ) \over {\Delta t} }\)
and
\({\partial^2 \over \partial x^2} C(x,t) \approx {C(x+\Delta x\ , \ t) - 2 \ C(x,\ t) + C(x-\Delta x\ ,\ t) \over (\Delta x)^2}\)
Substituting the above in the 1-D diffusion equation \({\partial \over \partial t} C(x, t) = D {\partial^2 \over \partial x^2} C(x,t)\) , and switching to discretized values:
\({ C(x_i \ ,\ t_{\mu+1}) - C(x_i \ , \ t_\mu ) \over {\Delta t} } = D [ {C(x_{i+1} , \ t_\mu) - 2 \ C(x_i,\ t_\mu) + C(x_{i-1} ,\ t_\mu) \over (\Delta x)^2}]\)
i.e.
Given the knowledge of the values of the concentration C at time tμ , at the 3 positions xi-1 , xi and xi+1 , the above equation lets us compute the value of C at the next time step tμ+1 , at the position xi
The Software: the above equation is implemented in the diffuse_step_single_species method of the class BioSim1D.
Reposity file.
(Future releases will explore more efficient implementations, including parallelized computations.)
The Intuition Behind the Equation : advance 1 step in time
For any given bin, the above equation permits the computation of the concentration of a particular chemical species at time tμ+1 if we know it at time tμ
The new concentration is the old value, plus a corrective term (the one starting with D)
Why does the concentration change? Well, because of an influx in or out from the neighboring bins (xi-1 to the left and xi+1 to the right).
The part in square brackets can be regard as the sum of 2 terms; in simplified notation: [C(right) - C] + [C(left) - C] The 1st one is the influx from the right, while the 2nd one is the influx from the left. That's where the -2C arises from. The "influxes" depend on the concentration imbalance, and of course may be negative if the original bin has higher concentrations that its neighbors.
If you imagine the concentration values as money, the middle bin receives (possibly negative) deposits from its left and right neighbors; adding those up will give the new "bank balance."
The D factor is the Diffusion rate; the higher, the more change per unit time.
Likewise, if Δt is larger, i.e. we consider a longer time span, more change will occur.
The (Δx)2 in the denominator may be the least intuitive part... What happens if the bins get smaller, i.e. a smaller Δx ? That means they're packed more closely. For starters, the concentration gradient become steeper (more change per unit distance.) On top of that, the diffusing molecules have less distance to travel to spread across to the next bin. Hence, the more rapid change - on both counts - of the concentration values. Close-by regions simply mix up much better!
Boundary conditions:
All the concentration values at t=0 are assumed known: that's our initial setup. From that, we can compute the values at all later times, step by step.
But what do we do when i=0 or i=MAXi , i.e. when we're at the edges of the system, and there is no xi-1 nor xi+1 ?
Conceptually, we can imagine an extra bin to the left of the leftmost bin, with the same concentration; and, likewise, an extra bin to the right of the righmost one, again with the same concentration. Because of the identical concentrations, there's be no flux of the chemical. That is, nothing "leaking" in or out of the container, as specified by our 2nd boundary condition.
Therefore, when i=0, the term \([ C(x_{i+1} , \ t_\mu) - 2 \ C(x_i,\ t_\mu) + C(x_{i-1} ,\ t_\mu)]\) becomes \([ C(x_1 , \ t_\mu) - 2 \ C(x_0,\ t_\mu) + C(x_0 ,\ t_\mu)]\) , i.e. \([ C(x_1 , \ t_\mu) - \ C(x_0,\ t_\mu) ]\)
Similarly, when i=MAXi, the term \([ C(x_{i+1} , \ t_\mu) - 2 \ C(x_i,\ t_\mu) + C(x_{i-1} ,\ t_\mu)]\) becomes \([ C(x_{MAXi} , \ t_\mu) - 2 \ C(x_{MAXi},\ t_\mu) + C(x_{{MAXi}-1} ,\ t_\mu) ]\) , i.e. \([ - \ C(x_{MAXi},\ t_\mu) + C(x_{{MAXi}-1} ,\ t_\mu) ]\)
Stability:
The above approach can lead to instability in the computed results if Δt is too large relative to Δx. The Von Neumann stability analysis shows possible instability when:
\(\Delta t > {(\Delta x)^2 \over 2 D}\)
For reasons discussed in the experiment 1D/diffusion/overly_large_single_timestep, we'll be using a slightly more strict ceiling on the value of Δt.
In multiple dimentions, the concentration C of a particular chemical species can be expressed as a function of 2 variables; but, this time, the linear position of the 1-dimensional case gets replaced with the position vector r :
\(C(\mathbf{r} , t)\)
In 2-D, r = (x, y).
The value of the function C is a scalar (a single number), just like in 1-D. To keep in mind that, as usual, C is for a particular chemical species. So, there's really a family of functions Ci , for each of the chemical species. We're just considering one of them.
The rate of change with time of the concentration function is given by:
\({\partial \over \partial t} C(\mathbf{r}, t) = D \ \nabla ^2C(\mathbf{r},t)\)
again, assuming for now that the Diffusion coefficient D is a constant. ∇2 is the Laplace operator.
In 2-D, using a less terse notation:
\({\partial \over \partial t} C(\mathbf{r}, t) = D [ {\partial^2 \over \partial x^2}(\mathbf{r},t) + {\partial^2 \over \partial y^2}(\mathbf{r},t) ]\)
Note that the right-hand side, just like the left-hand side, is a scalar function of x, y and t; for any given values of x, y and t, it will evaluate to a number.
We'll discretize C(r, t) as \(C(x_i, y_j, t_\mu)\) , where xi , yj and tμ are discrete values for position and time.
The current simulation, handled by the class BioSim2D, makes use of a set of "bins" (discrete x-positions and y-positions), numbered in the range from 0 to (n_bins_x - 1) and (n_bins_y - 1), inclusive.
We use the term "bin" rather than "cell" to avoid confusion with modeling of biological cells.
And we're using a Greek letter as an index for the discrete time, to avoid confusion with the spacial indices. As a mnemonic, tμ , "t sub mu" spells "time".
The following is the basis for the computations to be released at a later date in an early version of the method BioSim2D.diffuse_step
As mentioned before, Life123 will later move to existing open-source libraries.
The most straightforward approach is to use multivariate finite differences to approximate the partial derivatives.
In the 2-D diffusion equation, we'll use the forward difference to approximate the time derivative, and the central difference to approximate the 2nd spacial derivatives:
\({\partial \over \partial t} C(x,y, t) \approx { C(x \ , \ y,\ t + \Delta t) - C(x, \ y \ , \ t ) \over {\Delta t} }\)
and
\({\partial^2 \over \partial x^2} C(x,y, t) \approx {C(x+\Delta x\ , \ y, \ t) - 2 \ C(x,\ y,\ t) + C(x-\Delta x\ ,\ y, \ t) \over (\Delta x)^2}\)
\({\partial^2 \over \partial y^2} C(x,y, t) \approx {C(x, \ y+\Delta y\ , \ t) - 2 \ C(x, \ y,\ t) + C(x, \ y-\Delta y\ ,\ t) \over (\Delta y)^2}\)
Substituting the above in the 2-D diffusion equation, and switching to discretized values:
\({ C(x_i \ ,\ y_j \ , \ t_{\mu+1}) - \ C(x_i \ \ ,\ y_j \ , \ t_\mu ) \over {\Delta t} } = D \left[ {C(x_{i+1} ,\ y_j \ , \ t_\mu) - 2 \ C(x_i ,\ y_j \ , \ t_\mu) + C(x_{i-1} ,\ y_j \ , \ t_\mu) \over (\Delta x)^2} + {C(\ x_i \ , \ y_{j+1} , \ t_\mu) - 2 \ C(x_i ,\ y_j \ , \ t_\mu) + C( x_i \ , \ y_{j-1} ,\ t_\mu) \over (\Delta y)^2} \right]\)
i.e.
If Δx = Δy (which we'll call Δs, where s stands for "spacial"), the above simplifies to:
\(\begin{align} C(x_i \ , \ y_j \ , \ t_{\mu+1}) =\ &C(x_i \ , \ y_j \ , \ t_{\mu}) \\ + D { \ \Delta t \over (\Delta s)^2 } & \ [ C(x_{i+1} ,\ y_j \ , \ t_\mu) + C(x_{i-1} ,\ y_j \ , \ t_\mu) \\ &+ C(x_i \ , \ y_{j+1} , \ t_\mu) + C( x_i \ , \ y_{j-1} ,\ t_\mu) - 4 \ C(x_i ,\ y_j \ , \ t_\mu) ] \end{align}\)
Given the knowledge of the values of the concentration C at time tμ , at the positions (xi, yj) , (xi-1, yj) , (xi+1, yj), (xi, yj-1), (xi, yj+1) , the above equation lets us compute the value of C at the next time step tμ+1 , at the position (xi , yj)
Note: the above is a coarse approximation that weighs in the concentration at (xi, yj), and at its neighbors to the left, right, top and bottom - but not at any diagonal positions. This is called a 5-point stencil.
The Software: the above equation is implemented in the diffuse_step_single_species method of the class BioSim2D. (Future releases will explore more accurate and efficient implementations; in particular, switching to existing open-source libraries.)
The last formula, inside the large square brackets, combines 5 terms: -4 times the concentration at the point of interest (xi, yi), plus the concentrations of its neighbors to the left, right, top and bottom.
That can be expressed more concisely as a convolution with the matrix:
\(\mathbf{L}_{xy}^2 = \begin{pmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1& 0 \end{pmatrix}\)
where L stands for the Discrete Laplacian.
The above notation makes is easier to see the 5-point stencil.
We're assuming for now that Δx = Δy, and we're calling that quantity Δs.
A finer approximation will use the diagonal elements, too, i.e. all 9 values - and hence called a 9-point stencil. Its coefficients are:
\(\mathbf{L}_{xy}^2 = {1 \over 6} \begin{pmatrix} 1 & 4 & 1 \\ 4 & -20 & 4 \\ 1 & 4 & 1 \end{pmatrix}\)
The above approximation is more accurate, though still of second order with respect to Δs.