Degenerate area preserving surface Allen-Cahn equationMichal Beneš, Miroslav
Kolář, Jan M. Sischka, and Axel Voigt

# Degenerate area preserving surface Allen-Cahn equation and its sharp
interface limit

Michal Beneš Department of Mathematics, Faculty of Nuclear Sciences and
Physical Engineering, Czech Technical University in Prague, 12000 Prague,
Czech Republic<EMAIL_ADDRESS><EMAIL_ADDRESS>Miroslav Kolář11footnotemark: 1 Jan Magnus Sischka Institut für
Wissenschaftliches Rechnen, Technische Universität Dresden, 01062 Dresden,
Germany (, ) JMS and AV were supported by the German Research Foundation (DFG)
through EXC PoL. MB and MK were partly supported by the project 21-09093S of
the Czech Science Foundation. We further acknowledge computing resources
provided by ZIH at TU Dresden within projects WIR. jan_magnus.sischka@tu-
dresden.de<EMAIL_ADDRESS>Axel Voigt22footnotemark: 2

###### Abstract

We consider formal matched asymptotics to show the convergence of a degenerate
area preserving surface Allen-Cahn equation to its sharp interface limit of
area preserving geodesic curvature flow. The degeneracy results from a surface
de Gennes-Cahn-Hilliard energy and turns out to be essential to numerically
resolve the dependency of the solution on geometric properties of the surface.
We experimentally demonstrate convergence of the numerical algorithm, which
considers a graph formulation, adaptive finite elements and a semi-implicit
discretization in time, and uses numerical solutions of the sharp interface
limit, also considered in a graph formulation, as benchmark solutions.

###### keywords:

motion by geodesic curvature, surface Allen-Cahn equation, de Gennes-Cahn-
Hilliard energy, matched asymptotic expansion, graph formulation

35K57, 53E10

##  1 Introduction

The connection between phase field approximations and geometric partial
differential equations is well established and can formally be justified by
matched asymptotics, see [1]. Geometric partial differential equations are
evolution equations that evolve curves or surfaces according to their
curvature. Similarly to these curvature driven flows in 2D or 3D one can
consider the evolution of curves on surfaces. The evolution of these curves is
governed by geodesic curvature and thus in addition strongly depending on the
local geometric properties of the underlying surface. First analytical
attempts to connect these geodesic evolution laws to surface phase field
models have been considered in [2, 3, 4, 5]. We here show this connection for
a degenerate area preserving surface Allen-Cahn equation and an area
preserving geodesic curvature flow. The surface Allen-Cahn equation provides
the mathematical basis to study curvotaxis of cells in epithelia tissue [6].

We consider a surface de Gennes-Cahn-Hilliard energy

(1)
${\cal{F}}_{dGCH}(\phi)=\tilde{\sigma}\int_{\cal{S}}\frac{1}{G(\phi)}\left(\frac{\epsilon}{2}\|\nabla_{\cal{S}}\phi\|^{2}+\frac{1}{\epsilon}W(\phi)\right)d{\cal{S}}$

with surface ${\cal{S}}$, phase field variable $\phi$, surface gradient
$\nabla_{\cal{S}}$, double well potential
$W(\phi)=\frac{1}{4}(\phi^{2}-1)^{2}$, rescaled surface tension
$\tilde{\sigma}$ and small parameter $\epsilon>0$ determining the thickness of
the diffuse interface. The factor $1/G(\phi)$ is called the de Gennes
coefficient in polymer science. We consider $G(\phi)=\frac{3}{2}(\phi^{2}-1)$
or a regularized version
$G_{\eta}(\phi)=(\frac{9}{4}(\phi^{2}-1)^{2}+\eta^{2}\epsilon^{2})^{1/2}$ with
$\eta>0$. The scaling coefficient is such that the sharp interface limit
equals the one obtained from the usual Cahn-Hilliard energy without the de
Gennes coefficient [7, 8]. Evolution equations based on this energy, at least
in flat space, have been shown numerically advantageous, as the singularity,
$G(\phi)$ or $G_{\eta}(\phi)$, helps to keep solutions confined in $[-1,1]$.
Even though the theoretical foundation of this argument remains open, several
numerical studies confirm this [9, 10]. The resulting degenerate area
preserving surface Allen-Cahn equation on ${\cal{S}}$ reads

(2) $\displaystyle\epsilon\tilde{\beta}G(\phi)\partial_{t}\phi$
$\displaystyle=$
$\displaystyle\tilde{\sigma}\left(\epsilon\Delta_{\cal{S}}\phi-\frac{1}{\epsilon}W^{\prime}(\phi)\right)+G(\phi)\lambda,$
(3) $\displaystyle\lambda$ $\displaystyle=$
$\displaystyle\tilde{\sigma}\frac{1}{|{\cal{S}}|}\int_{\cal{S}}\left(-\frac{\epsilon}{G(\phi)}\Delta_{\cal{S}}\phi+\frac{1}{\epsilon
G(\phi)}W^{\prime}(\phi)\right)\,d{\cal{S}},$

with initial condition $\phi(0)=\phi_{0}$. $\Delta_{\cal{S}}$ denotes the
Laplace-Beltrami operator and $\tilde{\beta}>0$ is a rescaled kinetic
coefficient. $\lambda$ is a Lagrange multiplier for the area constraint
$\frac{1}{|{\cal{S}}|}\int_{\cal{S}}\phi\;d{\cal{S}}=\alpha$ with
$\alpha\in[-1,1]$. The model follows from eq. Eq. 1 as a constrained
$L^{2}$-gradient flow and considers the asymptotic approximation
$\frac{\epsilon}{2}|\nabla_{\cal{S}}\phi|^{2}\approx\frac{1}{\epsilon}W(\phi)$,
see [7]. The idea for this approximation was first used for phase field
approximations of surface diffusion in [11], where $G(\phi)$ was introduced as
a stabilizing function. The zero-levelset of the solution of eqs. (2) and (3)
provides an approximation of a curve $\gamma$ evolving on the surface
${\cal{S}}$.

The area preserving geodesic curvature flow for the curve $\gamma$ on
${\cal{S}}$ reads

(4)
$\beta{\cal{V}}=-\sigma{\cal{H}}_{\gamma}+\sigma\frac{1}{|\gamma|}\int_{\gamma}{\cal{H}}_{\gamma}\;d\gamma$

with initial condition $\gamma(0)=\gamma_{0}$. ${\cal{V}}$ is the velocity of
$\gamma(t)$ in the direction of the co-normal $\mu$ and ${\cal{H}}_{\gamma}$
is the geodesic curvature of $\gamma$. We consider
${\cal{S}}={\cal{S}}^{1}(t)\cup\gamma(t)\cup{\cal{S}}^{2}(t)$ and enforce the
constraint $|{\cal{S}}^{1}(t)|-|{\cal{S}}^{2}(t)|+\alpha|{\cal{S}}|=0$ for
each $t\in[0,T]$. The kinetic coefficient and line tension are related by
$\beta=\frac{4\sqrt{2}}{5}\tilde{\beta}$ and
$\sigma=\frac{2\sqrt{2}}{3}\tilde{\sigma}$.

Besides the connection between eqs. (2) and (3) and eq. (4) by formal matched
asymptotics, we use numerical solutions of eq. (4) in a graph formulation, see
[12], as benchmark problems for a numerical approach to eqs. (2) and (3),
again using a graph formulation. We use adaptive finite elements to discretize
in space and a semi-implicit time-stepping scheme.

##  2 Matched asymptotic analysis

We closely follow [5] in the analysis of a phase field model in the context of
two-phase biomembranes, and use the tools introduced in [2] to extend the
formal matched asymptotics for the area preserving Allen-Cahn equations in
flat space [13] and for the de Gennes-Cahn-Hilliard energy in flat space [7]
to surfaces. We demonstrate that eqs. (2) and (3) formally converges to eq.
(4) for $\epsilon\to 0$. We therefore require $\gamma(t)$ to be a $C^{1}$
closed curve.

###  2.1 Expansions and matching conditions

By $(\phi_{\epsilon},\lambda_{\epsilon})$, we denote a family of solutions of
eqs. (2) and (3) that converge formally to some limit denoted by
$(\phi,\lambda)$. We assume that $\phi=\chi_{\gamma}$ with
$\chi_{\gamma}:BV({\cal{S}})\to\mathbb{R}$ with $\chi_{\gamma}=-1$ on
${\cal{S}}^{1}$ and $\chi_{\gamma}=1$ on ${\cal{S}}^{2}$ for some smooth curve
$\gamma$ that separates the regions
${\cal{S}}^{1}=\\{(x,t)\in{\cal{S}}\times[0,T]:\phi(x,t)=-1\\}$ and
${\cal{S}}^{2}=\\{(x,t)\in{\cal{S}}\times[0,t]:\phi(x,t)=1\\}$. We consider an
outer and an inner expansion

(5) $\displaystyle f_{\epsilon}(x,t)$ $\displaystyle=$ $\displaystyle
f_{0}(x,t)+\epsilon f_{1}(x,t)+\epsilon^{2}f_{2}(x,t)+\ldots$ (6)
$\displaystyle F_{\epsilon}(x,t)$ $\displaystyle=$ $\displaystyle
F_{0}(s,z,t)+\epsilon F_{1}(s,z,t)+\epsilon^{2}F_{2}(s,z,t)+\ldots,$

respectively, with $f_{\epsilon}(x,t)=F_{\epsilon}(s,z,t)$, $z=r/\epsilon$ and
$\Theta(s,r,t)$ a parametrization such that $s\to\Theta(s,0,t)$ is a
paramtetrization of $\gamma(t)$ on ${\cal{S}}$ and $r$ is the signed geodesic
distance of $x=\Theta(s,r,t)\in{\cal{S}}$ to $\gamma(t)$. Thereby,
$f_{\epsilon}=\phi_{\epsilon},\lambda_{\epsilon}$ and
$F_{\epsilon}=\Phi_{\epsilon},\Lambda_{\epsilon}$ and the outer expansion
holds away from $\gamma(t)$ and the inner expansion near $\gamma(t)$. In
regions where both expansions are valid the matching conditions hold

(7) $\displaystyle F_{0}(s,\pm\infty,t)$ $\displaystyle=$ $\displaystyle
f_{0}^{\pm}(x,t),$ (8) $\displaystyle\partial_{z}F_{0}(s,\pm\infty,t)$
$\displaystyle=$ $\displaystyle 0,$ (9)
$\displaystyle\partial_{z}F_{1}(s,\pm\infty,t)$ $\displaystyle=$
$\displaystyle\nabla_{\cal{S}}f_{0}^{\pm}(x,t)\cdot\mu(x,t),$

with $f_{0}^{\pm}(x,t)=\lim_{\delta\to 0}f(\Theta(s,\pm\delta,t),t)$.

###  2.2 Outer solution

Considering the terms of $O(\epsilon^{-1})$ in eq. (2) leads to
$W^{\prime}(\phi_{0})=0$ and thus

(10) $\phi_{0}=\pm 1.$

###  2.3 Inner solution

The Laplace-Beltrami operator and the time derivative of $F_{\epsilon}(s,z,t)$
gives

(11) $\displaystyle\Delta_{\cal{S}}F_{\epsilon}$ $\displaystyle=$
$\displaystyle\frac{1}{\epsilon^{2}}\partial_{zz}F_{\epsilon}+\frac{{\cal{H}}_{\gamma}}{\epsilon}\partial_{z}F_{\epsilon}+\partial_{ss}F_{\epsilon}+O(\epsilon)$
(12) $\displaystyle\frac{d}{dt}F_{\epsilon}$ $\displaystyle=$
$\displaystyle-\frac{1}{\epsilon}{\cal{V}}\partial_{z}F_{\epsilon}+\partial_{t}F_{\epsilon}+O(\epsilon)$

with $\partial_{ss}$ the second derivative along $\gamma$. Considering terms
of $O(\epsilon^{-1})$ in eq. (2) leads to
$0=\tilde{\sigma}\left(\partial_{zz}\Phi_{0}-W^{\prime}(\Phi_{0})\right)$.
Using the outer expansion eq. (10) and the matching condition eq. (7) shows
that $\Phi_{0}(z,s,t)$ is a solution of
$\partial_{zz}\Phi_{0}=W^{\prime}(\Phi_{0})$ with $\Phi_{0}(\pm\infty)=\pm 1$
and thus

