May 21, 2018, Monday

# Approximating Diffusion as a "Flux"

## Introduction

We are interested in approximating diffusion as a "flux" for the easy analysis of the biological system. Currently, using systems biology techniques, we are only able to apply fundamental subspace analysis to chemical reactions of a particular biological system, and not to the inevitable diffusive aspects. Mathematically, we have a biological model that follows the form:

${\frac {\partial \phi }{\partial t}}=S{\textbf {v}}_{{1}}+D\phi$,

where S is the stoichiometric matrix, D is the diffusion matrix, ${\textbf {v}}_{{1}}$ is the vector of reaction fluxes, and $\phi$ is the vector of concentrations. We wish to derive a matrix (S:D) that includes information on both the stoichiometry and the diffusion of the system by approximating diffusion as a "flux", such that:

${\frac {\partial \phi }{\partial t}}=(S:D){\textbf {v}}_{{2}}$.

(We denote the reaction-diffusion "flux" vector as ${\textbf {v}}_{{2}}$ to differentiate it from ${\textbf {v}}_{{1}}$.)

## Approximation Using the Second Derivative

Fick's laws of diffusion provide us with a general basis for the construction of a matrix that includes information on the diffusion of a particular biological system. Namely, we have:

Fick's First Law, in one dimension: $J=-D{\frac {\partial \phi }{\partial x}}$,

where $J$ is the diffusion flux, $D$ is the diffusion coefficient, and $x$ is the position. In other words, the flux of a diffusive system moves from regions of higher concentration to regions of lower concentration proportional to the concentration gradient, which is a function of space.

Fick's Second law, in one dimension: ${\frac {\partial \phi }{\partial t}}=D{\frac {\partial ^{{2}}\phi }{\partial x^{{2}}}}$

with the same definitions as above. This law defines how diffusion causes the concentration at a particular position to change with time.

Using Fick's Second Law in particular, we see that the general form of the partial differential equations in a reaction-diffusion model is as follows:

${\frac {\partial \phi }{\partial t}}=S{\textbf {v}}_{{1}}+D{\frac {\partial ^{{2}}\phi }{\partial x^{{2}}}}$.

Thus, the key to approximating diffusion in this system lies in the second derivative. Using the formal limit definition of the second derivative, where

$f''(x)=\lim _{{h\to 0}}{\frac {f(x+h)-2f(x)+f(x-h)}{h^{2}}}.$

after writing the second derivative as a difference quotient of two first derivatives, that is

${\frac {f(x+h)-2f(x)+f(x-h)}{h^{2}}}={\frac {{\frac {f(x+h)-f(x)}{h}}-{\frac {f(x)-f(x-h)}{h}}}{h}}$,

We see that it becomes logical to approximate diffusion by examining the behavior at a finite number of evenly spaced points, using this limit definition for the second derivative in accordance with Fick's Second Law. Specifically, to approximate the concentration change at a particular point $x_{{i}}$, we utilize information regarding the current concentration at three points: $x_{{i-1}},x_{{i}},x_{{i+1}}$. Of course, in general, the approximation improves as the distance between the discrete points shrinks to zero, by definition. Here, we will use three points in the spatial domain, because this number is appropriate for our C5a production example. An extension to a greater number of approximation points within the spatial domain may be found here.

Additionally, it may be helpful to view diffusion as effectively a system of chemical "reactions", where

$u_{{1}}\leftrightarrow u_{{2}}$
$u_{{2}}\leftrightarrow u_{{3}}$

and the behavior of $u_{{1}}$ and $u_{{2}}$ is determined by the boundary conditions at those points.

## Does This Approximation Work?

Figure 1: Diagram of the four points for the approximation of diffusion.

Again, we are approximating diffusion at three discrete, evenly spaced points. Again, by the nature of diffusion and our model, this necessarily requires a point to the left and to the right of each point at which we are interested in approximating diffusion at. Then we have five discrete, evenly spaced points at which we require concentration information to apply Fick's laws to our approximation: 0, 1, 2, 3, 4 (Figure 1). From our second derivative approximation in the previous section, we have, in general, the following system of differential equations:

${\frac {\partial \phi _{{1}}}{\partial t}}=D({\frac {\phi _{{0}}-2\phi _{{1}}+\phi _{{2}}}{h^{{2}}}})$
${\frac {\partial \phi _{{2}}}{\partial t}}=D({\frac {\phi _{{1}}-2\phi _{{2}}+\phi _{{3}}}{h^{{2}}}})$
${\frac {\partial \phi _{{3}}}{\partial t}}=D({\frac {\phi _{{2}}-2\phi _{{3}}+\phi _{{4}}}{h^{{2}}}})$

The question that remains is how to convert this information into stoichiometric matrix form and thereby conduct analysis via SVD decomposition and the like. We know the "reactants" that must be present in this simple diffusion matrix: the three "reactants" really correspond to a single reactant at our chosen three equally-spaced points within the spatial domain. We are concerned with how to define the "fluxes" in such a matrix.

It is logical to separate fluxes by diffusive component -- for example, ${\frac {D}{h^{{2}}}}(\phi _{{2}}-\phi _{{1}})$ and ${\frac {D}{h^{{2}}}}(\phi _{{0}}-\phi _{{1}})$ in the first differential equation above, corresponding to the separate diffusion gradients on the "left hand side" and the "right hand side" of the spatial domains relative to point 1. It is worth noting that we are not explicitly specifying whether there is flux "into" or "out of" the system in terms of diffusion by Fick's First Law in our stoichiometric matrix. That is, if $\phi _{{1}}>\phi _{{0}}$ (or $\phi _{{2}}$), a -1 should be entered into the stoichiometric matrix, but if $\phi _{{1}}<\phi _{{0}}$ (or $\phi _{{2}}$), a +1 should be entered into the stoichiometric matrix. Still, as these are effectively equal but opposite fluxes, we can simply modify them for whatever biological situation we are interested in.

We are now interested in examining the effectiveness of our approximation for various boundary condition types. For our boundary conditions to "work", this approximation must blur the line between point and domain, such that the "boundary" conditions are taken across an entire area. Since our spatial domain we are interested in is from points 1 through 3, our left and right boundary conditions are applied in the space from point 0 to point 1 and point 3 to point 4, respectively. Therefore, in context of our diffusion approximation, we must specify "boundary" conditions as differences between concentrations at the two points specifying the edges of the boundary (not the spatial domain). For example, from point 0 to point 1 in the diffusion approximation for point 1, we define:

• Dirichlet boundary condition: $\phi _{{0}}-\phi _{{1}}=0$, as we are setting the concentration at the boundary at point 1 to be constant -- stated in its typical form, $\phi _{{1}}=\alpha$ for some concentration $\alpha$. This condition effectively states that the system is "well-mixed" past the boundary at point 1.
• Neumann boundary condition: $D({\frac {\phi _{{0}}-\phi _{{1}}}{h^{{2}}}})=\beta$, for some diffusive flux $\beta$ as we are specifying that the flux through point 1 is constant. Note that there may be a flux from point 2 in to point 1. A biological case of this might be that there exists a solid boundary at point 1 such that there is no net flux coming from within the solid matter, nor is there any flux entering the solid matter, in which case $\beta =0$ and the Neumann boundary condition is represented identically to the Dirichlet boundary condition.

We will explore toy examples for different boundary condition scenarios:

1. Dirichlet-Dirichlet
2. Dirchlet-Neumann
3. Neumann-Neumann

by constructing diffusion-only matrices such that ${\frac {\partial \phi }{\partial t}}=D{\textbf {v}}$ and exploring their biological validity via SVD decomposition and interpretation.

Dirichlet-Dirichlet Boundary Conditions

If we define Dirchlet boundary conditions at both the left and right boundaries of our spatial domain, we have:

$\phi _{{0}}-\phi _{{1}}=0$ and
$\phi _{{4}}-\phi _{{3}}=0$.

