Daphne implementation

ECM Laplacian

(ManifoldRing.NodeInterpolator)ScalarField Laplacian(ScalarField sf)
(ManifoldRing.Trilinear3D)LocalMatrix[][] laplacianMatrix()

Array representations

In the current implementation of Daphne, all scalar fields are represented by arrays which will be designated by upper case letters. In the case of a moment-expansion scalar field \(a({\bf y}) = a^{(0)} + {\bf a^{(1)}} \cdot {\bf y}\), the array representation has 4-elements \(A = \left[ A_0, A_1, A_2, A_3 \right]\), where \(A_0 = a^{(0)}\) and \(A_i\) is the \(i\text{th}\) component of \({\bf a^{(1)}}\).

Interpolated scalar fields are represented by \(N\)-element arrays \(A = \left[ A_0, A_1, \ldots, A_N \right]\), where \(A_i\) is the value of the scalar field at the \(i\text{th}\) lattice point.

It is sometimes more convenient to refer to the lattice points in a matrix format \(A(i,j,k) \equiv A_n\), where \(i \in \{0,1,\ldots, n_1-1\}, j \in \{0,1,\ldots, n_2-1\} , k\in \{0,1,\ldots, n_3-1\}\) and \(n= i + j n_1 + k n_1 n_2\)

ECS:

Using the matrix form of the interpolated scalar field,

(1)\[\begin{split} \mathcal{L}(S(i,j,k)) & = & \frac{1}{\Delta^2} \sum_{p=1}^{3} \left[ S(i\!+\!\delta_{p1},j\!+\!\delta_{p2}, k\!+\!\delta_{p3}) + S(i\!-\!\delta_{p1},j\!-\!\delta_{p2}, k\!-\!\delta_{p3}) - 2 S(i,j,k) \right] + \frac{2}{\Delta} F(i,j,k),\end{split}\]

where

\(\Delta\) is the grid size, \(\delta_{ij}\) is the Kroenecker delta, and \(F\) is the flux.

Cytoplasm (TinyBall):

(2)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} \mathcal{L}(s) & = & -\frac{5}{\rho^2} {\bf s}^{(1)} = -\frac{5}{\rho^2} [0, S_1, S_2, S_3 ], \label{eq.laplacian TinyBall} \end{array}\end{split}\]

with moment-expansion flux \(\phi\) and array flux \(F\equiv D \phi\) at the boundaries according to

(3)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} \phi & = & \frac{3}{D \rho} c^{(0)} + \frac{5}{D \rho} {\bf c^{(1)}} \\ & \equiv & F = \frac{1}{D \rho} [ 3 F_0 , 5 F_1, 5 F_2, 5 F_3 ] \label{eq.laplacian TinyBall flux} \end{array}\end{split}\]

and array flux (\(F=D\phi\))

Plasma membrane (TinySphere):

(4)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} \mathcal{L}(s) & = & \frac{2}{\rho^2} {\bf c^{(1)}} = \frac{2}{\rho^2} [0, S_1, S_2, S_3 ] \label{eq.laplacian TinySphere} \end{array}\end{split}\]
(5)\[\begin{split} \renewcommand{\arraystretch}{2.0} \begin{array}{rcl} \nabla^2 c({\bf x}) & \approx & l^{(0)} + {\bf l^{(1)} } \cdot {\bf x} \\ l^{(0)} & = & 0\\ {\bf l^{(1)} } & = & -\frac{2}{\rho^2} {\bf c^{(1)}}, \end{array}\end{split}\]

Reactions

Reactions are described in terms of reactant, product, and modifier molecules, where modifier molecules participate in, but are not changed by, the reaction (eg, enzymes).

Bulk reactions

Reactions whose reactants, products, and modifiers are in the same space/manifold \(V\) (or \(B\) or \(S\)) are designated as Bulk Reactions.

Let us use catalyzed association as an example of a Bulk Reaction. Reactant molecules \(A\) and \(B\) in \(V\) bind to form product molecule \(C\), also in \(V\)

(6)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} A + B & \xrightarrow{M} & C . \end{array}\end{split}\]

The reaction requires the presence of modifier molecule \(M\), which is not altered by the reaction.

Boundary reactions