(13) $\Phi_{0}(z)=\tanh\left(\frac{z}{\sqrt{2}}\right)$

independent of $s$ and $t$. Using this in $O(\epsilon^{0})$ of eq. (2) leads
to

(14)
$-\tilde{\beta}{\cal{V}}G(\Phi_{0})\partial_{z}\Phi_{0}=\tilde{\sigma}\left({\cal{H}}_{\gamma}\partial_{z}\Phi_{0}-W^{\prime\prime}(\Phi_{0})\Phi_{1}+\partial_{zz}\Phi_{1}\right)+G(\Phi_{0})L_{0}.$

Multiplying by $\partial_{z}\Phi_{0}$ and integrating leads to

(15)
$\displaystyle-\tilde{\beta}{\cal{V}}\int_{-\infty}^{+\infty}G(\Phi_{0})(\partial_{z}\Phi_{0})^{2}\;dz$
$\displaystyle=$
$\displaystyle\tilde{\sigma}\int_{-\infty}^{+\infty}{\cal{H}}_{\gamma}(\partial_{z}\Phi_{0})^{2}-\partial_{z}W^{\prime}(\Phi_{0})\Phi_{1}+\partial_{zz}\Phi_{1}\partial_{z}\Phi_{0}\;dz$
$\displaystyle+L_{0}\int_{-\infty}^{+\infty}G(\Phi_{0})\partial_{z}\Phi_{0}\;dz.$

It follows
$\int_{-\infty}^{+\infty}-\partial_{z}W^{\prime}(\Phi_{0})\Phi_{1}+\partial_{zz}\Phi_{1}\partial_{z}\Phi_{0}\;dz=\int_{-\infty}^{+\infty}\Phi_{1}\partial_{z}(-W^{\prime}(\Phi_{0})+\partial_{zz}\Phi_{0})\;dz=0$
and thus

(16)
$-\tilde{\beta}{\cal{V}}\int_{-\infty}^{+\infty}\\!\\!G(\Phi_{0})(\partial_{z}\Phi_{0})^{2}\;dz=\tilde{\sigma}{\cal{H}}_{\gamma}\int_{-\infty}^{+\infty}\\!\\!(\partial_{z}\Phi_{0})^{2}\;dz+L_{0}\int_{-\infty}^{+\infty}\\!\\!G(\Phi_{0})\partial_{z}\Phi_{0}\;dz.$

With eq. (13) we obtain
$\partial_{z}\Phi_{0}=\frac{1}{\sqrt{2}}(1-\Phi_{0}^{2})$ and thus
$\int_{-\infty}^{+\infty}G(\Phi_{0})(\partial_{z}\Phi_{0})^{2}\;dz=\frac{4\sqrt{2}}{5}$,
$\int_{-\infty}^{+\infty}(\partial_{z}\Phi_{0})^{2}\;dz=\frac{2\sqrt{2}}{3}$
and $\int_{-\infty}^{+\infty}G(\Phi_{0})\partial_{z}\Phi_{0}\;dz=2$ and
therefore

(17) $\beta{\cal{V}}=-\sigma{\cal{H}}_{\gamma}-2L_{0}.$

In order to determine $L_{0}$ we consider the constraint
$\frac{1}{|{\cal{S}}|}\int_{\cal{S}}\phi_{\epsilon}\;d{\cal{S}}=\alpha$. Using
eq. (12) in $O(\epsilon^{-1})$ gives
$0=\int_{\gamma}{\cal{V}}\partial_{z}\Phi_{0}\;d\gamma$ and as
$\partial_{z}\Phi_{0}$ is independent of $s$ also
$\int_{\gamma}{\cal{V}}\;d\gamma=0$. Integrating eq. (17) we thus obtain

(18)
$2L_{0}=-\sigma\frac{1}{|\gamma|}\int_{\gamma}{\cal{H}}_{\gamma}\,d\gamma$

which leads to the desired eq. (4). This analysis is not affected by
considering $G_{\eta}$ instead of $G$.

##  3 Graph formulations

Before we numerically solve both models, the degenerate area preserving
surface Allen-Cahn equation (2) and (3) and the area preserving geodesic
curvature flow (4), we reformulate them in a graph formulation. We represent
the closed curve $\gamma(t)$ on ${\cal{S}}$ as a graph of a function
$h:\mathbb{R}^{2}\to\mathbb{R}$, such that

(19) $\gamma(t)=\\{(\mathbf{X},h(\mathbf{X}))^{T}:\mathbf{X}\in g(t)\\}$

where $g(t)$ is a planar curve in $\mathbb{R}^{2}$ and
$\mathbf{X}=\mathbf{X}(l,t)$ its 1-periodic parameter form such that
$g(t)=\\{\mathbf{X}(l,t)=(X_{1}(l,t),X_{2}(l,t))^{T},l\in[0,1]\\}$. In [12]
eq. (4) is analysed by means of the flow of $g(t)$:

(20) $\beta V=-aH_{g}+b+c\frac{1}{\int_{g}\sqrt{1+(\nabla
h\cdot\mathbf{t}_{g})^{2}}\;dg}\int_{g}{\cal{H}}_{\gamma}\sqrt{1+(\nabla
h\cdot\mathbf{t}_{g})^{2}}\;dg$

with normal velocity $V$, coefficients $a>0$, $b$ and $c$ given as