We generate a diffusion matrix after having eliminated these fluxes at the boundaries. Note that we have three "reacting species" (really, the same reactant at three points in our spatial domain) and four fluxes, as two have been canceled by the Dirichlet boundary condition.

But first, note that there are two equal and opposite fluxes included in our model; namely:

$D({\frac {\phi _{{2}}-\phi _{{1}}}{h^{{2}}}})$ in ${\frac {\partial \phi _{{1}}}{\partial t}}=-D({\frac {\phi _{{1}}-\phi _{{2}}}{h^{{2}}}})$ in ${\frac {\partial \phi _{{2}}}{\partial t}}$.
$D({\frac {\phi _{{3}}-\phi _{{2}}}{h^{{2}}}})$ in ${\frac {\partial \phi _{{2}}}{\partial t}}=-D({\frac {\phi _{{2}}-\phi _{{3}}}{h^{{2}}}})$ in ${\frac {\partial \phi _{{3}}}{\partial t}}$.

This eliminates two fluxes ${\textbf {v}}$ from our diffusion matrix. Thus, we can construct a $3\times 2$ matrix:

${\textbf {D}}_{{Dirichlet-Dirchlet}}={\begin{pmatrix}1&0\\-1&1\\0&-1\end{pmatrix}}$
where
$\phi =(\phi _{{1}},\phi _{{2}},\phi _{{3}})$
and
${\textbf {v}}=(D{\frac {\phi _{{2}}-\phi _{{1}}}{h^{{2}}}},D{\frac {\phi _{{3}}-\phi _{{2}}}{h^{{2}}}})^{{{\textbf {T}}}}$

Let $m=3$, corresponding to the number of rows of our matrix, and $n=2$, corresponding to the number of columns of our matrix.

We now perform SVD analysis on this matrix using MATLAB's SVD command, generating three matrices $U,\Sigma ,V^{{T}}$, with $U$, where $U$ is an $m\times m$ unitary matrix, $\Sigma$ is an $m\times n$ diagonal matrix containing the singular values, and $V^{{T}}$ is an $n\times n$ unitary matrix. Furthermore, using MATLAB's rank operation for $(S:D)$ gives a rank $r=2$.

Now, first, note that the dimension of the left null space is $m-r=3-2=1$, which means that the left null space has one spanning vector. Through our SVD analysis, we see from the matrix $U$ that an orthonormal basis for the left null space is composed of the vectors $U_{{r+1}},...,U_{{m}}$, which in this case is just $U_{{3}}$.

$U_{{3}}=({\frac {1}{{\sqrt {3}}}},{\frac {1}{{\sqrt {3}}}},{\frac {1}{{\sqrt {3}}}})$

This corresponds to a conserved quantity, where $\phi _{{1}}+\phi _{{2}}+\phi _{{3}}$ is time invariant, which is indeed the case, since in this approximation we have no net flux ath the left and right boundaries.

Meanwhile, the dimension of the right null space is $n-r=2-2=0$, meaning that the right null space has no spanning vectors. Biologically, there are no steady-state flux distributions, besides the trivial state in which the initial reactant concentrations are equal and no diffusive gradients form.

Dirichlet-Neumann Boundary Conditions

If we define a Dirchlet boundary condition at the left boundary of our spatial domain and a Neumann boundary condition at the right boundary, we have:

$\phi _{{0}}-\phi _{{1}}=0$ and
$D({\frac {\phi _{{4}}-\phi _{{3}}}{h^{{2}}}})=\beta$.

The presence of the constant flux at the right boundary gives us an additional flux in our system. Otherwise, this matrix is constructed analogously to the Dirichlet-Dirichlet case

Thus, we can construct a $3\times 3$ matrix:

${\textbf {D}}_{{Dirichlet-Neumann}}={\begin{pmatrix}1&0&0\\-1&1&0\\0&-1&1\end{pmatrix}}$
where
$\phi =(\phi _{{1}},\phi _{{2}},\phi _{{3}})$
and
${\textbf {v}}=(D{\frac {\phi _{{2}}-\phi _{{1}}}{h^{{2}}}},D{\frac {\phi _{{3}}-\phi _{{2}}}{h^{{2}}}},\beta )^{{{\textbf {T}}}}$