Reactions between molecules in \(V\) (or \(B\)) and molecules in \(S\) are designated as Boundary Reactions.

Let us use receptor/ligand binding as an example of a Boundary Reaction. We have a cell centered at position \({\bf Y_c}\) in the ECS. Ligand \(L\) in \(V\) binds to receptor \(R\) in \(S\) to form complex \(C\) in \(S\).

(7)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} R + L & \xrightarrow{} & C . \end{array}\end{split}\]

In the local (cell-centered) coordinate system, the rate of change of complex \(C\) is given by

(8)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} \frac{d c({\bf x})}{dt} & = & k r({\bf x}) l({\bf x}) , \end{array}\end{split}\]

where the moment-expansion representation of the receptor and complex concentrations on the cell surface are given by

(9)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} r({\bf x}) & \approx & r^{(0)} + {\bf r^{(1)}} \cdot {\bf x} \\ c({\bf x}) & \approx & c^{(0)} + {\bf c^{(1)}} \cdot {\bf x}. \end{array}\end{split}\]

and, using the results of (?), the concentration of the ligand in the ECS at the surface of the cell \(L({\bf Y(x)})\) is given by

(10)\[ \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} l({\bf x}) \approx l^{(0)} + {\bf l^{(1)}} \cdot {\bf x}, \end{array}\]

where \(l^{(0)}=L({\bf Y_c})\) and \({\bf l^{(1)}}=\nabla L({\bf Y_c})\).

In the local (cell-centered) coordinate system, the rate of change of ligand is given by

(11)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} \frac{d l({\bf x})}{dt} & = & -k r({\bf x}) l({\bf x}). \end{array}\end{split}\]

The evolution equation for the zeroth moment is

(12)\[\begin{split} \renewcommand{\arraystretch}{2.0} \begin{array}{rcl} \frac{d l^{(0)} }{dt} & = & - \frac{1}{4\pi \rho^2} \int_{S} dS\; k r({\bf x}) l({\bf y(x)}) \\ & \approx & - \frac{1}{4\pi \rho^2} \int_{S} dS\; k \left( r^{(0)} + {\bf r^{(1)}} \cdot {\bf x} \right) \left( l^{(0)} + {\bf l^{(1)}} \cdot {\bf x} \right) \\ & \approx & - \frac{1}{4\pi \rho^2} \int_{S} dS\; k \left[ r^{(0)} l^{(0)} + l^{(0)} {\bf r^{(1)}} \cdot {\bf x} + r^{(0)} {\bf l^{(1)}} \cdot {\bf x} + \left( {\bf r^{(1)}} \cdot {\bf x} \right) \left( {\bf l^{(1)}}\cdot {\bf x} \right) \right] \\ & = & - k r^{(0)} l^{(0)} - \frac{1}{4\pi \rho^2} \int_{S} dS\; k \left( {\bf r^{(1)}} \cdot {\bf x} \right) \left( {\bf l^{(1)}} \cdot {\bf x} \right) \\ & = & - k r^{(0)} l^{(0)} , \end{array}\end{split}\]

where terms of quadratic order have been neglected.

The evolution equation for the first moment is

(13)\[\begin{split} \renewcommand{\arraystretch}{2.0} \begin{array}{rcl} \frac{d {\bf l^{(1)}} }{dt} & = & - \frac{3}{4\pi \rho^4} \int_{S} dS\; k r({\bf x}) l({\bf y(x)}) \,{\bf x} \\ & \approx & - \frac{3}{4\pi \rho^4} \int_{S} dS\; k \left( r^{(0)} + {\bf r^{(1)}} \cdot {\bf x} \right) \left( l^{(0)} + {\bf l^{(1)}} \cdot {\bf x} \right) \,{\bf x} \\ & \approx & - \frac{3}{4\pi \rho^4} \int_{S} dS\; k \left[ r^{(0)} l^{(0)} + l^{(0)} {\bf r^{(1)}} \cdot {\bf x} + r^{(0)} {\bf l^{(1)}} \cdot {\bf x} + \left( {\bf r^{(1)}} \cdot {\bf x} \right) \left( {\bf l^{(1)}}\cdot {\bf x} \right) \right] {\bf x} \\ & \approx & - \frac{3}{4\pi \rho^4} \int_{S} dS\; k \left[ l^{(0)} {\bf r^{(1)}} \cdot {\bf x} + r^{(0)} {\bf l^{(1)}} \cdot {\bf x} \right] {\bf x} \\ & \approx & - \frac{3}{4\pi \rho^2} \int_{S} dS\; k \left[ l^{(0)} {\bf r^{(1)}} \cdot {\bf \hat{x}} + r^{(0)} {\bf l^{(1)}} \cdot {\bf \hat{x}} \right] {\bf \hat{x}} \\ & \approx & - \frac{3}{4\pi \rho^2} \frac{4\pi}{3} \rho^2 k \left[ l^{(0)} {\bf r^{(1)}} + r^{(0)} {\bf l^{(1)}} \right] \\ & \approx & - k \left[ l^{(0)} {\bf r^{(1)}} + r^{(0)} {\bf l^{(1)}} \right] \end{array}\end{split}\]