$a=\sigma\frac{1}{1+(\nabla
h\cdot\mathbf{t}_{g})^{2}},\;\;b=\sigma\frac{\mathbf{t}_{g}^{T}\nabla^{2}h\mathbf{t}_{g}(\nabla
h\cdot\mathbf{n}_{g})}{(1+(\nabla h\cdot\mathbf{t}_{g})^{2}(1+|\nabla
h|^{2})},\;\;c=\sigma\sqrt{\frac{1+(\nabla
h\cdot\mathbf{t}_{g})^{2}}{1+|\nabla h|^{2}}}.$

The unit tangent $\mathbf{t}_{g}$, unit outer normal $\mathbf{n}_{g}$ and
curvature $H_{g}$ of the curve $g(t)$ given as

(21)
$\mathbf{t}_{g}=\frac{\partial_{h}\mathbf{X}}{|\partial_{h}\mathbf{X}|},\quad\mathbf{n}_{g}=\frac{1}{|\partial_{h}\mathbf{X}|}(\partial_{h}X_{2},-\partial_{h}X_{1})^{T},\quad
H_{g}=-\frac{1}{|\partial_{h}\mathbf{X}|}\partial_{h}\left(\frac{\partial_{h}\mathbf{X}}{|\partial_{h}\mathbf{X}|}\right)\cdot\mathbf{n}_{g}$

and the geodesic curvature ${\cal{H}}_{\gamma}$ of $\gamma(t)$

(22) ${\cal{H}}_{\gamma}=-\frac{1}{\sqrt[3]{1+(\nabla
h\cdot\mathbf{t}_{g})^{2}}}\left(\sqrt{1+|\nabla
h|^{2}}H_{g}-\frac{\mathbf{t}_{g}^{T}\nabla^{2}h\mathbf{t}_{g}}{\sqrt{1+|\nabla
h|^{2}}}(\nabla h\cdot\mathbf{n}_{g})\right).$

For details and the numerical realization we refer to [12, 14, 15].

The graph formulation for eq. (2) and (3) can be formulated as

(23) $\displaystyle\hskip
28.45274pt\epsilon\tilde{\beta}G_{\eta}(\phi)\partial_{t}\phi=\tilde{\sigma}\left(\epsilon\nabla\cdot\left(\left(I-\frac{(\nabla
h)^{2}}{1+|\nabla
h|^{2}}\right)\nabla\phi\right)-\frac{1}{\epsilon}W^{\prime}(\phi)\right)$
$\displaystyle+\frac{\tilde{\sigma}G_{\eta}(\phi)}{\int_{\Omega}\sqrt{1+|\nabla
h|^{2}}\,d\Omega}\int_{\Omega}\frac{\sqrt{1+|\nabla
h|^{2}}}{G_{\eta}(\phi)}\left(\frac{1}{\epsilon}W^{\prime}(\phi)-\epsilon\nabla\cdot\left(\left(I-\frac{(\nabla
h)^{2}}{1+|\nabla h|^{2}}\right)\nabla\phi\right)\right)\,d\Omega,$

where $\Omega\subset\mathbb{R}^{2}$ and $h(\Omega)={\cal{S}}$. For the
numerical realization we follow [16] and introduce a relaxation rate for the
Lagrange multiplier, which can also be interpreted as an additional
penalization of the area of the form
$c\left(\frac{1}{|\cal{S}|}\int_{\cal{S}}\phi\;d{\cal{S}}-\alpha\right)$, with
penalization parameter $c$. The resulting semi-discrete graph formulation with
$\phi^{n}=\phi(t^{n})$. reads

$\displaystyle\epsilon\tilde{\beta}G_{\eta}(\phi^{n})\frac{\phi^{n+1}-\phi^{n}}{\tau}=\tilde{\sigma}\left(\epsilon\nabla\cdot\left(\left(I-\frac{(\nabla
h)^{2}}{1+|\nabla
h|^{2}}\right)\nabla\phi^{n+1}\right)-\frac{1}{\epsilon}W^{\prime}(\phi^{n+1})\right)$
$\displaystyle+\frac{\tilde{\sigma}G_{\eta}(\phi^{n})}{\int_{\Omega}\sqrt{1+|\nabla
h|^{2}}\,d\Omega}\int_{\Omega}\frac{\sqrt{1+|\nabla
h|^{2}}}{G_{\eta}(\phi^{n})}\left(\frac{1}{\epsilon}W^{\prime}(\phi^{n})-\epsilon\nabla\cdot\left(\left(I-\frac{(\nabla
h)^{2}}{1+|\nabla h|^{2}}\right)\nabla\phi^{n}\right)\right)\,d\Omega$
$\displaystyle-c\left(\frac{\int_{\Omega}\sqrt{1+|\nabla
h|^{2}}\phi^{n}\,d\Omega}{\int_{\Omega}\sqrt{1+|\nabla
h|^{2}}\,d\Omega}-\alpha\right),$

with $W^{\prime}(\phi^{n+1})\approx
W^{\prime}(\phi^{n})+W^{\prime\prime}(\phi^{n})(\phi^{n+1}-\phi^{n})$. The
resulting equation for $\phi^{n+1}$ is linear and discretized in space by
standard $P^{1}$ finite elements. The problem is implemented in AMDiS [17, 18]
and the linear system is solved with the direct solver of UMFPACK.

##  4 Numerical results

We consider the 4 examples provided in [12] and take the numerical solutions
provided in [12] (with $M=200$ finite volumes) as benchmark solutions. The
problem settings are provided in Table 1.

Table 1: Numerical examples, initial parametrization of curve and height profile. $r(l)=1+0.65\cos{(10\pi l})$. |  $\mathbf{X}_{0}$, $l\in[0,1]$ | $h$
---|---|---
Problem 1 | $\mathbf{X}_{0}=(\frac{1}{4}+r(l)\cos{2\pi l},-\frac{1}{4}+r(l)\sin{2\pi l})^{T}$ | $h(x,y)=\sqrt{4-x^{2}-y^{2}}$
Problem 2 | $\mathbf{X}_{0}=(\cos{2\pi l},\frac{1}{10}+\sin{2\pi l})^{T}$ | $h(x,y)=y^{2}$
Problem 3 | $\mathbf{X}_{0}=(\cos{2\pi l},-\frac{1}{5}+\sin{2\pi l})^{T}$ | $h(x,y)=\sin{\pi y}$
Problem 4 | $\mathbf{X}_{0}=(\frac{1}{2}\cos{2\pi l},\sin{2\pi l})^{T}$ | $h(x,y)=x^{2}-y^{4}$

Figure 1 shows selected time instances of the solutions. The spatial
resolution considers at least 10 mesh points across the projected interface in
$\Omega$ and the time step is chosen to ensure the CFL condition with
$\tau\approx k^{2}$, where $k$ is the corresponding mesh size within the
diffuse interface. The mesh is adaptivly refined to ensure these conditions.
Other numerical parameters are chosen as $\eta=0.01$ and $c=2000$. The
physical parameters are set as $\beta=1$ and $\sigma=1$.

Problem 1hxyProblem 2hxyProblem 3hxyProblem 4hxy Figure 1: Comparison of the
numerical solutions of the degenerate area preserving surface Allen-Cahn
equation (solid) with the numerical solutions of the corresponding sharp
interface limit as reference solutions (dotted). Shown is the zero contour of
the phase field function $\phi_{\epsilon}$ for $\epsilon=0.025$. The evolution
in time is shown in color, running from light green to dark blue. The problem
numbers correspond to Table 1 and the time instances shown are Problem 1:
$t=0,0.05,0.1,0.2,0.4,1$, Problem 2: $t=0,0.25,0.5,1,2$, Problem 3:
$t=0,0.5,1,2,4$, and Problem 4: $t=0,0.25,0.5,1,2.4$.

We also measure the space-time error of the Hausdorff distance for different
$\epsilon$ of the projected curves $\gamma_{200}(t)$ (benchmark solutions) and
$\gamma_{\epsilon}(t)$ ($\phi_{\epsilon}=0$ level-sets) onto $\Omega$. We
consider the $L^{2}$ norm in time of this distance. In addition we provide the
Hausdorff distance for the reached equilibrium solutions, see Table 2. As the
benchmark solution is also just a numerical approximation, we only discuss
convergence qualitatively and do not consider any order of convergence. The
values at least indicate a reduction of the considered errors. A more detailed
and analytically supported convergence study of the numerical solutions
requires to extend results of [19] to surfaces, which is beyond the scope of
this paper.

Additionally we analyse the quality of area conservation and the evolution of
the system energy Eq. 1. The results are shown in Fig. 2 and demonstrate the
desired properties.

Table 2: Hausdorff distance of the projected curves $\gamma_{200}(t)$ (benchmark solution) and $\gamma_{\epsilon}(t)$ ($\phi_{\epsilon}=0$ level-set) onto $\Omega$ in space. Shown are the $L^{2}$-norm in time of the Hausdorff distances and the values at the equilibrium state for the four problems shown in Table 1 | Problem 1 | Problem 2 | Problem 3 | Problem 4
---|---|---|---|---
$\epsilon$ |  $L^{2}$-norm | equil |  $L^{2}$-norm | equil |  $L^{2}$-norm | equil |  $L^{2}$-norm | equil
0.1 | 0.0075 | 0.0577 | 0.0014 | 0.0096 | 0.0033 | 0.0211 | 0.0031 | 0.0175
0.05 | 0.0043 | 0.0332 | 0.0009 | 0.0044 | 0.0019 | 0.0105 | 0.0011 | 0.0069
0.025 | 0.0032 | 0.0311 | 0.0007 | 0.0026 | 0.0015 | 0.0074 | 0.0007 | 0.0056
Figure 2: Time evolution of the area enclosed by the zero contour of the phase
field function $\phi_{\epsilon}$ (left) and the system energy Eq. 1 (right)
for $\epsilon=0.025$. The time is normalized for all four problems. In the
left plot we show the reference areas for the four problems according to [12].
For $\epsilon=0.025$, the maximum absolute deviation from the reference area
is $0.082$ for Problem 1, $0.0180$ for Problem 2, $0.0745$ for Problem 3, and
$0.0344$ for Problem 4. The mean absolute deviations are $0.0646$ for Problem
1, $0.0059$ for Problem 2, $0.0633$ for Problem 3, and $0.0133$ for Problem 4.

We would like to remark that the equilibrium shapes in Problem 2 and Problem
3, which are considered on ruled surfaces, are isotropic. Unrolling the
surfaces provide the circular shapes of the final curves. This is consistent
with known results for (reaction-)diffusion problems on curved surfaces. E.g.
for the surface heat equation it is known that the heat kernel to lowest order
only depends on the Gaussian curvature of the underlying surface [20, 21]. As
this is zero for ruled surfaces, the surface should not have any effect on the
evolution and a circular equilibrium shape on the surface, as in flat space,
can be expected.

##  5 Conclusions

We propose a phase field approximation for area preserving geodesic curvature
flow. The considered equation is a degenerate area preserving surface Allen-
Cahn equation. The connection between both models is established by formal
matched asymptotic analysis and confirmed by numerical solutions for different
problems in a graph formulation.

The degeneracy in the surface Allen-Cahn equation results from the de Gennes
factor $G(\phi)$ in the energy (1), see [7]. While the formal matched
asymptotic analysis also holds for $G(\phi)=1$ the numerical results relay on
the de Gennes factor. It ensures $\phi\in[-1,1]$ much better than without it.
This is a desired feature also in 2D and 3D, where the factor $G(\phi)$ is
used in various applications [9, 10]. However, on curved surfaces it is even
more essential as deviations have a more dramatic effect due to the spatially
varying geometric properties of the surface which can enhance the resulting
errors. These geometric properties also need to be considered in the mesh
resolution in the graph formulation. We need to ensure a desired resolution of
the projected diffuse interface. Also the additional penalization of the area
[16] helps to obtain the shown convergence results. If all these aspects are
considered the proposed phase field approximation provides an appropriate way
to solve the highly non-linear problem of area preserving geodesic curvature
flow by standard tools for solving partial differential equations in 2D.

## References

  * [1] P. C. Fife, O. Penrose, Interfacial dynamics for thermodynamically consistent phasefield models with nonconserved order parameter, Electr. J. Diff. Eq. 16 (1995) 1–49.  doi:10877/7579.
  * [2] C. M. Elliott, B. Stinner, A surface phase field model for two-phase biological membranes, SIAM J. Appl. Math. 70 (2010) 2904–2928.  doi:10.1137/090779917.
  * [3] H. Garcke, J. Kampmann, A. Rätz, M. Röger, A coupled surface-Cahn-Hilliard bulk-diffusion system modeling lipid raft formation in cell membranes, Math. Models Meth. Appl. Sci. 26 (2016) 1149–1189.  doi:10.1142/S0218202516500275.
  * [4] A. Rätz, A benchmark for the surface Cahn–Hilliard equation, Appl. Math. Lett. 56 (2016) 65–71.  doi:10.1016/j.aml.2015.12.008.
  * [5] C. M. Elliott, L. Hatcher, B. Stinner, On the sharp interface limit of a phase field model for near spherical two phase biomembranes, Interf. Free Bound. 24 (2022) 263–286.  doi:10.4171/IFB/473.
  * [6] L. Happel, D. Wenzel, A. Voigt, Effects of curvature on epithelial tissue —coordinated rotational movement and other spatiotemporal arrangements, Europ. Phys. Lett. (EPL) 138 (2022) 67002.  doi:10.1209/0295-5075/ac757a.
  * [7] M. Salvalaglio, A. Voigt, S. M. Wise, Doubly degenerate diffuse interface models of surface diffusion, Math. Meth. Appl. Sci. 44 (2021) 5385–5405.  doi:10.1002/mma.7116.
  * [8] S. Dai, J. Renzi, S. M. Wise, Gamma Convergence for the de Gennes-Cahn-Hilliard energy, arXiv (2022) 2210.16492arXiv:2210.16492.
  * [9] M. Naffouti, R. Backofen, M. Salvalaglio, T. Bottein, M. Lodari, A. Voigt, T. David, A. Benkouider, I. Fraj, L. Favre, A. Ronda, I. Berbezier, D. Grosso, M. Abbarchi, M. Bollani, Complex dewetting scenarios of ultrathin silicon films for large-scale nanoarchitectures, Sci. Adv. 3 (11) (2017) eaao1472.  doi:10.1126/sciadv.aao1472.
  * [10] W. B. Andrews, K. L. M. Elder, P. W. Voorhees, K. Thornton, Effect of transport mechanism on the coarsening of bicontinuous structures: A comparison between bulk and surface diffusion, Phys. Rev. Materials 4 (2020) 103401.  doi:10.1103/PhysRevMaterials.4.103401.
  * [11] A. Rätz, A. Ribalta, A. Voigt, Surface evolution of elastically stressed films under deposition by a diffuse interface model, J. Comput. Phys. 214 (2006) 187–208.  doi:10.1016/j.jcp.2005.09.013.
  * [12] M. Kolář, M. Beneš, D. Ševčovič, Area preserving geodesic curvature driven flow of closed curves and a surface, Disc. Contin. Dyn. Sys. B 22 (2017) 3671–3689.  doi:10.3934/dcdsb.2017148.
  * [13] J. Rubenstein, P. Sternberg, Nonlocal reaction—diffusion equations and nucleation, IMA J. Appl. Math. 48 (1992) 249–264.  doi:10.1093/imamat/48.3.249.
  * [14] V. Minárik, M. Beneš, J. Kratochvíl, Simulation of dynamical interaction between dislocations and dipolar loops, J. Appl. Phys. 107 (2010) 061802\.
  * [15] M. Beneš, J. Kratochvíl, J. Křišt’an, V. Minárik, P. Pauš, A parametric simulation method for discrete dislocation dynamics, Europ. Phys. J. ST 177 (2009) 177–192.
  * [16] Q. Du, C. Liu, X. Wang, Simulating the deformation of vesicle membranes under elastic bending energy in three dimensions, J. Comput. Phys. 212 (2006) 757–777.  doi:10.1016/j.jcp.2005.07.020.
  * [17] S. Vey, A. Voigt, AMDiS: adaptive multidimensional simulations, Comput. Vis. Sci. 10 (2007) 57–67.  doi:10.1007/s00791-006-0048-3.
  * [18] T. Witkowski, S. Ling, S. Praetorius, A. Voigt, Software concepts and numerical algorithms for a scalable adaptive parallel finite element method, Adv. Comput. Math. 41 (2015) 1145–1177.  doi:10.1007/s10444-015-9405-4.
  * [19] X. Feng, A. Prohl, Numerical analysis of the Allen-Cahn equation and approximation for mean curvature flows, Num. Math. 94 (2003) 33–65.  doi:10.1007/s00211-002-0413-1.
  * [20] J. H. P. McKean, I. M. Singer, Curvature and the eigenvalues of the Laplacian, J. Diff. Geom. 1 (1967) 43 – 69.  doi:10.4310/jdg/1214427880.
  * [21] J. Faraudo, Diffusion equation on curved surfaces. I. Theory and application to biological membranes, J. Chem. Phys. 116 (2002) 5831–5841.  doi:10.1063/1.1456024.


Biomedical Research and Innovation Institute of Cadiz (INiBICA) Hospital
Universitario Puerta del Mar, Cádiz, Spain 1,2 1,2 2Department of Pediatric
Hematology and Oncology, Hospital de Jerez, Cádiz, Spain 1,2 The edpsmath
documentclass users guide Version 2

# Mathematical modeling of leukemia chemotherapy in bone marrow

Ana Niño-López<EMAIL_ADDRESS>Department of Mathematics, Universidad de
Cádiz, Puerto Real, Cádiz, Spain;  ,  Salvador Chulián  ,  Álvaro Martínez-
Rubio  ,  Cristina Blázquez-Goñi  and  María Rosa

###### Abstract.

Acute Lymphoblastic Leukemia (ALL) accounts for the $80\%$ of leukemias when
coming down to pediatric ages. Survival of these patients has increased by a
considerable amount in recent years. However, around $15-20\%$ of treatments
are unsuccessful. For this reason, it is definitely required to come up with
new strategies to study and select which patients are at higher risk of
relapse. Thus the importance to monitor the amount of leukemic cells to
predict relapses in the first treatment phase.

In this work we develop a mathematical model describing the behavior of ALL,
examining the evolution of a leukemic clone when treatment is applied. In the
study of this model it can be observed how the risk of relapse is connected
with the response in the first treatment phase. This model is able to simulate
cell dynamics without treatment, representing a virtual patient bone marrow
behavior. Furthermore, several parameters are related to treatment dynamics,
therefore proposing a basis for future works regarding childhood ALL survival
improvement.

###### Key words and phrases:

Mathematical model, Leukemia, Treatment, Cancer

###### 1991 Mathematics Subject Classification:

92-10, 34A12

This work has been partially supported by the Fundación Española para la
Ciencia y la Tecnología (FECYT project PR214), the Asociación Pablo Ugarte
(APU, Spain), Junta de Andalucía (Spain) group FQM-201, Ministry of Science
and Technology, Spain (grant number PID2019-110895RB-I00, funded by MCIN/AEI/
10.13039/501100011033). This work was also subsidized by a grant for the
research and biomedical innovation in the health sciences within the framework
of the Integrated Territorial Initiative (ITI) for the province of Cadiz
(grant number ITI-0038-2019). ITI is $80\%$ co-financed by the funds of the
FEDER Operational Program of Andalusia 2014-2020 (Council of Health and
Families).

La leucémie aiguë lymphoblastique (LAL) représente $80\%$ des leucémies
diagnostiquées chez l’enfant. La survie de ces patients a considérablement
augmenté ces dernières années. Malgré cela, environ $15$ à $20\%$ des
traitements échouent. Pour cette raison, il est crucial de réaliser de
nouvelles stratégies d’étude et de sélectionner les patients présentant le
risque de rechute le plus élevé. D’où l’importance de contrôler le nombre de
cellules leucémiques pour prédire les rechutes dans la première phase du
traitement.

Dans ce travail, un modèle mathématique est développé qui décrit le
comportement de la LAL, en étudiant l’évolution du clone leucémique lorsque le
traitement est appliqué. Dans l’étude de ce modèle, nous pouvons observer
comment le risque de rechute est lié à la réponse à la première phase du
traitement. Ce modèle est capable de simuler la dynamique cellulaire sans
traitement, représentant le comportement de la moelle osseuse d’un patient
virtuel. De plus, les différents paramètres liés à la dynamique du traitement
sont étudiés. Par conséquent, la base des travaux futurs visant à améliorer la
survie de la LAL chez l´enfant est proposée.

###### Contents

  1. 1 Material and methods
    1. 1.1 Healthy bone marrow model
    2. 1.2 Modeling leukemia cells
    3. 1.3 Modeling treatment for leukemia
    4. 1.4 Parameters estimation
  2. 2 Results
    1. 2.1 Leukemic and healthy cells dynamics without treatment are properly represented by the model.
    2. 2.2 The first blood extraction in day $+8$ of treatment allows to approximate the influence of prednisone ($\delta_{P}$).
    3. 2.3 Prednisone and Vincristine influence values analysis implies a classification of patients
    4. 2.4 A new model about bone marrow behavior in presence of treatment is proposed to study.
  3. 3 Discussion

## Introduction

Leukemia is a malignant disease originating in the bone marrow. Particularly,
it arises from a disruption in hematopoiesis, the process in charge of blood
cells production [1, 2]. Hematopoiesis is usually depicted as hierarchical
tree, in which a hematopoietic stem cell (HSC) can differentiate into other
cells from the lymphoid or myeloid line. Depending on the branch, platelet,
red blood cells and lymphocytes are produced. Leukemia is not only
distinguished by the linage but also depending on the maturation stage of the
cancer cells. Acute Lymphoblastic Leukemia (ALL) is caused by cells with fast
growth in the lymphoid branch, with special incidence in pediatric patients
[3, 4].

Each type and subtype of leukemia has an associated protocol that specifies
the therapeutic recommendations, and that varies from country to country [5].
In particular, SEHOP-PETHEMA-2013 protocol is used in Spain to treat pediatric
patients with ALL. This protocol was designed by Spanish Society of Pediatric
Hematology and Oncology along with Spanish Hematology Treatment Program and it
is regularly reviewed [6, 7].

There are different treatments and medications depending on the risk group
assigned to a patient: standard, intermediate or high. In such protocols,
patients go through successive phases: induction, consolidation, re-induction
and maintenance. Each phase consists of the combination of chemotherapeutic
agents and corticoids with varying schedules. As reflected in other works [8,
9], the progressive improvement and modification of these protocols has
increased survival rates to around $80$ to $85\%$. Nonetheless, $15$ to $20\%$
of patients fail to achieve long term remission and therefore relapse [10].
Apart from adding to and modifying these protocols, recent reviews have
highlighted the need for new therapies and approaches that would tackle the
disease differently [11, 12]. This also includes the improvement of risk
assignment, the management of secondary effects and the study of personalized
medicine, whose power lies not only in treatment but also in prevention [13].

In this sense, mathematical models have recently found a wide variety of
applications in medicine [14, 15]. In particular, the discipline of
mathematical oncology has been able to approach recent oncological issues [16,
17]. To understand cell dynamics, several works have described biological
processes related to cellular behaviour [18, 19] and also related to
hematopoiesis and lymphopoiesis [20, 21, 22]. This has allowed the
mathematical characterization of the hierarchical structure of blood cell
lineages [23, 24, 25, 26] and their development in bone marrow. Specifically,
ALL has been considered in previous mathematical models [27, 28, 29] with an
structural vision of leukemic cells similar to that of healthy cells. In order
to fulfill their potential, these models should be validated against
appropriate treatment data [30, 31, 32]. Based on the current protocols to
treat ALL patients [33, 34], some mathematical models describe the drugs’
effect regardless of the leukemia type [35, 36], while other focus explicitly
on the myeloid linage [37, 38, 39]. Other mathematical models describe
alternative strategies to chemotherapy, such as CAR T-cell therapy [40, 41,
42, 43], which is successful in patients that do not respond to the standard
treatment.

Inspired by the examples above, we set out in this study to build a
mathematical model that describes the evolution of ALL and its response to
treatment in the bone marrow. This type of leukemia has received less
attention from the mathematical community, and can benefit just as much. The
mathematical model would allow the derivation of theoretical results (like
stability analysis) but also clinical information like the relative importance
of parameters, the estimation of response times or the influence of dosage and
schedule. It could eventually be used to fit individual patient data and
obtain insights into the possibility of personalized therapy for this disease
[5, 28].

In this paper, we focus on the simulation of standard risk ALL pediatric
patients in the first phase of the treatment, or induction phase, which lasts
28 days. In terms of drugs, this phase comprises the corticosteroid
Prednisone, the chemotherapeutic agents Vincristine and Daunorubicine, and
Asparaginase, each with their own dosages and infusion times. In addition,
triple intrathecal is given to prevent central nervous system disease
involvement. The modelling of this treatment stage requires an analysis of the
bone marrow environment in addition to the assumptions related to the growth
of leukemic cells. We study this scenario and simulate how treatment affects
both leukemia and the lymphoid linage in bone marrow. Besides, we propose a
classification of the patients depending on their leukemic cells development
to predict which patients will respond to treatment correctly.

The work has been organised as follows: In Sec. 1 we describe an existing
healthy lymphopoiesis mathematical model and extend it to include leukemia.
This extended model will also include treatment administration. In Sec. 2 we
show model simulations of leukemia growth and its evolution under therapy,
obtaining realistic ranges for parameters and highlighting the most relevant
ones in terms of response. In Sec. 3 we discuss the results obtained.

##  1\. Material and methods

In this section we describe a mathematical model of healthy, homeostatic
lymphopoiesis in the bone marrow. We then include the appearance of a leukemic
cell and study its evolution. We finally approximate a treatment function
following the current treatment scheme during induction phase and model its
effect in the previous models. An overview of all such models is shown in
Figure 1.

Figure 1. Diagram of cell development along with the appearance of leukemia
and its treatment. a)  Immature healthy B cells from the bone marrow can
develop into more mature cells as follows: In an early maturation stage, a
Pro-B cell can either proliferate in its own compartment with a rate
$\rho_{1}$ or progress to the next compartment (Pre-B) with rate $\alpha_{1}$.
Analogously, a Pre-B cell can remain in a Pre-B stage with a proliferation
rate $\rho_{2}$ and develop to the transition compartment with rate
$\alpha_{2}$. Finally we consider that transition cells do not proliferate in
their own compartment and they arrive to the blood flux with rate
$\alpha_{3}$. Growth rate is regulated by the feedback signal $s$. b) Leukemic
cells could come from a Pro-B stage or Pre-B stage, inheriting the
proliferation rate $\rho_{L}$ of the related compartment. A regulatory signal
$s_{L}$ controls leukemic cell development. Furthermore, for the case of ALL
we assume that leukemic cells do not differentiate into a next stage and
therefore $\gamma_{L}$ indicates the blood exit rate for leukemic cells into
blood. c)  When drug is administered, both healthy and leukemic cells are
affected by the treatment, with cells dying at a rate proportional to their
proliferation rate.

