Life123 : Open-Source Engine for Quantitative Reactomics

Dynamical Modeling of Biological Systems ♦ In-Silico Experiments

Dynamical Modeling of Biological Systems ♦ In-Silico Experiments

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!

“But what is a partial differential equation?”

“Solving the heat equation” (note: same form as the diffusion equation)

“But what is a partial differential equation?”

“Solving the heat equation” (note: same form as the diffusion equation)

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 x_{L} and x_{R} are, respectively, the leftmost and righmost coordinates of our system, we'll require:

\(C(x, t) = C(x_L, t)\) for any x < x_{L} and any t

and

\(C(x, t) = C(x_R, t)\) for any x > x_{R} and any t

An alternate second boundary condition that we might explore later on, is to specify a specific fixed concentration (possibly, slowly varying in time) outside of our system; i.e. imagine it "immersed in a bath" of fixed concentration. The Diffusion coefficient for flux into/out of the "container" might be different than the value in the bulk medium; i.e. semi-permeable membranes at the boundary.

We'll discretize C(x, t) as \(C(x_i, t_\mu)\) , where x_{i }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 **m**u" spells "time".

IMPORTANT: the following is for general background. Life123 will soon switch to existing libraries that also handle variable meshes (bin sizes.)

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 2^{nd} 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.

\(\begin{align} C(x_i \ ,\ t_{\mu+1}) =\ &C(x_i \ , \ t_\mu ) \\ &+ D \ {\Delta t \over (\Delta x)^2 } \ [ C(x_{i+1} , \ t_\mu) - 2 \ C(x_i,\ t_\mu) + C(x_{i-1} ,\ t_\mu)] \end{align}\)

Given the knowledge of the values of the concentration **C** at time t_{μ} , at the 3 positions x_{i-1} , x_{i} and x_{i+1} , the above equation lets us compute the value of C at the next time step t_{μ+1} , at the position x_{i}

**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 (x_{i-1} to the left and x_{i+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=MAX_{i} , i.e. when we're at the edges of the system, and there is no x_{i-1} nor x_{i+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=MAX _{i}**, 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 C_{i} , 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 x_{i }, y_{j }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 **m**u" spells "**t**i**m**e".

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 2^{nd} 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.

\(\begin{align} C(x_i \ , \ y_j \ , \ t_{\mu+1}) =\ &C(x_i \ , \ y_j \ , \ t_{\mu}) \\ + D \ \Delta t & \ [ {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} ] \end{align}\)

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 (x_{i}, y_{j}) , (x_{i-1}, y_{j}) , (x_{i+1}, y_{j}), (x_{i}, y_{j-1}), (x_{i}, y_{j+1}) , the above equation lets us compute the value of C at the next time step t_{μ+1} , at the position (x_{i} , y_{j})

Note: the above is a **coarse approximation** that weighs in the concentration at (x_{i}, y_{j}), 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 (x_{i}, y_{i}), 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.