Let $m=3$, corresponding to the number of rows of our matrix, and $n=3$, corresponding to the number of columns of our matrix. MATLAB's rank operation for $(S:D)$ gives a rank $r=3$. It turns out that we do not need the SVD decomposition:

The dimension of the left null space is $m-r=3-3=0$, which means that the left null space has no spanning vectors. Biologically, there are no conserved quantities. This makes sense, due to the existence of a flux on our right hand boundary. Note that in the case of the Neumann boundary condition that specifies zero flux, the system will look exactly the same as a Dirichlet boundary condition under this approximation, and so will have a conserved quantity of the sum of all three fluxes.

Meanwhile, the dimension of the right null space is $n-r=3-3=0$, meaning that the right null space has no spanning vectors. Biologically, there are no steady-state flux distributions, besides the trivial state in which the initial reactant concentrations are equal and no diffusive gradients form.

Neumann-Neumann Boundary Conditions

If we define Neumann boundary conditions at both the left and right boundaries of our spatial domain, we have::

$D({\frac {\phi _{{0}}-\phi _{{1}}}{h^{{2}}}})=\beta _{{1}}$ and
$D({\frac {\phi _{{4}}-\phi _{{3}}}{h^{{2}}}})=\beta _{{2}}$,

with the subscripts 1 and 2 to distinguish the two potentially different fluxes at the left and right boundaries.

We have added yet another flux to the system. Thus, we can construct a $3\times 4$ matrix:

${\textbf {D}}_{{Neumann-Neumann}}={\begin{pmatrix}1&1&0&0\\0&-1&1&0\\0&0&-1&1\end{pmatrix}}$
where
$\phi =(\phi _{{1}},\phi _{{2}},\phi _{{3}})$
and
${\textbf {v}}=(\beta _{{1}},D{\frac {\phi _{{2}}-\phi _{{1}}}{h^{{2}}}},D{\frac {\phi _{{3}}-\phi _{{2}}}{h^{{2}}}},\beta _{{2}})^{{{\textbf {T}}}}$

Let $m=3$, corresponding to the number of rows of our matrix, and $n=4$, corresponding to the number of columns of our matrix. MATLAB's rank operation for $(S:D)$ gives a rank $r=3$.

The dimension of the left null space is $m-r=3-3=0$, which means that the left null space has no spanning vectors. Biologically, there are no conserved quantities. This makes sense, due to the existence of a flux on both the left and right hand boundaries. Again, if the Dirichlet boundary condition specifies zero flux at both of these ends, the Dirichlet-Dirichlet conserved quantity arises.

Meanwhile, the dimension of the right null space is $n-r=4-3=1$, meaning that the right null space has one spanning vector. Indeed, through our SVD analysis, we see from taking the transpose of the $V^{{T}}$ matrix (not reproduced here) that an orthonormal basis is composed of the vectors $V_{{r+1}},V_{{r+2}},...,V_{{n}}$, or in this case, $V_{{4}}$.

We reproduce this vector here using MATLAB's svd function:

$V_{{4}}=(-{\frac {1}{2}},{\frac {1}{2}},{\frac {1}{2}},{\frac {1}{2}})$

Remember that the Neumann boundary conditions were specified as +1 in the diffusive matrix, which corresponds to fluxes into the spatial domain. Biologically, our 'V' vector corresponds to the steady state distribution:

$-\beta _{{1}}=D{\frac {\phi _{{2}}-\phi _{{1}}}{h^{{2}}}}=D{\frac {\phi _{{3}}-\phi _{{2}}}{h^{{2}}}}=\beta _{{2}}$

This makes biological sense, as it corresponds to the situation in which each point on the spatial domain is balanced by equal and opposite diffusion fluxes. Specifically, this steady-state case really corresponds to the reactant in question diffusing across the entirety of the spatial domain at a uniform rate.