###  1.1. Healthy bone marrow model

Following the immunophenotypical characterization of B cells [21], previous
mathematical models, have considered B lymphopoiesis as a sequential process
dividing their maturation stages in Pro-B (early stage characterized by
immunophenotypic markers CD10+/CD45-), Pre-B (intermediate stage with
CD10+/CD45+) and Transition cells (last stage with CD10-/CD45+) [25, 27].
Mathematically, this can be regarded as three compartments respectively
depending on each stage: $C_{1}(t)$, $C_{2}(t)$ and $C_{3}(t)$, along with
their associated parameters related to proliferation and differentiation. This
was translated in [26] into the compartmental model presented below:

$\displaystyle\dfrac{dC_{1}}{dt}$
$\displaystyle=\bar{s}\rho_{1}C_{1}-\alpha_{1}C_{1},$ (1a)
$\displaystyle\dfrac{dC_{2}}{dt}$
$\displaystyle=\bar{s}\rho_{2}C_{2}+\alpha_{1}C_{1}-\alpha_{2}C_{2},$ (1b)
$\displaystyle\dfrac{dC_{3}}{dt}$
$\displaystyle=\alpha_{2}C_{2}-\alpha_{3}C_{3},$ (1c)

where $\rho_{i}$, for $i=1,2$, is the proliferation rate and $\alpha_{i}$,
$i=1,2,3$, is the transition rate related to each compartment $i$.

In this model the authors also took into account a signal in charge of
regulating cell dynamics, $\bar{s}=\bar{s}(t)$, with the following expression:

$\bar{s}=\bar{s}(t)=\dfrac{1}{1+kN},$ (2)

where $k$ is the inhibitory parameter and $N$ the cells subpopulation in
charge of the feedback signalling. According to the results in [26], we
consider here $N=\sum_{i=1}^{3}C_{i}$. This signalling was shown to affect
predominately the proliferation rate due to biological reasons such as bounded
growth or mathematical ones such as positivity of the solutions. Having
reviewed this model, we expand it by including the appearance of a leukemic
clone.

###  1.2. Modeling leukemia cells

Now, in addition to the three compartments $C_{1},C_{2},C_{3}$, we consider a
leukemic cells compartment, $L$. We assume that leukemic cells grow as a
logistic function where the curve’s maximum value is $L_{\max}$ and $\rho_{L}$
is the logistic growth rate. Just like cells in a healthy hematopoiesis, a
leukemic clone is characterized by a maturation stage [23, 27]. Therefore, we
assume $L$ to depend on the stage where the leukemic clone originates, i.e.,
from the Pro-B stage ($\rho_{L}=\rho_{1}$) or from Pre-B stage
($\rho_{L}=\rho_{2}$). Finally, $\gamma_{L}$ would denote the blood exit rate
of leukemic cells. Here we do not consider clones coming from the transition
stage due to the fact that ALL is produced by immature, proliferating cells.

We consider then the following equation to define the flux in the leukemic
cells compartment:

$\dfrac{dL}{dt}=s_{L}\rho_{L}L\left(1-\dfrac{L}{L_{\max}}\right)-\gamma_{L}L.$
(3)

Finally, signalling $s_{L}$ in this compartment has a similar behavior to the
healthy bone marrow compartments as in [26], and we assume it to affect the
proliferation rate as in Model (1). Therefore,

$s_{L}=\bar{s}=\dfrac{1}{1+k\sum_{i=1}^{3}C_{i}},$ (4)

assuming leukemic cells should not be affected by their own signalling due to
the evasion of growth suppressors [44]. Nevertheless, healthy stages are
influenced by cells in compartment $L$, due to the fact that the more cells in
$L$, the less healthy cells develop due to the invasion of the bone marrow.
Consequently, healthy signalling $\bar{s}$ in Model (1) should be modified to:

$s=\dfrac{1}{1+k\left(L+\sum_{i=1}^{3}{C_{i}}\right)}.$ (5)

In addition to these assumptions, we can also include a constant influx of
healthy stem cells $c_{0}$ in our initial compartment $C_{1}$ [40], which
results in the following model, including both leukemic and healthy B cells:

$\displaystyle\dfrac{dC_{1}}{dt}$
$\displaystyle=c_{0}+s\rho_{1}C_{1}-\alpha_{1}C_{1},$ (6a)
$\displaystyle\dfrac{dC_{2}}{dt}$
$\displaystyle=s\rho_{2}C_{2}+\alpha_{1}C_{1}-\alpha_{2}C_{2},$ (6b)
$\displaystyle\dfrac{dC_{3}}{dt}$
$\displaystyle=\alpha_{2}C_{2}-\alpha_{3}C_{3},$ (6c)
$\displaystyle\dfrac{dL}{dt}$
$\displaystyle=s_{L}\rho_{L}L\left(1-\dfrac{L}{L_{\max}}\right)-\gamma_{L}L,$
(6d)

with $s$ as in Eq. (5) and $s_{L}$ as Eq. (4), where $\rho_{L}$ depends on the
origin of the leukemic clone [24, 27].

This model allows us to study different treatment regimes, to compare with
actual results and to include and test what is known about the behaviour and
effect of these drugs [5, 34, 45].

###  1.3. Modeling treatment for leukemia

As explained previously, there are different treatment protocols that can be
applied [28], depending on the drugs used, their doses, or even the time of
infusion. We consider the Induction I’A phase, as described by the SEHOP-
PETHEMA-2013 protocol, for a standard risk patient. This phase consists of
four drugs and lasts 37 days. Triple intrathecal is not included in the study
since it does not affect the lymphoid branch. Each drug has a schedule with
corresponding doses and days of administration, depending on the patient’s
body surface, as indicated in Figure 2.

Figure 2. Schedule for Induction I’A treatment for a Standard Risk Patient.
Data from SEHOP-PETHEMA protocol: Prednisone (PRED) is administered at
$60mg/m^{2}/day$ for the first $28$ days. Then, the dose is reduced to
$30mg/m^{2}/day$ for days $29,30$ and $31$, to $15mg/m^{2}/day$ for days
$32,33$ and $34$ and to $7.5mg/m^{2}/day$ for days $35,36$ and $37$.
Vincristine (VCR): $1.5mg/day$ on days $8,15,22$ and $29$. Daunorubicine
(DNR): $30mg/m^{2}/day$ on days $8$ and $15$. Asparaginase (ASP):
$10000U/m^{2}/day$ on days $12,15,18,21,24,27,30$ and $33$.

To simplify, we consider treatment days as a scale which corresponds leukemia
detection day with day $+0$ and treatment beginning with day $+1$. The most
important days for monitoring the patient are treatment days $+8$, $+15$ and
$+33$. In those days, bone marrow sample extractions are carried out and
lymphocyte levels are studied. In day $+8$ a blood extraction is done to
assess response to prednisone. According to SEHOP-PETHEMA-2013 protocol, there
must be less than $10^{6}$ leukemic cells per milliliters of blood in day
$+8$. The Minimal Residual Disease (MRD) is studied in bone marrow, whose
positivity is defined by the presence of $0.01\%$ or more leukemic cells in
bone marrow [46]. Therefore, patient responds to the treatment if leukemic
cells involve at most $0.01\%$ in bone marrow in day $+15$, and there are no
blasts in day $+33$.

Following studies on the behaviour of each drug [47, 48, 49, 50, 51], we
consider the treatment $\mu_{j}(t)=\mu_{j}$ as a function that describes the
amount of drug $j$ in the bone marrow at time $t$, being $j\in
J=\\{P,V,D,A\\}$: Prednisone $(\mu_{P})$, Vincristine $(\mu_{V})$,
Daunorubicin $(\mu_{D})$ and Asparaginase $(\mu_{A})$. Once the dose is
administered, the drug has an exponential decrease associated with that
medicine half-life.

We then define
$\mu_{j}:\mathbb{R\textsuperscript{+}}\rightarrow\mathbb{R\textsuperscript{+}}$
as

$\dfrac{d\mu_{j}}{dt}=-\lambda_{j}\mu_{j},\\\ $ (7)

with $\lambda_{j}$ related to each drug $j$ half-life $\tau_{j}$ measured in
days:

$\lambda_{j}=\dfrac{\log(2)}{\tau_{j}}.$ (8)

Each drug $j$ has a different dose $q_{j}$ in several days
$\mathcal{D}_{j}\subset\mathbb{N}$, as shown in Figure 2. We define the dose
administered as:

$Q_{j}(t)=\begin{cases}q_{j}&t\in\mathcal{D}_{j},\\\ &\\\
0&t\notin\mathcal{D}_{j}.\\\ \end{cases}$ (9)

We finally consider the total treatment function $\mu$, represented in Figure
3, as the weighted sum of all drugs:

$\mu=\sum_{j\in J}\delta_{j}(\mu_{j}+Q_{j}),$ (10)

where $\delta_{j}$ indicates the influence of drug $j$ on the total effect of
the treatment.

Figure 3. Total treatment influence in bone marrow. Percentage of total drug
dose for Induction phase, following SEHOP-PETHEMA-2013 (Figure 2) and with
parameters from Table 1.

Finally, we include a term in Model (6) that indicates cellular death by
chemotherapy. We consider this does not affect transition cells, owing to the
fact that they do not have any associated proliferation. The resulting model
is shown below (11).

$\displaystyle\dfrac{dC_{1}}{dt}$
$\displaystyle=c_{0}+s\rho_{1}C_{1}-\alpha_{1}C_{1}-\mu\rho_{1}C_{1},$ (11a)
$\displaystyle\dfrac{dC_{2}}{dt}$
$\displaystyle=s\rho_{2}C_{2}+\alpha_{1}C_{1}-\alpha_{2}C_{2}-\mu\rho_{2}C_{2},$
(11b) $\displaystyle\dfrac{dC_{3}}{dt}$
$\displaystyle=\alpha_{2}C_{2}-\alpha_{3}C_{3},$ (11c)
$\displaystyle\dfrac{dL}{dt}$
$\displaystyle=s_{L}\rho_{L}L\left(1-\dfrac{L}{L_{\max}}\right)-\gamma_{L}L-\mu\rho_{L}L,$
(11d) $\displaystyle\mu$ $\displaystyle=\sum_{j\in
J}\delta_{j}\left(\mu_{j}+Q_{j}\right),$ (11e) $\displaystyle Q_{j}$
$\displaystyle=\begin{cases}q_{j}&t\in\mathcal{D}_{j},\\\ &\\\
0&t\notin\mathcal{D}_{j},\\\ \end{cases},$ (11f)
$\displaystyle\dfrac{d\mu_{j}}{dt}$ $\displaystyle=-\lambda_{j}\mu_{j}(t).$
(11g)

Particularly, when $\mu=0$, (i.e. no treatment is applied) Model (11)
transforms into Model (6).

###  1.4. Parameters estimation

We consider literature data for healthy bone marrow properties [26] along with
parameters related to leukemia [40]. Drugs half lives are taken from
pharmacokinetics and pharmacodynamics studies. All parameter values are
included in Table 1. The estimations for dosage assume a $25-30$ $kg$ body
mass for a child, which in turn implies a body surface area equivalent to
$1\;m^{2}$ [52].

Parameter Meaning Value Unit Source $c_{0}$ Influx cells from HSC $10^{7}$
$cell$ Estimated from [40] $\rho_{1}$ Pro-B proliferation rate $\log(2)$
$day^{-1}$ [26] $\rho_{2}$ Pre-B proliferation rate $\log(2)/1.5$ $day^{-1}$
[26] $\alpha_{1}$ Transition rate: Pro-B to Pre-B $0.168$ $day^{-1}$ [26]
$\alpha_{2}$ Transition rate: Pre-B to Transition $0.144$ $day^{-1}$ [26]
$\alpha_{3}$ Blood exit rate $0.288$ $day^{-1}$ [26] $k$ Signal Intensity
$10^{-10}$ $cell^{-1}$ [40] $\gamma_{L}$ Leukemic cells blood exit rate
$0.288\times 10^{-3}$ $day^{-1}$ [40] $L_{\max}$ Leukemic cells carrying
capacity $10^{12}$ $cell$ [40] $q_{P}$ Prednisone dose $60$ $mg/day$ SEHOP-
PETHEMA-2013 $q_{V}$ Vincristine dose $1.5$ $mg/day$ SEHOP-PETHEMA-2013
$q_{D}$ Daunorubicin dose $30$ $mg/day$ SEHOP-PETHEMA-2013 $q_{A}$
Asparaginase dose $10000$ $U/day$ SEHOP-PETHEMA-2013 $\delta_{P}$ Prednisone
influence $\left[\dfrac{1}{60},\dfrac{1}{6}\right]$ $day/mg$ Estimated
$\delta_{V}$ Vincristine influence
$\left[\dfrac{1}{1.5},\dfrac{10}{1.5}\right]$ $day/mg$ Estimated $\delta_{D}$
Daunorubicin influence $\left[\dfrac{1}{30},\dfrac{1}{3}\right]$ $day/mg$
Estimated $\delta_{A}$ Asparaginase influence
$\left[\dfrac{1}{10^{4}},\dfrac{1}{10^{3}}\right]$ $day/U$ Estimated
$\lambda_{P}$ Prednisone decaying rate $9.6\log(2)$ $day^{-1}$ [47]
$\lambda_{V}$ Vincristine decaying rate $0.28\log(2)$ $day^{-1}$ [48]
$\lambda_{D}$ Daunorubicin decaying rate $1.17\log(2)$ $day^{-1}$ [49]
$\lambda_{A}$ Asparaginase decaying rate $0.8\log(2)$ $day^{-1}$ [50]
$C_{1}(0)$ Pro-B cells $3.52211\times 10^{9}$ $cell$ [26] $C_{2}(0)$ Pre-B
cells $1.84911\times 10^{10}$ $cell$ [26] $C_{3}(0)$ Transition cells
$9.24555\times 10^{9}$ $cell$ [26] $L(0)$ Leukemic cells $1$ $cell$ Assumption

Table 1. Parameter values. $L_{\max}$ is measured in number of cells and $k$
in $cell^{-1}$ while the rest of the presented parameters related to the
proliferation and transition rates are measured in $day^{-1}$. On the other
hand, medicine dosage $q_{j}$ is measured in $mg/day$ and the influence of
each drug is considered $day/mg$. Asparaginase is a particular case since is
measured in $\mathcal{U}/day$ and its influence in $day/\mathcal{U}$. Values
of $\lambda_{j}$ are measured in $day^{-1}$ due to the fact that they come
from half-life values. Values used to initialize each case $C_{i}(0)$,
$i=1,2,3$, have been obtained from the steady values in healthy bone marrow
models [26] along with our parameters. A single leukemic cell is supposed at
the beginning, from which Model (6) starts, and the influx of cells $c_{0}$ is
also obtained from Model (1) steady states.

The study of all estimated parameters is presented in Supplementary
Information document. Control parameters related to drug, $\delta_{j}$ with
$j\in J$, have been selected by searching for values that implies that the
patient responds to treatment correctly, i.e. considering an amount of blasts
less than $10^{6}$ $blasts/ml$ in blood in day $+8$ and MRD$<0.01\%$ in bone
marrow in day $+15$ [53]. In preliminary simulations we find that if
$\delta_{j}q_{j}<10^{0}$, the dose administered is ineffective. We further
consider that $\delta_{j}q_{j}>10^{1}$ involves the death of the patient due
to drug toxicity. We then set all $\delta_{j}q_{j}$ to the same order of
magnitude $\left[10^{0},10^{1}\right]$, $\forall j\in J$. In this way we
ensure Model (11) can represent a virtual patient whose behavior is similar to
real patients according to the amount of cells in each bone marrow aspiration.

##  2\. Results

###  2.1. Leukemic and healthy cells dynamics without treatment are properly
represented by the model.

