Accuracy Testing of Contact Reactions
cellular_raza 0.1.12 introduced the new simulation aspect of
ReactionsContact which can couple the
intracellular reactions of individual cells.
Testing the individual Adams-Bashforth solver which is used in the contact reactions update function
reactions_contact_adams_bashforth_2nd
is part of a different testing approach.
In contrast, these tests take a high-level view and aim to match known results for a complete system
onto results generated by cellular_raza.
Two-Component Uncoupled System
To get a general understand of what contact reactions are aiming to solve, we refer to the Reactions section. In this test, we test that a fully uncoupled system is correctly solved and no unspecified couplings are introduced during the process. Furthermore, this also tests for the accuracy of the used solver.
Theoretical Formulation
The first test ensures that a system of two components given by
$$\begin{align} \dot{x}(t) &= f(x) = \alpha\\ \dot{y}(t) &= g(y) = \alpha y \left(1 - \frac{y}{y_\text{max}} \right) \end{align}$$
is correctly solved and no unintended coupling is introduced. These equations correspond to linear growth for $x(t)$ and the well-known logistic curve [1] in $y(t)$. The analytical solutions to these equations are given by
$$\begin{align} x(t) &= \alpha(t-t_0)\\ y(t) &= y_\text{max}\left(1+ \frac{y_\text{max}-y_0}{y_0} e^{-\alpha(t-t_0)} \right)^{-1}. \end{align}$$
Implementation in cellular_raza
To model them with cellular_raza we define a new agent type ContactCell.
|
|
Notice that the mechanics: NewtonDamped1DF32 component only serves as the positional information
such that we can use the predefined
CartesianCuboid as the
simulation domain.
Furthermore, we implement the required concepts
Intracellular,
ReactionsContact as follows.
|
|
The ReactionsContact trait requires
that agents interact with each other, meaning we have to at least insert 2 agents into the
simulation to obtain any effect.
Since we employ no restraints on the range of interaction by not using the values _own_pos or
_ext_pos in the calculate_contact_increment function, every cell will interact with each other.
This also results in an increased reaction speed when more than 2 cells are present, meaning our
variable $\alpha$ for the exact solution is actually dependent on the number of species.
$$\begin{equation} \alpha = (N_\text{agents}-1)\alpha_0 \end{equation}$$
We have named the variable of the individual agents accordingly.
Solving the System
Now we are ready to solve our system with the
run_simulation! macro.
To test various configurations, we write a function which takes in all needed parameters to solve
the system.
|
|
Comparing Results
In order to meaningfully compare numerical results, we need an estimate for the local and global truncation error [2] which is introduced by our numerical solver. For a given ODE in the form of
$$\begin{equation} \dot{y} = f(t, y) \end{equation}$$
which can be solved by an algorithm $A$ in the form of
$$\begin{equation} y_n = y_{n-1} + \Delta t A(t_{n-1}, y_{n-1}, \Delta t, f) \end{equation}$$
the local truncation error $\tau_n$ is given by
$$\begin{equation} \tau_n = y(t_n) - y(t_{n-1}) - \Delta t A(t_{n-1}, y_{n-1}, \Delta t, f) \end{equation}$$
and is related to the global truncation error $e_n$ via
$$\begin{align} e_n &= y(t_n) - y_n\\ &= y(t_n) - (y_0 + \Delta A(t_0, y_0, \Delta t, f) + \dots + \Delta t A(t_{n-1},y_{n-1},\Delta t, f). \end{align}$$
Under the additional condition that $f$ is Lipschitz with Lipschitz-constant $L$, we can derive a bound for the global error.
$$\begin{equation} |e_n| \leq \frac{\text{max}_j \tau_j}{\Delta t L}\left(e^{L(t-t_0)} + 1 \right) \end{equation}$$
where $\tau_n$ is the local truncation error at time step $n$.
By finding an upper bound for it, we can further simplify this formula.
The local truncation error $\tau_n$ depends on the type of solver used.
In our case, we employ the Adams-Bashforth [3] solver with 3rd order.
The executed function in the chili backend is
contact_reactions_adams_bashforth_3rd.
The local truncation error is bound
$$\begin{equation} |\tau_n| \leq \frac{3\Delta t^4}{8}\sup\limits_{t\in(t_n-2\Delta t,x_n+\Delta t)}|y^{(4)}(t)| \end{equation}$$
where $y^{(4)}(t)$ is the fourth-order derivative of the analytical solution $y$. We can calculate this value for both ODEs by differentiating equations $(3)$ and $(4)$.
$$\begin{align} x^{(n)}(t) &= 0\\ y^{(n)}(t) &= n! y_\text{max} \alpha^n\left(\frac{y_\text{max}-y_0}{y_0}\right)^n \left(1+\frac{y_\text{max}-y_0}{y_0}e^{-\alpha(t-t_0)}\right)^{-(n+1)}\\ \end{align}$$
With this upper bound on the local truncation error we can construct a new function which tests that this upper bound is fulfilled. We begin by writing down a general formula for the nth drivative of the exact solution to the logistic curve ODE problem as given above. Notice that we deliberately use arguments for initial values and formulate the solution for the most general case.
|
|
In the next step we calculate the Lipschitz-constants $L_0,L_1$ and use them together with the fourth derivative bound to calculate the local and global truncation errors. The Lipschitz condition is given by
$$\begin{equation} |f(s) - f(t)| \leq L_0 |s-t| \end{equation}$$
In our case, the Lipschitz constant for $f(x)$ could be taken to be indentiaclly $0$. This would introduce a constant global error in the limit $L_0\rightarrow 0$. However, due to numerical uncertainties in even addition operations, we opt to take the more conservative approach and choose $L_0 = \alpha$.
We calculate Lipschitz-constant $L_1$ for $g(y)$ given by equation $(2)$.
$$\begin{align} |g(y)-g(z)| &= \left|\alpha y\left(1-\frac{y}{y_\text{max}}\right) - \alpha z\left(1-\frac{z}{y_\text{max}}\right) \right|\\ &=\alpha |y-z|\left|\left(1 - \frac{y+z}{y_\text{max}}\right) \right|\\ &\leq|y-z|\text{max}\left(1-\frac{2y_0}{y_\text{max}},1-\frac{2y(t_\text{max})}{y_\text{max}}\right)\\ L_1 &= \text{max}\left(1-\frac{2y_0}{y_\text{max}},1-\frac{2y(t_\text{max})}{y_\text{max}}\right) \end{align}$$
|
|
In the final step, we use the already defined function `` to generate results from cellular_raza
and compare them with the exact known results.
Their difference has to be within the margin of the calculated global truncation error $e$.
Note that we only start comparing after the first 3 initial steps which are performed by solvers
of lower orders (Euler and Adams-Bashforth 2nd) due to insufficient information about increments
of previous integration steps.
Afterwards, we use the numerically calculated value at $t=t_0 + 2\Delta t$ as the new initial value.
|
|
In the last step, we also added a function save_results which exports all generated results to a
.csv file which.
|
|
Equipped with this function compare_results, we can now use it for a collection of configurations
to test the solver.
|
|
The generated files can then be used by a small python script to create the plots seen below.
|
|
The full code can be found under
cellular_raza/tests/contact_reactions.rs.
To run it yourself, clone the git repo git clone https://github.com/jonaspleyer/cellular_raza and
execute cargo test -- two_component_contact_reaction in the cellular_raza repository directory.
Results
The following plots she the calculated results from cellular_raza and the analytical solution of
the second ODE described above in equation $(2)$ for various configurations of parameters.
They have been chosen in such a way to capture different dynamics and time-scales.
The global truncation error is used as errorbars for the analytical solution.
In the first image, we can clearly see the exponential growth of the global truncation error $e$
over time.
The size of the errorbars is mainly determined by the time interval $\Delta t$ chosen to solve the
equations.
The second studied example config1 shows how a low step-size $\Delta t$ yields results which
follow the trajectory of the analytical solution to an even higher margin of error than before.
Discussion
We have shown how to derive useful bounds for the local and global truncation error of a result
with a known analytical solution.
Furthermore, these mathematical results have been implemented in an automated testing scheme to
measure and verify the solvers behind the ReactionsContact simulation aspect of cellular_raza.
Results generated initially are of lower accuracy. This can be explained by the internals of our implementation of the Adams-Bashforth [3] solver. The numerical solver does not have any knowledge about the increments $\Delta y_n$ before the simulation has started (ie. $n<0$). It is thus forced to assume nothing and reverts back to the 1st order case. For the second simulation step, the 2nd order Adams-Bashforth solver can be used and only afterwards do we have enough information about previous time increments such that the full 3rd order solver can take over. This behaviour only occurrs in the very initial steps of our numerical routine.
References
[1] P. F. Verhulst, “Recherches mathématiques sur la loi d’accroissement de la population.” Nouveaux mémoires de l’Académie Royale des Sciences et Belles-Lettres de Bruxelles, vol. 18, pp. 14–54, 1845, Available: http://eudml.org/doc/182533
[2] E. Süli and D. F. Mayers, “An Introduction to Numerical Analysis.” Cambridge University Press, Aug. 28, 2003. doi: 10.1017/cbo9780511801181
[3] Bashforth, Francis and Adams, J. Couch An attempt to test the theories of capillary action : by comparing the theoretical and measured forms of drops of fluid. Cambridge [Eng.]: Cambridge University Press, 1883
[4] G. Fasshauer, “Numerical Methods for Differential Equations/Computational Mathematics II,” 2007, Available: http://www.math.iit.edu/~fass/478578_Chapter_2.pdf