where terms of quadratic order have been neglected.

Same as moment-expansion multiplication.

At one time-step \(\delta t\) the changes to the receptor and complex concentrations due to the boundary association reaction are given by

(14)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} \delta r(t) & = & - \delta t \, k r(t) l(t) = - \delta t \, k \left[ R_0 L_0, R_0 L_1+R_1L_0, R_0 L_2+R_2L_0, R_0 L_3+R_3L_0,\right] \\ \delta c(t) & = & \delta t \, k r(t) l(t) = \delta t \, k \left[ R_0 L_0, R_0 L_1+R_1L_0, R_0 L_2+R_2L_0, R_0 L_3+R_3L_0,\right] \end{array}\end{split}\]

where \(r\), \(c\), and \(l\) are represented by 4-element arrays and the rules for moment-expansion addition (?) and multiplication (?) hold. Note that a non-uniform ligand concentration around the position of the cell (non-zero gradient) will result in non-uniform receptor and complex concentrations even if they were initially uniform, due to the rules governing moment-expansion multiplication and addition (see section ref{sec.moment expansion}).

The above reaction (7) also produces a flux of ligand at the ECS/sphere surface

(15)\[\begin{split} \renewcommand{\arraystretch}{1.0} \begin{array}{rcl} \phi(t) & = & \delta t \, D^{-1} k \, r(t) l(t) \\ %\delta t \, D^{-1} k \, r(t) l(t) \\ f & = & \delta t \, k \left[ R_0 L_0, R_0 L_1+R_1L_0, R_0 L_2+R_2L_0, R_0 L_3+R_3L_0,\right] \label{eq.reaction flux} \end{array}\end{split}\]

where \(D\) is the diffusion coefficient of the ligand in the ECS, and \(f\) is the representation of the flux in a 4-element moment-expansion array. The units of \(f\) are