We now consider no treatment is applied, Model (6), i.e., Model (11) with
$\mu_{j}=0$, $\forall j\in J$. This model combines the healthy dynamics of
immature B-cells and how leukemic cells stability values as initial inputs for
the model. We simulate the appearance of a leukemic cell (at $t=0$), depending
on the leukemic clone originates, either in a Pro-B
$\left(\rho_{L}=\rho_{1}\right)$ or a Pre-B stage
$\left(\rho_{L}=\rho_{2}\right)$. Results of solving in a range to
$\left[0,300\right]$ days are shown in Figure 4, displaying the range of days
with more information in each case.

We next show the simulations for each case. Figure 4 represents Model (6)
cells evolution when $\rho_{L}=\rho_{1}$ is considered for the first case, and
secondly, model for $\rho_{L}=\rho_{2}$ is simulated. In addition, for each
time $t$ measured in days, the bone marrow cells proportions are shown for the
different maturation compartments.

Figure 4. Dynamics of cell compartments proportion in the bone marrow in the
presence of the leukemic clone. Initial data and parameters from Table 1. (a)
Case $\rho_{L}=\rho_{1}=0.6931\ day^{-1}$. Leukemic cells, $L$, (solid purple
line) have a logistic growth and reach $80\%$ around day $150$. At the same
time, Pro-B cells ($C_{1}$, solid blue line), Pre-B cells ($C_{2}$, dash-
dotted blue line) and Transition cells ($C_{3}$, dotted blue line) decrease.
(b) Case $\rho_{L}=\rho_{2}=0.4621\ day^{-1}$. Dynamics is similar to the
first case but leukemia is detected around the day 223.

As it can be observed in Figure 4, both graphics represent the total bone
marrow invasion by leukemia. Depending on the stage from which the first
leukemic cell comes, leukemic cells proliferate more or less quickly. The
appearance of a leukemic clone is set to be detected when the leukemic
population accounts for $80\%$ of blasts from the total population of B
lymphocytes, [54]. When the first leukemic cell is originated in Pro-B stage,
$\rho_{L}=\rho_{1}$, the disease is detected around day $150$. However, if the
original leukemic cell comes from the Pre-B stage, leukemia is appreciated
from day $223$. Both two cases imply a progressive decreasing in healthy cells
due to the fact that the invasive ability of leukemic cells prevents the
normal development of cells in the bone marrow.

###  2.2. The first blood extraction in day $+8$ of treatment allows to
approximate the influence of prednisone ($\delta_{P}$).

Treatment begins according to treatment schedule in SEHOP-PETHEMA-2013
protocol, after $80\%$ of blasts are detected, producing a drop in both
leukemic and healthy cells.

In order to find the optimal parameters to represent Model (11), we pay
attention to the first blood extraction (day $+8$). To do so, we estimate
relations between blasts in bone marrow and blood, owing to the fact that
there is a correlation between them [55, 56].

SEHOP-PETHEMA-2013 protocol defines a good response to prednisone if there is
fewer than $10^{6}$ blasts per milliliter of blood. We assume a pediatric
patient with $25-30kg$ body mass, and due to the blood - body mass correlation
[57, 58], this assumption entails approximately $2.1-2.5$ liters of blood. We
consider a patient with $2.3l$ of blood to our study, and consequently,
patients should have less than $2.3\times 10^{9}$ blasts in blood. It implies
that leukemic cells limit in the bone marrow should be less than $2.3\times
10^{10}$ blasts at day $+8$ to obtain a good response to prednisone. We
analyze all $\delta_{P}$ possible values within our range in Table 1, taking
into account leukemic cells evolution for the first eight days of treatment
(see Figure 5).

Figure 5. Variation of $\delta_{P}$ values to analyze results in day $+8$ of
treatment. As presented in Table 1,
$\delta_{P}\in\left[1/60,1/6\right]=\left[\delta_{P_{\min}},\delta_{P_{\max}}\right]$.
Five values of $\delta_{P}$ are highlighted: the midpoint between the maximum
and minimum values $\delta_{P_{0}}\approx 0.09$, $\delta_{P_{1}}\approx 0.05$
the middle between $\delta_{P_{\min}}$ and $\delta_{P_{0}}$, and
$\delta_{P_{2}}\approx 0.13$ the center of the interval
$\left[\delta_{P_{0}},\delta_{P_{\max}}\right]$ . (a) Model (11) is solved for
$50$ values of prednisone influence. For each $\delta_{P}$ value, amount of
blasts in day +8 is obtained and associated with a patient who responses
(green points) or not (gray points) to prednisone. (b) Leukemic cells
evolution depending on several $\delta_{P}$ values. Green-dotted lines denote
patients who response to prednisone, while gray-dashed lines represent bone
marrow dynamics which do not reach any response to prednisone from data in day
$+8$ (yellow vertical line). Horizontal solid red line represents the limit of
blasts assumed by SEHOP-PETHEMA-2013 protocol to consider a good response of
less than $2.3\times 10^{10}$ blasts in bone marrow.

As it is shown in Figure 5 (a), the amount of blasts in bone marrow in day
$+8$ of treatment implies a direct consequence related to possible
$\delta_{P}$ values. Values lower than $\delta_{P_{0}}$ imply that patients do
not response to prednisone, while values higher than $\delta_{P_{0}}$ entail a
good prognosis for the patient due to the fact that in day $+8$ of treatment
reach less than the limit of blasts assumed.

On the other hand, in Figure 5 (b) the development of leukemic cells according
to Model (11) is depicted for five different values of $\delta_{P}$. Values
which allow to reach blasts levels lower than $2.3\times 10^{10}$ total blasts
in day $+8$ corresponds to the green area in Figure 5 (a).

###  2.3. Prednisone and Vincristine influence values analysis implies a
classification of patients

Let us how consider all parameters $\delta_{j}$ in order to analyze their
behaviour in Model (11), taking into account the number of blasts in bone
marrow at day $+15$.

Based on the parameter sensitivity analysis (Supplementary Information), main
influential parameter is $\delta_{P}$. To study their influence, we vary
parameters $\delta_{P}$ and $\delta_{V}$, related, respectively, to prednisone
and vincristine, and check how many leukemic cells are there day $+15$ of
treatment schedule in comparison to the healthy population (MRD). We fix
parameters $\delta_{D}=2.5/30day/mg$ and $\delta_{A}=2.5/10^{4}day/U$ given
their minimal influence in the sensitivity analysis. We focus on the variation
of $\delta_{P}$ values in comparison with $\delta_{V}$ values to obtain the
minimal value which implies the patient responds to treatment. We assume
patient responds to treatment correctly if day $+15$ of treatment leukemic
cells represent less than $0.01\%$ of the whole lymphocyte population.

Values ranges considered for $\delta_{j}$, $j\in J$ are those in Table 1. In
Supplementary Information, we show a heat map which takes parameters in the
widest range. That range of parameter for such heat map has been adapted to a
more meaningful range. We show in Figure 6 a heat map which connects
parameters variation with the blasts percentage in day $+15$.

Figure 6. Blasts percentage in bone marrow day $+15$ of treatment. We consider
$21$ values up to $0.167day/mg$ for Prednisone and $21$ values up to
$4.22day/mg$ for Vincristine. For each values combination, Model (11) is
solved for parameters in Table 1 with $\delta_{D}=2.5/30day/mg$,
$\delta_{A}=2.5/10^{4}day/U$. The value in each box corresponds to the
percentage of leukemic cells in day $+15$ depending on values $\delta_{P}$ and
$\delta_{V}$.

This heat map leads to associate all values with two semiplanes which will be
able to predict if the patient responds correctly depending on the values
which match his data. Regarding Figure 6, we can approximate the influence of
each drug in a patient and predict if the patient responds to the treatment
correctly or not. Firstly, parameters related to the influence of drugs in
cells death, $\delta_{j}$, can be approximated for a set of data from a
specific patient in several days of evolution of the treatment. Thus, we
classify these parameters according to Figure 6 in purple or white zone to
predict which percentage of leukemic cells will be there in bone marrow day
$+15$. There will be a good prognosis if leukemic compartment entails less
than $0.01\%$ (MDR) in bone marrow, white zone, while if those values are in
the purple zone, the patient will not respond correctly to treatment.

According to choose the minimal values which implies a good prognosis of the
patient, $\delta_{P}=0.092day/mg$ and $\delta_{V}=2.11day/mg$ will be
parameters used in the next results.

###  2.4. A new model about bone marrow behavior in presence of treatment is
proposed to study.

After proposing Model (11), we simulate the situation using data in Table 1
along with parameters associated to the $\delta_{j}$ values study in previous
sections. As proposed in Sec. 1.2 it is possible to observe that leukemic
cells have a higher proliferation than the rest of lymphocytes. Moreover, as
said before, it is known that the blasts level in bone marrow when leukemia is
detected is around $80\%$. We assume that treatment starts the next day and
present results in Figure 7.

Therefore, according to the available data about different patients, we have
been able to approximate the evolution of bone marrow cells in the presence of
treatment as it is shown in Figure 7. It can be observed that the treatment
affects both healthy and leukemic cells and that once the disease is
eliminated, healthy cells compartments return to normal levels.

Figure 7. Dynamics of cell compartments in the bone marrow in the presence of
the leukemic clone along with therapy. Initial data, parameters related to
leukemia and parameters about treatment are those presented in Table 1 except
for $\delta_{P}=0.092day/mg$, $\delta_{V}=2.11day/mg$,
$\delta_{D}=2.5/30day/mg$, $\delta_{A}=2.5/10^{4}day/U$ and
$\rho_{L}=\rho_{1}=0.6931\ day^{-1}$ case which we consider without loss of
generality. (a) Cells proportion evolution in bone marrow from the first
leukemic cell appearance along with the treatment application. (b) Bone marrow
dynamics centered between day $150$ and day $180$, corresponding to $+0$ and
$+30$ from therapy, respectively. (c) Compartment values measured in amount of
cells in ten first days of treatment. It is observed that day $+8$ is the most
aggressive of them. (d) Treatment finishes in day $+37$. The extinction of
leukemic cells implies healthy cells back to normal.

Generally, in Figure 7, a leukemic bone marrow with treatment is represented.
From usual cells compartments proportions, leukemia develops until leukemic
cells exceed $80\%$ of total occupancy in bone marrow. From those conditions,
the therapy is applied for $37$ days, and then, bone marrow normal levels are
recovered. In a reduced range of days (Figure 7(b)) it can be appreciated how
cells proportions are modified according to drugs administered each day. In
fact, day $+8$ is the first aggressive day in which leukemic cells do not grow
so fast, therefore, in following days this population continues its reduction.
It is for that reason that the highest cells proportion in bone marrow for
about $30$ days of treatment is the transition compartment, $C_{3}$, which are
not influenced by treatment since they do not proliferate. Paying attention to
the amount of cells (Figure 7(c)-(d)), it is possible to verify that all of
the cells in bone marrow decrease gradually and finally, return to normal
levels.

##  3\. Discussion

Leukemia is the cancer with the highest incidence in pediatric age,
particularly, Acute Lymphoblastic Leukemia and $20\%$ of these diagnoses do
not respond correctly to treatment, and as a consequence, there is a relapse.
This problem is the main reason of our study and the principal question to
answer is why the treatment does not work for those patients. Our intention is
to find patterns in bone marrow along with medicines applied, so that, it is
possible to predict relapses and improve current treatments. We approximate
treatment influence in cell death to analyze the behavior of each bone marrow
based on flow cytometry data throughout treatment time. Modeling current
treatment in leukemia B leads to advances towards optimizing the application
of medicines. After studying leukemic cells behavior it is possible to analyze
and advance in this area. Available clinical data are important to determinate
parameters values according to different bone marrow samples from the same
patient.

