# Numerical implementation of the deterministic model resolution
We present here the numerical implementation of the deterministic model presented in L. Darrigade and al. All notations are presented in the article. We recal porous media type equations use to model cell densities and reaction-diffusion equations for metabolites concentrations.
## Porous media type equations
We aim at solving the porous media type equations :
\begin{equation}
\begin{cases}
\partial_t \rho_l - W \partial_z \left(\phi(z) \rho_l \partial_z \rho \right) = \sum_k \epsilon_{k,l} q_k^\infty q_k(z, l, D \rho, c_b), \\
\partial_t \rho_l - W \partial_z \left(\phi(z) \rho_l \partial_z \rho \right) = \sum_k \eta_{k,l} q_k^\infty q_k(z, l, D \rho, c_b), \\
\rho_l(0,t) = 0 \text{ if $l\neq sc$,}\\
\rho_{sc}(0,t) = \rho_{sc}^{bot}, \\
\partial_z \rho_l(z_{max}, t) = 0.
\end{cases}
\label{eq:PME}
\end{equation}
for $ l \in \{'sc','pc','gc','ent'\}$. However this type of equations are not trivial to solve due to the non-linear cross diffusion terms. To be able to simulate we use the kinetic diffusive explicit scheme [(D.Aregba-Driollet & all 2003)](https://www.ams.org/journals/mcom/2004-73-245/S0025-5718-03-01549-7/) and transform the previous equation in the following form :
with $\hat{\rho_l} = \sum_{k \ne l} \rho_k$ and $G(z) = \frac{1}{2} W \phi(z) $. This new transport-diffusion form is adapted to the use of the kinetic diffusive scheme to obtain numerical simulations.
### Definition of $\phi(z)$, $G(z)$ and space regulation term of cell fate events
We recal that $\phi(z)$ allows to take into account the impact of curvature on the top and at the bottom of the crypte on the interaction force, its value is given by :
\begin{equation}
\label{eq:phi}
\mathcal{\phi}(z):=\begin{cases}
\frac{\frac{f(z)}{r_0}-\frac{f(0)}{r_0}}{1-\frac{f(0)}{r_0}}\quad &\text{ if }z\leq r_0-\epsilon\\
1 \quad &\text{ if }r_0-\epsilon<z<z_{max}-r_0+\epsilon\\
\frac{\frac{f(z_{max})}{r_0}-\frac{f(z)}{r_0}}{\frac{f(z_{max})}{r_0}-1}\quad &\text{ if }z\geq z_{max}-r_0+\epsilon
\end{cases}
\end{equation}
with $f(z)$ the crypt form :
\begin{equation}\label{eq:def_f}
f(z)=\begin{cases}
r_0\sqrt{\frac{z+\epsilon}{r_0}\left(2-\frac{z+\epsilon}{r_0}\right)}\quad &\text{ if }z\leq r_0-\epsilon\\
r_0 \quad &\text{ if }r_0-\epsilon<z<z_{max}-r_0+\epsilon\\
r_0\left(2-\sqrt{\frac{z_{max}-z+\epsilon}{r_0}\left(2-\frac{z_{max}-z+\epsilon}{r_0}\right)}\right)\quad &\text{ if }z\geq z_{max}-r_0+\epsilon
\end{cases}
\end{equation}
In the implementation of the resolution we supposed the existence of a correction term to take into acount the impact of curvature on the local density of cells : $C(z)$. This would give $G(z) = \frac{W \phi(z) C(z)}{2}$.
In the article (\bold{REF} )this approach is not explored so we assume here that C(z) = 1.
concerning space regulation term of cell fat events we use the piecewise polynomial function $R(y,K,l)$ such as :
\begin{equation}
R(y,K,\ell) :=
\label{form:regul}
\begin{cases}
0 & \text{if } y \leq K - \ell,\\
-\frac{1}{4\ell^3}y^3 + \frac{3K}{4\ell^3}y^2 - \frac{3K^2 - 3\ell^2}{4l^3} y + \frac{K^3 + 2\ell^3 - 3 K\ell^2}{4\ell^3} & \text{if } K - \ell \leq y \leq K + \ell,\\
1 & \text{if } K + \ell \geq y.
\end{cases}
\end{equation}
Several regulations are not time-dependent but only space-dependent, as in the case of Wnt regulations (see table 1)(\bold{REF})
Space variables are stored in a dictionary : dict_space.
Time resolution until steady-state. At the final time get the resolution variable $y$, the number of iterations necessary for resolution, densities of cell living a cell fate event (div, diff, ext) and the final time $t$.
rho_cell_fate is configured as : \[Div_sc, Diff_sc_pc, Div_pc, Diff_pc_gc, Diff_pc_ent, Ext_gc_ent\]