(16)\[\begin{split} \renewcommand{\arraystretch}{2.0} \begin{array}{rcl} [f_0] & = & \left[ \delta t \, k \, R_0 L_0 \right] = T \frac{L^3}{T \#} \frac{\#^2}{L^5} = \frac{\#}{L^2} \\ \left[ f_i \right] & = & \frac{\#}{L^3} \;\; (i=1,2,3) . \end{array}\end{split}\]

Note that the flux is positive for this reaction, which is consistent with the fact that it produces a net loss of ligand from the ECS. (see section ref{subsec.flux})

At each time step \(\delta t\):

The ligand concentration \(L\) in the ECS is updated accoring to

(17)\[\begin{split} \renewcommand{\arraystretch}{1}\begin{array}{rcl} L(i,j,k) & += & \frac{2\delta t}{\Delta} F[0] - \frac{D\delta t}{\Delta^2} \sum_{p=1}^{3} \left[ L(i\!+\!\delta_{p1},j\!+\!\delta_{p2}, k\!+\!\delta_{p3}) + L(i\!-\!\delta_{p1},j\!-\!\delta_{p2}, k\!-\!\delta_{p3}) - 2 L(i,j,k) \right] . \end{array}\end{split}\]

Checking units:

(18)\[\begin{split} \renewcommand{\arraystretch}{2}\begin{array}{rcl} \left[ L(i,j,k) \right] & = & \frac{\#}{L^3} \\ \left[ \frac{2\delta t}{\Delta} F[0] \right] & =& \frac{T}{L} \frac{\#}{T L^2} = \frac{\#}{L^3} \\ \left[ \frac{D\delta t}{\Delta^2} L(i,j,k) \right] & = & \frac{L^2}{T} \frac{T}{L^2} \frac{\#}{ L^3} = \frac{\#}{L^3} \end{array}\end{split}\]

Plasma membrane (TinySphere):

(19)\[\begin{split} \renewcommand{\arraystretch}{2.0} \begin{array}{rcl} R & += & - \frac{2D\delta t }{\rho^2} [0, R_1, R_2, R_3 ] - \delta t \, k \left[ R_0 L_0, R_0 L_1+R_1L_0, R_0 L_2+R_2L_0, R_0 L_3+R_3L_0,\right] \\ C & += & - \frac{2D\delta t }{\rho^2} [0, C_1, C_2, C_3 ] + \delta t \, k \left[ R_0 L_0, R_0 L_1+R_1L_0, R_0 L_2+R_2L_0, R_0 L_3+R_3L_0,\right] \end{array}\end{split}\]

Checking units:

(20)\[\begin{split} \renewcommand{\arraystretch}{2}\begin{array}{rcl} \left[ R_0 \right] & = & \frac{\#}{L^2} \;\;\mbox{ and } \left[ R_i \right] = \frac{\#}{L^3} \;\; (i=1,2,3) \\ \left[ \frac{2D\delta t }{\rho^2} R_i \right] & =& \frac{L^2 T}{T L^2} \frac{\#}{ L^3} = \frac{\#}{L^3} \\ \left[ \delta t \, k R_0 L_0 \right] & = & T \frac{L^3}{T\#} \frac{\#^2}{ L^5} = \frac{\#}{L^2} \end{array}\end{split}\]

MolecularPopulation.Step()

concentration += dt * Molecule.DiffusionCoefficient * concentration.Laplacian();

\([\mbox{dt*Molecule.DiffusionCoefficient}] = L^2\)

  1. ECS:

    1. \([\mbox{concentration}] = \frac{\#}{L^3}\)

    2. \([\mbox{concentration.Laplacian()}] = \frac{\#}{L^5}\)

    3. :math:`[mbox{dt * Molecule.DiffusionCoefficient * concentration.Laplacian()}]

      = frac{#}{L^3}`

  2. TinyBall

    1. \([\mbox{concentration}] = \frac{\#}{L^3}, \frac{\#}{L^4}\)

    2. \([\mbox{concentration.Laplacian()}] = \frac{\#}{L^5}, \frac{\#}{L^6}\)

    3. :math:`[mbox{dt * Molecule.DiffusionCoefficient * concentration.Laplacian()}]

      = frac{#}{L^3}, frac{#}{L^4}`

  3. TinySphere

    1. \([\mbox{concentration}] = \frac{\#}{L^2}, \frac{\#}{L^3}\)
    2. \([\mbox{concentration.Laplacian()}] = \frac{\#}{L^4}, \frac{\#}{L^5}\)
    3. \([\mbox{dt * Molecule.DiffusionCoefficient * concentration.Laplacian()}] = \frac{\#}{L^2}, \frac{\#}{L^3}\)
concentration += -dt * concentration.DiffusionFluxTerm(kvp.Value, compartment.BoundaryTransforms[kvp.Key]);

  1. ECS:

    1. \([\mbox{concentration}] = \frac{\#}{L^3}\)
    2. \([\mbox{concentration.DiffusionFluxTerm()}] = \frac{\#}{TL^3}\)
    3. \([\mbox{dt * concentration.DiffusionFluxTerm.Laplacian()}] = \frac{\#}{L^3}\)

  2. TinyBall

    1. \([\mbox{concentration}] = \frac{\#}{L^3}, \frac{\#}{L^4},\)
    2. \([\mbox{concentration.DiffusionFluxTerm()}] = \frac{\#}{TL^3}, \frac{\#}{TL^4}\)
    3. \([\mbox{dt * concentration.DiffusionFluxTerm()}] = \frac{\#}{L^3}, \frac{\#}{L^4}\)

  3. TinySphere - not implemented