First of all, we have reviewed a healthy bone marrow model in Sec. 1.1. That
study allows us to include the assumption of the appearance of a leukemic
cell. It is found that leukemic cells have a strong development due to their
properties of malignant disease. The leukemic cell develops as it has been
described in Model (6), and it could be more or less accelerated depending on
the inherited proliferation rate for $\rho_{L}$, as shown in Figure 4. In
fact, available data [30, 54] of hematologic patients reinforce that the
disease is diagnosed when there are about $80\%$ leukemic cells in bone
marrow. As it is mentioned in Sec. 2.1, six steady states are obtained by
studying Model (6). The unique stable steady state is the one that represents
the total occupancy of the bone marrow by leukemic cells. The rest of
solutions are unstable. There are several possible explanations for this
result. On one hand, both two solutions which have negative compartments, are
biologically senseless results. In spite of the fact that there are two more
solutions which are non negative, they lack biological meaning. Moreover,
there is a solution which reproduces healthy bone marrow behavior without
leukemia, Sec. 1.1. Therefore, Model (6) replicates reality and its study
allows to adjust parameters related to leukemic population growth. In the case
in which leukemia appears in the Pro-B stage, $\rho_{L}=\rho_{1}$ for Model
(6), the amount of leukemic cells is the same as Pro-B cells around 127 days;
it is equal to Transition cells in the day 132; and Pre-B cells are reached by
leukemic cells in 135 days. Simulations for model which assume the appearance
of the leukemic cell in Pre-B stage, $\rho_{L}=\rho_{2}$, have a similar
behavior, but the process is slower. In general, both cases have the same
behavior. Every simulation in the presence of the leukemic clone allow for
deducing that healthy and leukemic cells coexist for a period of time which
depends on the case, and then, leukemic cells take up all the bone marrow
capacity, therefore, there is no space for healthy cells. According to the
Model (6) and its simulations, it is important to take into account that cells
development is supposed without any treatment which is responsible for
reducing the leukemic cells proportion. Hence, it explains why there are that
suddenly changes three months later.

From the information collected in the SEHOP-PETHEMA-2013 protocol, we model
the treatment behavior in the bone marrow when it is administered. Without
loss of generality, we consider a leukemic cell originated in Pro-B
compartment. From the first leukemic cell appearance, the proliferation of
that population takes about 150 days to be diagnosed with the result that when
it is detected, the disease reproduces more quickly. That is why treatment
begins to be administered in the following days, thus, there are no data
related to the same patient in different times without treatment. We have
focused on a standard risk patient in Induction I’A which is the first phase
of the treatment and it is a determining phase in the therapy. If there are
signs which do not correspond to proper treatment response, the protocol for
this patient is modified and associated to other risk level. Therefore, we
consider a therapy based on Prednisone, Vincristine, Daunorubicine and
Asparaginase, administered according to the schedule shown in Figure 2. Taking
into account the total amount of drugs in the body during the treatment, we
assume an exponential decay of each drug from its application in order to
different references about timing and dosage drugs effects [28, 59, 60]. In
addition, we include the term $\delta_{j}$ as the influence of each drug in
cells death, to obtain the function $\mu$ which imply leukemic and healthy
cells death by treatment depending on their proliferation.

When a sample from a patient is studied as a possible diagnosis of ALL is
owing to the patient has previously developed some of the hallmarks of cancer
[44]. Generally, bone marrow first extractions indicate that leukemia takes up
around $80\%$ of bone marrow. Therefore, protocol is immediately activated.
That treatment approximation implies an improvement about the perception
related to drugs actions in the body. Depending on each patient, the influence
of each drug will be more or less aggressive, hence, the importance of the
personalized medicine [13]. It is for that reason that a detailed study about
drugs influences must be contemplated to study those behaviors. The aim is to
provide each new patient with the best treatment. A patient is considered with
good prognosis if some conditions are reached. On the treatment day $+8$,
blood sample should show less than $10^{6}$ blasts per milliliters of blood,
on the day $+15$, leukemic cells must entail values smaller than $0.01\%$ and
the day $+33$ blasts in bone marrow must not be observed, considering that
this is the moment of the transplant.

A theoretical model of the bone marrow dynamics, Model (11), is presented to
study, based on previous works along with new assumptions about leukemic cells
behaviour and treatment. In order to study parameters related to prednisone,
we have focused on the first eight days of treatment. In those days only
prednisone is administered, therefore, it is possible to approximate optimal
$\delta_{P}$ value which implies the good prognosis condition in day $+8$.
From range of $\delta_{P}$ possible values in Table 1, Model (11) is solved
for 50 values in that range, obtaining amount of blasts in blood in day $+8$
related to each one. It is obtained a value $\delta_{P_{0}}\approx 0.09$ from
which, values imply blasts levels smaller than the limit of amount of cells,
consequently, we obtain a range of $\delta_{P}$ values in which patient
response to treatment. This result imply two important aspects to consider. On
one hand, a patient whose sample is extracted in day $+8$ could be associated
to the corresponding $\delta_{P}$ and predict the evolution of the disease. On
the other hand, if we analyze blood extraction in other day, $+5$ for
instance, we could be able to associated $\delta_{P}$ value before day $+8$
and predict results in that day. Thus, it would be possible to change
treatment protocol on time.

Apart from the study of the day $+8$ conditions, we analyze good response
conditions in day $+15$. According to sensitivity study, prednisone and
vincristine are the most significant parameters and hence, we propose a study
about influence parameters related to them. On day $+15$, MRD should be less
than $0.01\%$ to obtain good results in patients, therefore, we search for
values which lead to those leukemic cells levels when Model (11) is solved for
those values. Leukemic compartment percentages in bone marrow on day $+15$ are
obtained for each case and presented in Figure 6 in a heat map. As it is
observed, there are two semiplanes corresponding to patients who response and
patients who do not response correctly to treatment. With the purpose of
minimizing both $\delta_{P}$ and $\delta_{V}$, we take average minimal values
which imply a good response. What is interesting in this result is the fact
that the $\delta_{P}$ value obtained by this method is the same value found by
the study of the blasts levels in day $+8$. This finding reinforces that
studies must be based on fitting influence parameters to obtain more
predictive information of each patient.

Finally, simulations supported by previous results in this work are presented.
Previously, parameter estimation is performed and the bone marrow behavior is
studied assuming a leukemic cell appearance. The disease is diagnosed when it
entails an important part of cells in bone marrow and it is treated according
to protocol. As it is observed in the Figure 7 and generally, treatment
application leads to patients recovery and the stabilization of healthy cells
levels from fifty days after treatment finish. Treatment is applied when
leukemia is diagnosed, that it is supposed in day $150$ from the first
leukemic cell appearance since there are about $80\%$ of blasts in bone
marrow. During eight days prednisone is applied and it is possible to observe
a drop in all leukemic and healthy cells but in day $+8$ is more significant
because vincristine and daunorubicin are added to the therapy. From day $+12$
in which asparaginase is administered, leukemic cells values are insignificant
in comparison to healthy cells. The fact that this situation occurs agree with
parameters choice which guarantee that MRD is less than $0.01\%$ on day $+15$.
Graphs interpretation is facilitated by the comparison between amount of cells
and proportions. In those simulations, it should be highlighted that leukemic
cells do not appear the following $150$ days. For that reason, it is
considered a partial patient recovery, not completely due to the logistic
growth of leukemia. Moreover, the disease always tend to appear again and that
idea explains why treatment should continue with other phases which drop the
possibility of any rest of leukemic cells develops.

Being limited to available data, this study lacks accuracy related to drug
influence estimated parameters. On one hand, an exhaustive study of each
patient would required to be able to model the leukemia growth without
treatment. To do so, two bone marrow must be extracted, to observe how cells
develop and obtain the parameter related to leukemic proliferation,
$\rho_{L}$, in each patient. This option is unfeasible since two bone marrow
extractions are rare and unnecessary, as drugs are to be administered to the
child as soon as possible. Hence, we propose a future study to relate blood
blasts and search for a possible correlation. Another limitation of our
analysis would be the simplicity of the treatment function. Drugs evolution
inside the body has a complicated and more detailed process which has been
reduced and it produces a lack of information. Absorption time or action
mechanism are some of the conditions to be taken into account in order to get
closer to real data. Moreover, only the Induction A phase has been formulated
instead of the full treatment. We have focused on this phase of treatment
since effects of the next stages (Induction B and Consolidation) are not
visible enough. In addition, Induction A is the most interesting phase to
study owing to its critical influence throughout treatment.

Notwithstanding with these limitations, the study suggests that it is possible
to approximate contemplated parameters related to treatment provided that
parameters values for leukemia are corroborated. Additionally, different
compartment levels in patients bone marrows with ALL corresponds to values
obtained from this previous study in which a virtual patient is simulated.
These findings suggest several courses of action for improving techniques used
by protocols. There is a definite need for personalizing medicine. We suggest,
therefore, that studies consider mechanisms of action of drugs along with an
exhaustive monitoring to prove how they inhibit healthy and leukemic cells
proliferation as well as how they destroy cells and the environment. We
propose a clinical study about model with those conditions. Comparisons
between models including new assumptions will be able to offer some
approximations to specific models capable of obtaining realistic parameters.
Several patients data are essential to validate the results presented. The
effectiveness of that resource will allow to predict relapses in patients who
have been diagnosed with ALL in a short period of time after treatment begins.
The challenge now is to model full treatment, which takes into account
available data and suitable parameters which lead to new research to improve
current protocols and, consequently, survival expectancy in childhood ALL.

In conclusion, in this paper, a mathematical model of the appearance of
leukemic cells in a bone marrow has been presented based on a mathematical
model of cells development in a healthy bone marrow. Therefore, we have
studied leukemic cells behavior to characterise a model describing the
lymphopoiesis process when a cancer cell appears.Furthermore, we have
simulated the model proposed and obtained results which coincide with
biological basis of cancer. Afterwards, treatment application is included
according to current protocol, meaning healthy and leukemic cells death.
Dynamics in bone marrow are modeled and studied based on parameters estimation
to simulate a virtual patient who responds to treatment.This model provides a
basis of future action protocols of treatments of ALL.

## References

  * [1] S. H. Orkin, L. I. Zon, Hematopoiesis: an evolving paradigm for stem cell biology, Cell 132 (4) (2008) 631–644.
  * [2] M. Jagannathan-Bogdan, L. I. Zon, Hematopoiesis, Development 140 (12) (2013) 2463–2467.
  * [3] C.-H. Pui, D. Campana, W. E. Evans, Childhood acute lymphoblastic leukaemia–current status and future perspectives, The lancet oncology 2 (10) (2001) 597–607.
  * [4] C.-H. Pui, L. L. Robison, A. T. Look, Acute lymphoblastic leukaemia, The Lancet 371 (9617) (2008).
  * [5] R. A. Egler, S. P. Ahuja, Y. Matloub, L-asparaginase in the treatment of patients with acute lymphoblastic leukemia, Journal of pharmacology & pharmacotherapeutics 7 (2) (2016) 62.
  * [6] M. Mesegué, A. Alonso-Saladrigues, S. Pérez-Jaume, A. Comes-Escoda, J. L. Dapena, A. Faura, N. Conde, A. Català, A. Ruiz-Llobet, E. Zapico-Muñiz, et al., Lower incidence of clinical allergy with peg-asparaginase upfront versus the sequential use of native e. coli asparaginase followed by peg-asp in pediatric patients with acute lymphoblastic leukemia, Hematological Oncology 39 (5) (2021) 687–696.
  * [7] A. Ruiz-Llobet, S. Gassiot, E. Sarrate, J. Zubicaray, J. L. Dapena, S. Rives, J. Sevilla, Á. Menárguez López, M. Panesso Romero, C. Montoya, et al., Venous thromboembolism in pediatric patients with acute lymphoblastic leukemia under chemotherapy treatment. risk factors and usefulness of thromboprophylaxis. results of lal-sehop-pethema-2013, Journal of Thrombosis and Haemostasis (2022).
  * [8] S. P. Hunger, X. Lu, M. Devidas, B. M. Camitta, P. S. Gaynon, N. J. Winick, G. H. Reaman, W. L. Carroll, Improved survival for children and adolescents with acute lymphoblastic leukemia between 1990 and 2005: a report from the children’s oncology group, Journal of clinical oncology 30 (14) (2012) 1663.
  * [9] H. Ma, H. Sun, X. Sun, Survival improvement by decade of patients aged 0–14 years with acute lymphoblastic leukemia: a seer analysis, Scientific reports 4 (1) (2014) 1–7.
  * [10] Z. J. Ward, J. M. Yeh, N. Bhakta, A. L. Frazier, R. Atun, Estimating the total incidence of global childhood cancer: a simulation-based analysis, The Lancet Oncology 20 (4) (2019) 483–493.
  * [11] D. Bhojwani, C.-H. Pui, Relapsed childhood acute lymphoblastic leukaemia, The Lancet Oncology 14 (6) (2013) e205–e217.
  * [12] T. Terwilliger, M. Abdul-Hay, Acute lymphoblastic leukemia: a comprehensive review and 2017 update, Blood Cancer Journal 7 (6) (2017) e577–e577.
  * [13] S. Mathur, J. Sutton, Personalized medicine could transform healthcare, Biomedical reports 7 (1) (2017) 3–5.
  * [14] G. Bocharov, V. Volpert, B. Ludewig, A. Meyerhans, et al., Mathematical immunology of virus infections, Vol. 245, Springer, 2018.
  * [15] A. Tosenberger, F. Ataullakhanov, N. Bessonov, M. Panteleev, A. Tokarev, V. Volpert, Modelling of thrombus growth in flow with a dpd-pde method, Journal of theoretical biology 337 (2013) 30–41.
  * [16] R. A. Gatenby, P. K. Maini, Mathematical oncology: cancer summed up, Nature 421 (6921) (2003) 321–321.
  * [17] P. M. Altrock, L. L. Liu, F. Michor, The mathematics of cancer: integrating quantitative models, Nature Reviews Cancer 15 (12) (2015) 730.
  * [18] M. C. Mackey, R. Rudnicki, Global stability in a delayed partial differential equation describing cellular replication, Journal of Mathematical Biology 33 (1) (1994) 89–109.
  * [19] N. L. Komarova, Principles of regulation of self-renewing cell lineages, PloS one 8 (9) (2013) e72847.
  * [20] N. Bessonov, L. Pujo-Menjouet, V. Volpert, Cell modelling of hematopoiesis, Mathematical Modelling of Natural Phenomena 1 (2) (2006) 81–103.
  * [21] E. G. Van Lochem, V. H. J. Van der Velden, H. K. Wind, J. G. Te Marvelde, N. A. C. Westerdaal, J. J. M. Van Dongen, Immunophenotypic differentiation patterns of normal hematopoiesis in human bone marrow: Reference patterns for age-related changes and disease-induced shifts, Cytometry Part B: Clinical Cytometry 60 (1) (2004) 1–13.
  * [22] T. Stiehl, A. Marciniak-Czochra, Characterization of stem cells using mathematical models of multistage cell lineages, Mathematical and Computer Modelling 53 (7-8) (2011) 1505–1517.
  * [23] D. Bonnet, J. E. Dick, Human acute myeloid leukemia is organized as a hierarchy that originates from a primitive hematopoietic cell, Nature medicine 3 (7) (1997) 730–737.
  * [24] K. Anderson, C. Lutz, F. W. Van Delft, C. M. Bateman, Y. Guo, S. M. Colman, H. Kempski, A. V. Moorman, I. Titley, J. Swansbury, et al., Genetic variegation of clonal architecture and propagating cells in leukaemia, Nature 469 (7330) (2011) 356–361.
  * [25] A. Marciniak-Czochra, T. Stiehl, A. D. Ho, W. Jäger, W. Wagner, Modeling of asymmetric cell division in hematopoietic stem cells regulation of self-renewal is essential for efficient repopulation, Stem cells and development 18 (3) (2009) 377–386.
  * [26] S. Chulián, A. Martínez-Rubio, A. Marciniak-Czochra, T. Stiehl, C. B. Goñi, J. F. R. Gutiérrez, M. R. Orellana, A. C. Robleda, V. M. Pérez-García, M. Rosa, Dynamical properties of feedback signalling in b lymphopoiesis: A mathematical modelling approach, Journal of Theoretical Biology 522 (2021) 110685.
  * [27] T. Lorenzi, A. Marciniak-Czochra, T. Stiehl, A structured population model of clonal selection in acute leukemias with multiple maturation stages, Journal of mathematical biology 79 (5) (2019) 1587–1621.
  * [28] G. Clapp, D. Levy, A review of mathematical models for leukemia and lymphoma, Drug Discovery Today: Disease Models 16 (2015) 1–6.
  * [29] A. Ducrot, V. Volpert, On a model of leukemia development with a spatial cell distribution, Mathematical Modelling of Natural Phenomena 2 (3) (2007) 101–120.
  * [30] A. Möricke, M. Zimmermann, A. Reiter, H. Gadner, E. Odenwald, J. Harbott, W.-D. Ludwig, H. Riehm, M. Schrappe, Prognostic impact of age in children and adolescents with acute lymphoblastic leukemia: data from the trials all-bfm 86, 90, and 95, Klinische Pädiatrie 217 (06) (2005) 310–320.
  * [31] Q. Dai, G. Zhang, H. Yang, Y. Wang, L. Ye, L. Peng, R. Shi, S. Guo, J. He, Y. Jiang, Clinical features and outcome of pediatric acute lymphoblastic leukemia with low peripheral blood blast cell count at diagnosis, Medicine 100 (4) (2021).
  * [32] J. A. Bull, H. M. Byrne, The hallmarks of mathematical oncology, Proceedings of the IEEE (2022).
  * [33] C.-H. Pui, W. E. Evans, Treatment of acute lymphoblastic leukemia, New England Journal of Medicine 354 (2) (2006) 166–178.
  * [34] M. Ronghe, G. Burke, S. Lowis, E. Estlin, Remission induction therapy for childhood acute lymphoblastic leukaemia: clinical and cellular pharmacology of vincristine, corticosteroids, l-asparaginase and anthracyclines, Cancer treatment reviews 27 (6) (2001) 327–337.
  * [35] D. Jayachandran, A. E. Rundell, R. E. Hannemann, T. A. Vik, D. Ramkrishna, Optimal chemotherapy for leukemia: a model-based strategy for individualized treatment, PloS one 9 (10) (2014) e109623.
  * [36] C. L. Mouser, E. S. Antoniou, J. Tadros, E. K. Vassiliou, A model of hematopoietic stem cell proliferation under the influence of a chemotherapeutic agent in combination with a hematopoietic inducing agent, Theoretical Biology and Medical Modelling 11 (1) (2014) 1–14.
  * [37] O. Nave, A new protocol applied to cancer treatment-mathematical model-singular perturbed vector field algorithm (2022).
  * [38] S. Rubinow, J. Lebowitz, A mathematical model of the chemotherapeutic treatment of acute myeloblastic leukemia, Biophysical journal 16 (11) (1976) 1257–1271.
  * [39] E. Pefani, N. Panoskaltsis, A. Mantalaris, M. C. Georgiadis, E. N. Pistikopoulos, Chemotherapy drug scheduling for the induction treatment of patients with acute myeloid leukemia, IEEE Transactions on Biomedical Engineering 61 (7) (2014) 2049–2056.
  * [40] Á. Martínez-Rubio, S. Chulián, C. Blázquez Goñi, M. Ramírez Orellana, A. Pérez Martínez, A. Navarro-Zapata, C. Ferreras, V. M. Pérez-García, M. Rosa, A mathematical description of the bone marrow dynamics during CAR T-cell therapy in B-cell childhood acute lymphoblastic leukemia, International Journal of Molecular Sciences 22 (12) (2021) 6371.
  * [41] V. M. Pérez-García, O. León-Triana, M. Rosa, A. Pérez-Martínez, CAR T cells for T-cell leukemias: Insights from mathematical models, Communications in Nonlinear Science and Numerical Simulation 96 (2021) 105684.
  * [42] O. León-Triana, S. Sabir, G. F. Calvo, J. Belmonte-Beitia, S. Chulián, Á. Martínez-Rubio, M. Rosa, A. Pérez-Martínez, M. Ramirez-Orellana, V. M. Pérez-García, CAR T cell therapy in B-cell acute lymphoblastic leukaemia: Insights from mathematical models, Communications in Nonlinear Science and Numerical Simulation 94 (2021) 105570\.
  * [43] G. J. Kimmel, F. L. Locke, P. M. Altrock, The roles of T cell competition and stochastic extinction events in chimeric antigen receptor T cell therapy, Proceedings of the Royal Society B 288 (1947) (2021) 20210229.
  * [44] D. Hanahan, R. A. Weinberg, Hallmarks of cancer: the next generation, Cell 144 (5) (2011) 646–674.
  * [45] M. Karon, E. J. Freireich, E. Frei III, R. Taylor, I. J. Wolman, I. Djerassi, S. L. Lee, A. Sawitsky, J. Hananian, O. Selawry, et al., The role of vincristine in the treatment of childhood acute leukemia, Clinical Pharmacology & Therapeutics 7 (3) (1966) 332–339.
  * [46] D. Campana, Minimal residual disease in acute lymphoblastic leukemia, Hematology 2010, the American Society of Hematology Education Program Book 2010 (1) (2010) 7–12.
  * [47] D. Czock, F. Keller, F. M. Rasche, U. Häussler, Pharmacokinetics and pharmacodynamics of systemically administered glucocorticoids, Clinical pharmacokinetics 44 (1) (2005) 61–98.
  * [48] M. Qweider, J. M. Gilsbach, V. Rohde, Inadvertent intrathecal vincristine administration: a neurosurgical emergency: case report, Journal of Neurosurgery: Spine 6 (3) (2007) 280–283.
  * [49] F. M. Balis, J. S. Holcenberg, W. A. Bleyer, Clinical pharmacokinetics of commonly used anticancer drugs, Clinical pharmacokinetics 8 (3) (1983) 202–232.
  * [50] B. Asselin, C. Rizzari, Asparaginase pharmacokinetics and implications of therapeutic drug monitoring, Leukemia & lymphoma 56 (8) (2015) 2273–2280.
  * [51] G. M. Armstrong, C. P. Midgley, Applications: The exponential-decay law applied to medical dosages., Mathematics Teacher 80 (2) (1987) 110–13.
  * [52] G. B. Haycock, G. J. Schwartz, D. H. Wisotsky, Geometric method for measuring body surface area: a height-weight formula validated in infants, children, and adults, The Journal of pediatrics 93 (1) (1978) 62–66.
  * [53] T. Szczepański, Why and how to quantify minimal residual disease in acute lymphoblastic leukemia?, Leukemia 21 (4) (2007) 622–626.
  * [54] H. Amin, Y. Yang, Y. Shen, E. Estey, F. Giles, S. Pierce, H. Kantarjian, S. O’Brien, I. Jilani, M. Albitar, Having a higher blast percentage in circulation than bone marrow: clinical implications in myelodysplastic syndrome and acute lymphoid and myeloid leukemias, Leukemia 19 (9) (2005) 1567–1572.
  * [55] E. Bianconi, A. Piovesan, F. Facchin, A. Beraudi, R. Casadei, F. Frabetti, L. Vitale, M. C. Pelleri, S. Tassani, F. Piva, et al., An estimation of the number of cells in the human body, Annals of human biology 40 (6) (2013) 463–471.
  * [56] J. Choi, S. J. Lee, Y. A. Lee, H. G. Maeng, J. K. Lee, Y. W. Kang, Reference values for peripheral blood lymphocyte subsets in a healthy korean population, Immune network 14 (6) (2014) 289–295.
  * [57] A. Raes, S. Van Aken, M. Craen, R. Donckerwolcke, J. V. Walle, A reference frame for blood volume in children and adolescents, BMC pediatrics 6 (1) (2006) 1–8.
  * [58] O. Linderkamp, H. Versmold, K. Riegel, K. Betke, Estimation and prediction of blood volume in infants and children, European journal of pediatrics 125 (4) (1977) 227–234.
  * [59] R. Acharya, M. K. Sundareshan, Development of optimal drug administration strategies for cancer-chemotheraphy in the framework of systems theory, International journal of bio-medical computing 15 (2) (1984) 139–150.
  * [60] K. Kay, I. M. Hastings, Improving pharmacokinetic-pharmacodynamic modeling to investigate anti-infective chemotherapy with application to the current generation of antimalarial drugs, PLoS computational biology 9 (7) (2013) e1003151.