# Influence of the age structure on the stability in a tumor-immune model for chronic myeloid leukemia

Kyriaki Dariva, Thomas Lepoutre

► To cite this version:

Kyriaki Dariva, Thomas Lepoutre. Influence of the age structure on the stability in a tumor-immune model for chronic myeloid leukemia. Mathematical Modelling of Natural Phenomena, 2023, 54 (2), pp.619-648. 10.1051/mmnp/2023034. hal-03128770

HAL Id: hal-03128770

https://inria.hal.science/hal-03128770

Submitted on 2 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Influence of the age structure on the stability in a
tumor- immune model for chronic myeloid leukemia

Kyriaki Dariva¹ and Thomas Lepoutre²

¹Univ Lyon, Université Claude Bernard Lyon 1, CNRS UMR5208,
Inria, Institut Camille Jordan, F-69603 Villeurbanne, France

²Univ Lyon, Inria, Université Claude Bernard Lyon 1, CNRS
UMR5208, Institut Camille Jordan, F-69603 Villeurbanne, France

February 2, 2021

**Abstract**

We propose and analyze a system of partial differential equations (PDEs) for chronic myeloid leukemia (CML), generalizing the ordinary differential equations' system (ODE) proposed in [3, 2]. This model describes the proliferation and differentiation of leukemic cells in the bone marrow and the interactions of leukemic and immune cells. We consider that the differentiation of cells can be described by a continuous variable which we use to structure our system. The model is based on a non-monotonic immune response. At low levels, immune response increases with the tumor load whereas for high levels, tumor is suppressing the effect of immune system (immunosuppression). We analyze the stability of the steady states of the model and compare it to the case of [2] where maturity was described as a discrete state. In particular, a steady state describing remission induced by a control due to the immune system is shown to be unstable in certain situations for the PDE model, whereas in [2] it was systematically stable.

# 1 Introduction

The Chronic Myeloid Leukaemia (CML) is a blood and bone marrow cancer. It is characterized by the presence of the fusion gene BCR-ABL1 which is the result of the translocation between the chromosomes 9 and 22. The gene ABL1 of chromosome 9 is juxtaposed onto the gene BCR of chromosome 22
to give the mutated gene BCR-ABL1 which is contained in the chromosome
of Philadelphie (Ph). The presence of chromosome Ph entails the synthesis
of a tyrosine kinase protein which makes the stem cells divide uncontrollably
[5].

Leukemic cells proliferation process begins in the bone marrow by the ac-
tivation of quiescent cancer stem cells. Once activated, these non-differentiated
cells divide, and then differentiate into blood cells. The small number and
the low basal activity of stem cells make any direct characterization of these
cells difficult. Once committed for differentiation, leukemic cells will divide
over twenty times before being released in the blood. They are then generi-
cally called mature cells (actually, in certain phases, immature cells can be
released in the blood). Circulating leukemic cells are not dividing anymore,
and are cleared from the blood within a few days, so the BCR-ABL/ABL
rate represents evolution of the disease in bone marrow, with a delay of few
weeks.

The Tyrosine kinase inhibitor (TKI) Imatinib is the first treatment to
have specifically targeted the BCR-ABL gene in leukemic cells. In the case
of treatment with standard doses, Imatinib limits the leukemic cells' prolif-
eration and promotes apoptosis (cell death) of leukemic cells with acceptable
effects. After the important trial StopImatinib [9] and others that followed
(see [12] for recent review), it has been established that a (small) subset of
patients who achieve deep response can safely stop the treatment. Inter-
estingly, there is increasing evidence that these patients might still harbor
leukemic cells, suggesting control of the disease rather than its eradication.

The immune system seems to play an important role in this control [11, 8, 13], but the effect of the immune system on the leukemic cells is yet to be thoroughly understood. A tumor-immune interaction model focusing on the mechanism of immunosuppression was introduced in [3]. This system was simplified and analyzed in [2]. With this approach, treatment free remission (TFR) is generically interpreted as a remission steady state (terminology introduced in [6]) where the disease is controlled by the immune system without complete eradication of stem cells. This approach has been also considered in [6] and [7]). We also notice that this approach appears very consistent with recent approaches of slower treatment cessation [4].

In this paper we propose a model of PDEs that describes the interac-
tions of immune and leukemic cells. Although we have a more theoretical
approach in this work, one should consider that parameters describing pro-
liferation and death rates of leukemic cells are treatment-dependent (effect
of treatment on other parameters is possible but that would correspond to
less understood mechanisms). Our model is a generalization of a model of
ODEs proposed in [2] where two levels of maturation of leukemic cells are considered: stem cells and fully differentiated (mature). The ODE system is the following:

$$ egin{cases} \frac{dc}{dt} = rc(1 - \frac{c}{K}) - \mu_c cz \\ \frac{dy_2}{dt} = a_1 c - d_2 y_2 - \mu y_2 z \\ \frac{dz}{dt} = s - S(y_2)z \end{cases} \qquad (1) $$

where $y_2$ is the population of mature leukemic cells and $c$ the leukemic stem cells. Our PDE model is the continuous analogue of (1):

$$ egin{cases} \frac{dc}{dt} = rc(1 - \frac{c}{K}) - \mu_c cz \\ \frac{\partial u}{\partial t}(x, t) + g \frac{\partial u}{\partial x}(x, t) = h(x)u(x, t) - \mu_u u(x, t)z \\ gu(0, t) = ac(t) \\ \frac{dz}{dt} = s - S(I)z \end{cases} \qquad (2) $$

with initial conditions $(c^0, u^0, z^0) \in \mathbb{R}_+ \times L_+^1(\mathbb{R}_+) \times \mathbb{R}_+$ where the variables of the model are:

* $c(t)$: the concentration of leukemic stem cells at time $t$.

* $u(x, t)$: the density of leukemic differentiated cells of maturity $x$ at time $t$ and $I = \int_0^\infty u(t, x)dx$.

* $z(t)$: the concentration of immune cells at time $t$.

For the leukemic stem cells we consider logistic growth of rate $r$ and carrying capacity $K$, because this population has self renewal capacities. Let $\mu_c$ be a parameter representing the efficacy of the immune system over the leukemic stem cells.

For the leukemic differentiated cells we suppose that the velocity of differentiation is constant and equal to $g$. The function $h(x) = p(x) - d(x)$ is the balance between the proliferation and natural death rates respectively, two positive real valued functions of the maturity. The parameter $\mu_u$ is the efficacy of the immune cells over the leukemic differentiated cells.

For the immune cells, we consider their natural source. The presence of leukemic cells in the blood on the one hand stimulates the immune system but on the other hand it inhibits the product of immune cells. So, we consider a function $S$ depending on the total leukemic population which describes the stimulation of the immune system. We make the biological assumption that $h$ is negative after a certain level of maturity. This translates into a death rate which exceeds the proliferation rate when maturity is big
enough simply because in practice there is no cells with infinite maturity. When h is constant, this system corresponds to the ODE of [2], but when h really depends on x the PDE model may be richer. In particular we are interested in the effect that a more complex maturity structure can have on the stability of the equilibrium points, especially on the remission steady state. For this purpose we will see the effect of a and h on the stability of the remission steady state, keeping the values of $\bar{c}, \bar{z}, \bar{I}$ constant (so that instability is driven by the shape of h and nothing else).

The stability of the remission steady state is fully characterized by the roots of the characteristic equation,

$$P(\lambda) = Q(\lambda) \int_{0}^{\infty} e^{-\lambda x} \bar{p}(x) dx.$$

where $\bar{p}$ is a probability measure on $\mathbb{R}_+$ corresponding to the normalised distribution of maturity of differentiated cells. The analysis of some extreme cases, where for example $\bar{p}$ takes the form of a Dirac mass, allows us to prove that the steady state can lose stability due to the shape of $\bar{p}(x)$. Moreover it gives us the order of magnitude of the mean age of the distribution that can destabilize the steady state. However, contrary to the linear case (where P is of degree 1 and Q is a constant, we can not establish the optimality of the Dirac mass for this property. Indeed in [1] the authors proved that when P has degree 1 and Q is some constant, then the Dirac is the less stable distribution, meaning that if the Dirac is stable, then all distributions with the same mean are stable too. This leads to a minimal value of the mean delay that can destabilize the steady state. In our case, this property is not satisfied which is established through explicit counterexamples.

The paper is organized as follows. First, we establish the well-posedness of the system under reasonable hypothesis. Second, we characterize the equilibrium points of the system and we study their stability. The results concern linear stability which entails stability of the non linear system (as shown in the Appendix A). Concerning the disease free equilibrium and the high equilibria the structure is shown to be qualitatively similar to [2]. On the other hand, the situation can be very different compared to [2] for the remission equilibrium (if it exists). Its stability can be affected by the distribution of the leukemic differentiated cells. We analyse conditions allowing this instability and finally we investigate the consequence of changing the distribution of maturity on the potential destabilization of steady states obtained for system (1) using realistic parameters.
# 2 Well-posedness and definitions

## 2.1 Immune window

In what follows we suppose that the function S satisfies the following conditions: $S \in C^2(\mathbb{R}^+)$ and that a Y exists such that $S' < 0$ in $(0, Y)$ and $S' > 0$ in $(Y, +\infty)$. As already mentioned, the function S represents the stimulation of the immune system by the leukemic differentiated cells. S can take negative values. We interpret a negative value of this function as the capacity of the immune system to counterbalance the inhibition of the production of immune cells because of leukemic cells.

**Definition 1.** We call immune window the set of I where the function $S(I)$ is negative. If it exists, it has the form of an interval $(y_{\min}, y_{\max})$.

The immune window is characterized as the region where the tumor load is such that the stimulation exceeds the natural loss, leading to expansion of the immune population. The tumor load is neither too high (immunossupression inhibits the growth) nor too low (in which case the immune population may not be stimulated at all).

**Remark** In the case of an immune window, no relevant (i.e. nonnegative) steady state can satisfy $I \in (y_{\min}, y_{\max})$. We are looking for equilibrium points in the zone of I where $S(I) \ge 0$. Consequently, we are looking for equilibrium points either in $\mathbb{R}^+$ or in an interval of the form $(0, y_{\min}) \cup (y_{\max}, +\infty)$.

The existence of a solution of the system is given in Appendix A.

## 2.2 Steady state of the system

In what follows, we will distinguish two sorts of steady states:

* The disease free equilibrium (DFE) $(0, 0, \frac{s}{S(0)})$ which always exists,

* The equilibrium points with positive leukemic load which we call positive equilibria. They are characterized in the following way:

Those that belong to the nondecreasing area, called high equilibria, and those in the nonincreasing area called remission equilibria.

**Proposition 1.** The system has always a disease free equilibrium $(0, 0, \frac{s}{S(0)})$. The positive equilibrium points, if they exist, are fully described by the value
of $\bar{I} = \int_{0}^{\infty} \bar{u}(x)dx$ in the following way :

$$ egin{cases} \bar{c} = \frac{K}{r} \left( r - \mu_c \frac{s}{S(I)} \right) \\ \bar{u}(x) = \frac{a}{g} \frac{K}{r} \left( r - \mu_c \frac{s}{S(I)} \right) e^{\frac{H(x)}{g}} - \mu_u \frac{s}{S(I)} \frac{x}{g} \\ \bar{z} = \frac{s}{S(I)} \end{cases} \quad (3) $$

As a consequence, a positive value of I leads by (3) to a non negative steady state if:

* $S(\bar{I}) > \frac{\mu_c s}{r}$,

* $\bar{I}$ is a solution of the equation $F(I) = I$ where F is defined by

$$ F(I) = \int_{0}^{\infty} \frac{a}{g} \frac{K}{r} \left( r - \mu_c \frac{s}{S(I)} \right) e^{\frac{H(x)}{g}} - \mu_u \frac{s}{S(I)} \frac{x}{g} dx $$

*Proof.* We are looking for stationary points of the system. These are the solutions of the following system:

$$ egin{cases} 0 = r\bar{c}\left(1 - \frac{\bar{c}}{K}\right) - \mu_c\bar{c}\bar{z} \\ g\frac{d\bar{u}}{dx}(x) = h(x)\bar{u}(x) - \mu_u\bar{u}(x)z \\ g\bar{u}(0) = a\bar{c} \\ 0 = s - S(\bar{I})\bar{z} \end{cases} $$

We solve the system and we obtain:

$$ egin{cases} 0 = \bar{c}\left(\left(r - \frac{r\bar{c}}{K}\right) - \mu_c\bar{z}\right) \\ \bar{u}(x) = \bar{u}(0)e^{\frac{1}{g}\int_0^x h(s)ds - \frac{1}{g}\mu_u\bar{z}x} \\ g\bar{u}(0) = a\bar{c} \\ \bar{z} = \frac{s}{S(I)} \end{cases} $$

We deduce that: Either $\bar{c} = 0$ and consequently $\bar{u}(0) = 0$. Hence $\bar{u}(x) = 0$ and $\bar{z} = \frac{s}{S(0)}$. In this case we have the disease free equilibrium point

$$ (\bar{c}, \bar{u}, \bar{z}) = (0, 0, \frac{s}{S(0)}) $$

If we have $\bar{c} \neq 0$ then immediately $\bar{c} = K - \frac{\mu_c \bar{z}}{r}$, $\bar{z} = \frac{s}{S(I)}$ which leads to the condition $S(I) > 0$ to make sure that $\bar{z} > 0$ and to the condition $S(\bar{I}) > \frac{\mu_c s}{r}$
to have $\bar{c} > 0$. The formulation for $\bar{u}$ is obtained after integration along characteristics and the fixed point formulation is then just a compatibility condition. $\square$

**Proposition 2.** • If $rS(0) - \mu_c s > 0$ there exists at least one positive equilibrium.

Moreover, if the function $S(I) - \mu_c s/r$ has two zeros, $I_{min}$, $I_{max}$, then there exists a unique positive equilibrium parametrized by a value of $I$ in the zone $(0, I_{min})$, that we call remission equilibrium. All other equilibria, if they exist, they are parametrized by value of $I$ in $(I_{max}, \infty)$ and we call them high equilibrium points.

• If $rS(0) - \mu_c s \le 0$ then no remission equilibrium can exist (high equilibria may exist in $(y_{max}, +\infty)$)

*Proof.* We reformulate the fixed point:
We introduce the function

$$f(I) = \frac{\mu_u}{\mu_c g} \left( r - \mu_c \frac{s}{S(I)} \right)$$

and the function $\psi$

$$\psi(x) = a \frac{\mu_c K}{r \mu_u} e^{\frac{1}{g} \int_0^x h(s) ds - \frac{r \mu_u}{\mu_c g} x}$$

so that $F$ is written in the condensed form:

$$F(I) = f(I) \int_{0}^{\infty} \psi(x) e^{f(I)x} dx$$

(a) We observe that $f$ starts by being positive: $f(0) = \frac{\mu_u}{\mu_c g} (r - \frac{\mu_c s}{S(0)}) = \frac{\mu_u S(0) r - \mu_c s}{\mu_c g S(0)} > 0$.

We compute the derivative of $f$:

$$f'(I) = \frac{s \mu_u}{g} \frac{S'(I)}{S(I)^2}$$

and the derivative of $F$:

$$F'(I) = f'(I) \int_{0}^{\infty} \psi(x) e^{f(I)x} dx + f(I) f'(I) \int_{0}^{\infty} x \psi(x) e^{f(I)x} dx$$

The function $f$ has the same monotony as $S$. We have two situations: Either $f$ is nonnegative on $\mathbb{R}^+$ and the existence of an equilibrium is simply
due to the fact that the continuous function $F(I) - I$ (defined then on $\mathbb{R}_+$) changes sign between 0 and $+\infty$ ($F(\infty)$ is finite). If $f$ vanishes, we note $I_{min}$ its first zero. By construction, $f$ is decreasing on $(0, Y)$. If it does not vanish before $Y$, it does not vanish at all. If it does, the first zero $I_{min}$ is located in $[O, Y]$. In that case $F(I) - I$ changes sign between 0 and $I_{min}$ and it is also decreasing, making the fixed point unique in $[0, I_{min}]$. Due to its shape, $f$ vanishes at most twice. We call $I_{max}$ the second zero (if it exists). By a symmetric argument, we have $Y \le I_{max}$ and $S'(I) \ge 0$ on $[I_{max}, +\infty[$. Since there cannot be any fixed points of $F$ in $[I_{min}, I_{max}]$ all the remaining equilibria can only correspond to values of $I$ located in $]I_{max}, +\infty[ \subset ]Y, +\infty[$ and are therefore high equilibria. If $I_{max} > F(0)$ there is no intersection with the bisector after $y_{max}$. Otherwise, the number of intersection points depends on the variations of $F$.

(b) In that case $f$ starts by being non positive: $f(0) = \frac{\mu_u S(0)r - \mu_c s}{g\mu_c S(0)} \le 0$. As $f$ is decreasing until $Y$ and increasing afterwards, we see that: If $f$ has no zeros then $f$ is always negative and there is no intersection with the bisector. If $f$ has a zero in $(y_{max}, +\infty)$ then the intersections with the bisector, if any, correspond to equilibrium points with high tumor concentration. $\square$

**Remark:** If there is no immune window (S(I) positive on $\mathbb{R}^+$) there is always an equilibrium but we do not call it remission equilibrium point anymore.

# 3 Stability analysis

## 3.1 From linear to nonlinear stability

We are now ready to study the stability of the identified equilibrium points.

**Theorem 1.** If a stationary point is linearly exponentially stable, then it is also non linearly stable.

For the proof see Appendix A.

The linear system around an equilibrium is written as:

$$ 
\left\{ 
\begin{aligned}
\frac{dc}{dt} &= (r - \frac{2r\bar{c}}{K} - \mu_c \bar{z})c - \mu_c \bar{c} z \\
\frac{\partial u}{\partial t}(x, t) + g \frac{\partial u}{\partial x}(x, t) &= h(x)u(x, t) - \mu_u \bar{u} z - \mu_u u \bar{z} \\
u(x, t) &= \frac{ac(t)}{g} \\
\frac{dz}{dt} &= -S(\bar{I})z - \bar{z} S'(\bar{I})I.
\end{aligned}
\right.
\quad (4) $$

We have proved that the nonlinear stability stems from the behavior of (4).
We consider the semi group $T(t)$ defined on $X = \mathbb{R}^{+} \times L^{1}(\mathbb{R}^{+}) \times \mathbb{R}^{+}$ by:

$$T(t) : X \to X$$

$$(c_0, u_0(x), z_0) \mapsto (c(t), u(x, t), z(t))$$

which to an initial condition associates a solution of (L).

Let A be its infinitesimal generator.

Then the linear system can be written like:

$$\begin{pmatrix} \frac{dc}{dt} \\ \frac{du}{dt} \\ \frac{dz}{dt} \end{pmatrix} = A \begin{pmatrix} c \\ u \\ z \end{pmatrix}$$

**Proposition 3.** Consider a steady state $(\bar{c}, \bar{u}, \bar{z})$ and a complex number $\lambda$ satisfying $Re(\lambda) > -\mu_u \bar{z}$. We introduce the notation

$$E_{\lambda}(x) = \frac{a}{g} e^{\frac{1}{g} \left( \int_{0}^{x} h(s) ds - \mu_{u} \bar{z} x \right)} e^{-\frac{\lambda x}{g}}$$

and the matrix $M(\lambda)$

$$M(\lambda) = \begin{bmatrix} \lambda + 2 \frac{r\bar{c}}{K} - r + \mu_c \bar{z} & \mu_c \bar{c} \\ S'(\bar{I}) \bar{z} \int_0^\infty E_\lambda(x) dx & \lambda + S(\bar{I}) - \bar{z} \mu_u S'(\bar{I}) \int_0^\infty \bar{u}(x) \int_0^x \frac{e^{-\frac{\lambda(x-y)}{g}}}{g} dy dx \end{bmatrix}$$

Then M is well defined and $\lambda$ satisfies necessarily one of these two properties

* $\det(M) \neq 0$ and $(A - \lambda)$ is invertible ($\lambda$ is in the resolvent set),

* $\lambda$ is an eigenvalue of A and $\det M = 0$.

*Proof.* We try to solve the resolvent operator for a given $\lambda$, that is we look for a nontrivial solution $(C, U, Z)$ of the problem with a source

$$\begin{cases} (\lambda + 2 \frac{r\bar{c}}{K} - r + \mu_c)C + \mu_c \bar{c} Z = c_1, \\ \lambda U(x) + g \frac{dU(x)}{dx} - h(x)U(x) + \mu_u \bar{z} U + \mu_u \bar{u} Z = u_1(x), \\ gU(0) = aC, \\ \lambda Z + S(\bar{I})Z + \bar{z} S'(\bar{I}) \int_0^\infty U(x) dx = z_1 \end{cases}$$

With our notations, then,

$$U(x) = C E_{\lambda}(x) - \frac{\mu_u Z}{g} \int_{0}^{x} \frac{E_{\lambda}(x)}{E_{\lambda}(y)} \bar{u}(y) dy + \int_{0}^{x} \frac{E_{\lambda}(x)}{E_{\lambda}(y)} u_{1}(y) dy$$
As $\Re(\lambda) > -\mu_u \bar{z}$, $u_1 \in L^1$ we have no integrability issues. Indeed, $\int_y^x h$ is bounded from above if $y < x$, so $E_\lambda(x)/E_\lambda(y)$ is controlled by $Ke^{-(Re(\lambda)+\mu_u \bar{z})(x-y)} = Ke^{-\epsilon(x-y)}$ with $\epsilon > 0$. In particular, it is straightforward to see that:

$$ \|E_\lambda\|_1 \le \frac{K}{\epsilon}, \quad \left\| \int_0^x \frac{E_\lambda(x)}{E_\lambda(y)} \bar{u}(y) dy \right\|_1 \le \frac{K}{\epsilon} \|\bar{u}\|_1, \quad \left\| \int_0^x \frac{E_\lambda(x)}{E_\lambda(y)} u_1(y) dy \right\|_1 \le \frac{K}{\epsilon} \|u_1\|_1. $$

So the problem has a solution iff we can solve the linear problem in $C, Z$

$$ \left\{ \begin{array}{l} (\lambda + 2 \frac{r \bar{c}}{K} - r + \mu_c)C + \mu_c \bar{c} Z = c_1 \\ \frac{\bar{z} S'(\bar{I})}{g} \int_0^\infty E_\lambda(y) dy C + (\lambda + S(\bar{I}) - \bar{z} S'(\bar{I}) \frac{\mu_u}{g}) \int_0^\infty \int_0^x \frac{E_\lambda(x)}{E_\lambda(y)} \bar{u}(y) dy dx Z \\ = z_1 - \bar{z} S'(\bar{I}) \frac{1}{g} \int_0^\infty \int_0^x \frac{E_\lambda(x)}{E_\lambda(y)} u_1(y) dy dx \end{array} \right. $$

The problem can be written in terms of matrices if we notice that:

$$ \frac{E_\lambda(x)}{E_\lambda(y)} \bar{u}(y) = \bar{u}(x) e^{-\lambda(x-y)}, $$

Hence the aforementioned matrix $M(\lambda)$ appears. Then we have the following alternative:

* if the matrix is not invertible, there exists a nontrivial solution associated to a source 0 which leads to the construction of an eigenvector,

* in the opposite case, the problem is uniquely solvable which makes $(A - \lambda)$ invertible.

This ends the proof of the proposition. □

**Proposition 4.** The growth bound $\omega_0$ of the linear semi group satisfies $\omega_0 \le \max(-\mu_u \bar{z}, s(A))$ where we have classically denoted $s(A) = \sup_{\lambda \in \sigma(A)} Re(\lambda)$. In particular if $s(A) < 0$ the steady state is linearly (and also non-linearly) stable.

*Proof.* Firstly, observe that the difficulty derives from some lack of compactness (if maturity lied in a bounded interval, we would be able to derive the result immediately from an eventual compactness). For the proof we will need some properties and definitions from spectral analysis given in Appendix B. First notice that all elements of $\sigma(T(t))$ have exponential form. For the point and the residual spectrum it is a result of spectral theorem. For the approximate spectrum it is immediate from its limit property. Indeed if
$\mu$ is an approximate eigenvalue of $T(t)$ and $x_0^n$ an approximate eigenvector then taking the limit:

$$ \lim_{n \to +\infty} \|T(t)x_0^n - \mu x_0^n\| = 0 $$

we finally get after some easy calculation that the limit of the approximate eigenvector $x_0$ is a proper vector hence $\mu$ has exponential form lets say $e^{\lambda t}$. Then the eigenvector and the eigenvalue define a solution with exponential profile of the linear problem, let's say: $c_n^0 e^{\lambda t}, u_n^0(x) e^{\lambda t}, z_n^0 e^{\lambda t}$. Back to our proof, it suffices to prove that the approximate spectrum of $T_L(t)$ does not contain any element $e^{\lambda t}$ with $Re(\lambda) > -\mu_u \bar{z}$. Therefore, we consider a sequence $x_n^0$ such that $\|T_L(t)x_n^0 - e^{\lambda t}x_0^n\| \to 0$. Let $c_n^0, u_n^0, z_n^0$ be the components of $x^0$ and define $c_n(s), u_n(s), z_n(s)$ as the solution of the linear problem with initial data $c_n^0, u_n^0, z_n^0$. It is straightforward to derive a bound on

$$ \sup_{[0,t]} |c_n(s)| + \int_0^\infty |u_n(s,x)|dx + |z_n(s)| \le M \max(1, e^{\omega_0 t}). $$

From the Arzela-Ascoli theorem in the components associated to $c, z$ we have the convergence (up to subsequence):

$$ \sup_{[0,t]} |c_n(s) - c^\infty(s)| + |z_n(s) - z^\infty(s)| \to 0. $$

Then, keeping the previous notations, we can solve the equations on $u_n$ leading to:

$$ u_n(t,x) = \begin{cases} v_n^0(x - gt) \frac{E_0(x)}{E_0(x-gt)} - \mu_u \bar{u}(x) \int_0^t z_n(s) ds & \text{if } x - gt > 0, \\ \frac{c_n(t-x/g)}{g} E_0(x) - \mu_u \bar{u}(x) \int_0^{x/g} z_n(t - x/g + s), & \text{if } x - gt < 0 \end{cases} $$

Hence we infer the (pointwise) convergence of $u_n^0$ on $[0, gt]$ (simply using the second line) towards :

$$ \begin{aligned} u_n^0(x) &\to e^{-\lambda t} \left( \frac{c^\infty\left(t - \frac{x}{g}\right)}{g} E_0(x) - \mu_u \bar{u}(x) \int_{t-x/g}^t z^\infty(s) ds \right) \\ &= \frac{e^{-\lambda(t-x/g)} c^\infty\left(t - \frac{x}{g}\right)}{g} E_\lambda(x) - e^{-\lambda t} \mu_u \bar{u}(x) \int_{t-x/g}^t z^\infty(s) ds \end{aligned} $$

Noting

$$ \frac{E_0(x)}{E_0(x - gt)} e^{-\lambda t} = \frac{E_\lambda(x)}{E_\lambda(x - gt)}, \quad \frac{E_0(x + gt)}{E_0(x)} \bar{u}(x) = \bar{u}(x + gt), $$
Since we have by construction, $c^\infty(t) = c^\infty(0)e^{\lambda t}$, we can naturally extend $s \mapsto c^\infty(s)e^{-\lambda s}$ into a (continuous) t periodic function $c_p^\infty(\cdot)$ on $\mathbb{R}$. Similarly, we define $z_p^\infty(\cdot)$. This writing is then quite convenient since it allows a condensed formula for the pointwise limit:

$$u^\infty(x) = \frac{c_p^\infty\left(t - \frac{x}{g}\right)}{g} E_\lambda(x) - \mu_u \bar{u}(x) \int_{t-x/g}^{t} e^{\lambda(s-t)} z_p^\infty(s) ds$$

If the above limits are nonzero, we have just built an eigenvector of $T(t)$ with eigenvalue $e^{\lambda t}$ (and we are dealing with the point spectrum instead of the approximate point spectrum). In this case, we have $\lambda \le s(A)$. On the other hand, if the pointwise limit is 0, it means that the sequence presents compactness problem. By construction $u_0^n \to u^\infty$ in $L^1$ on compact intervals. Since the total mass is preserved, it means that we have mass going to $\infty$. In this case, the $c, z$ components go to zero and hence so does $u_0^n$ on any compact interval. So for any fixed $R > 0$ we have:

$$\|T(t)x^n\| = \int_R^\infty |u^n(t, x)|dx + o(1) \le \|u_0^n\| \sup_{R, \infty} \frac{E_0(x + gt)}{E_0(x)} + o(1)$$

So

$$e^{\lambda t} \le \limsup_{x \to +\infty} \frac{E_0(x + gt)}{E_0(x)} = \limsup_{x \to +\infty} e^{\int_x^{x+gt} \frac{h-\mu_u \bar{z}}{g}} \le e^{-\mu_u \bar{z}t}.$$

In this case, the hypothesis $e^{\lambda t} \le e^{-\mu_u \bar{z}t}$ entails that $\lambda \le -\mu_u \bar{z}$. $\square$

<figure><img src="image_2.png" /></figure>

## 3.2 Disease free equilibrium and high steady states

First consequence of the former section concerns the disease free equilibrium.

**Proposition 5.** For the disease free steady state, the growth bound of the linearized semigroup satisfies $\omega_0 \le \max(-\mu_u \bar{z}, r - \mu_c s/S(0), -S(0))$. In particular, the DFE is linearly unstable if $r - \mu_c s/S(0) > 0$ and linearly stable if $r - \mu_c s/S(0) < 0$.

*Proof.* In this case we have $\bar{c} = 0, \bar{u} = 0, \bar{z} = s/S(0)$, so the matrix $M(\lambda)$ is upper triangular and $\det(M(\lambda)) = (\lambda + r + \mu_c s/S(0))(\lambda + S(0))$. If $\max(-\mu_u \bar{z}, r - \mu_c s/S(0), -S(0)) > 0$ then $s(A) \ge \omega_0$, hence $s(A) = \omega_0$. $\square$

We turn now to the high steady states for which $S$ is increasing.

**Proposition 6.** Stability of high equilibrium points is given by real eigenvalues.
The main reason is that in this case the induced linear semigroup applied to $(c, u, -z)$ (one has to change the sign for the variable $z$) is a positive semigroup. In any case, we can establish it directly from the equation on $det M(\lambda)$ from proposition 3 to have immediately a link with the alternated stability shown in [2].

For the positive steady states, it is more convenient to use $2^{r\bar{c}/K} + \mu_c\bar{z} - r = r\bar{c}/K$ so the determinant can be written as

$$ \det M(\lambda) = \left(\lambda + \frac{r\bar{c}}{K}\right) \left(\lambda + S(\bar{I}) - \bar{z}\mu_u S'(\bar{I}) \int_0^\infty \bar{u}(x) \int_0^{x/g} e^{-\lambda y} dy dx\right) \\ - \mu_c \bar{c} S'(\bar{I}) \bar{z} \int_0^\infty E_\lambda(x) dx \qquad (5) $$

For $\Re(\lambda) > -min(r\bar{c}/K, \mu_u\bar{z})$ and $S'(\bar{I}) > 0$ we can prove by using the following trivial results:

$$ \Re(E_{\lambda}(x)) \le |E_{\lambda}(x)| \le E_{\Re(\lambda)}(x), \quad \left|\frac{1}{(\lambda + a)}\right| \le \frac{1}{\Re(\lambda) + a}, \Re(\lambda) + a > 0. $$

that

$$ \left| \frac{\det M(\lambda)}{(\lambda + \frac{r\bar{c}}{K})(\lambda + S(\bar{I}))} \right| \ge 1 - \bar{z}\mu_u S'(\bar{I}) \int_0^\infty \bar{u}(x) \int_0^{x/g} \frac{e^{-\Re(\lambda)y}}{\Re(\lambda) + S(\bar{I})} dy dx \\ - \mu_c \bar{c} S'(\bar{I}) \bar{z} \int_0^\infty \frac{E_{\Re(\lambda)}(x)}{(\Re(\lambda) + \frac{r\bar{c}}{K})(\Re(\lambda) + S(\bar{I}))} dx \\ = \frac{\det M(\Re(\lambda))}{(\Re(\lambda) + \frac{r\bar{c}}{K})(\Re(\lambda) + S(\bar{I}))}. $$

The inequality is strict for $\lambda$ complex. In particular $det(M(\lambda)) = 0$ implies then $det(M(\Re(\lambda)) < 0$. Since the function

$$ \lambda \mapsto \frac{\det M(\lambda)}{(\lambda + \frac{r\bar{c}}{K})(\lambda + S(\bar{I}))} $$

(seen as function over the reals) is increasing on $]-\min(\frac{r\bar{c}}{K}, \mu_u\bar{z}, S(\bar{I})), +\infty[$ and becomes 0 for any (real) eigenvalue, the stability is determined by a real eigenvalue. In particular, this function changes sign at most once. If it remains positive on the whole interval, $det(M)$ cannot vanish for any $\lambda$ satisfying $Re(\lambda) \ge -min(\frac{r\bar{c}}{K}, \mu_u\bar{z}, S(\bar{I}))$ and consequently $\omega_0 < 0$ and $det(M) > 0$. Otherwise, $\omega_0$ is precisely the unique zero of the function and
therefore since the function is increasing, we have $det(M(0)) > 0$ if $\omega_0 < 0$ and vice versa.

We finish this section with an equivalent result as in theorem 2, Chapter 2 of [2].

**Proposition 7.** *The stability of high equilibrium points is alternated. The set of I associated to high equilibria has a minimum and the equilibrium associated to this minimum is unstable.*

*Proof.* With the same computations as above, we have that $det(M(0))$ is equal to:

$$ 
\begin{aligned}
\det(M(0)) &= 1 - \bar{z}\mu_u S'(\bar{I}) \int_0^\infty \bar{u}(x) \frac{x}{gS(\bar{I})} dydx \\
&\quad - \mu_c \bar{c} S'(\bar{I}) \bar{z} \int_0^\infty \frac{E_0(x)}{\frac{r\bar{c}}{K} S(\bar{I})} dx \\
&= 1 - F'(\bar{I})
\end{aligned}
$$

In particular, $\omega_0$ has the sign of $F'(\bar{I}) - 1$. Since the steady states are characterized as fixed points of $F$ and in the considered region $F$ is increasing, at the lowest (high) fixed point $F'(\bar{I}) > 1$, then at the next one $F'(\bar{I}) < 1$ etc (except of double roots where stability is not given by $\omega_0$). The last one (which exists because $F$ is bounded and hence we have a maximum number for solutions) is generically stable. The situation is qualitatively the same as in [2] (the number of high steady states might depend on the structure but this point is not as crucial as the potential destabilization of the remission steady state). $\square$

### 3.3 Stability of remission equilibrium

Since our goal is to investigate the effect of the maturity in the structure of our system, we will work with fixed values of $\bar{c}, \bar{I}, \bar{z}$ in order to see how the shape of $h$ can affect the result.

**Proposition 8.** *For a remission steady state, there cannot be a real eigenvalue such that*

$$ \lambda > -\min(\frac{r\bar{c}}{K}, \mu_u \bar{z}, S(\bar{I})) $$

*Proof.* The computation is the same as for high equilibrium points but now $S'(I) < 0$, so we have:

$$ \lambda > -\min(\frac{r\bar{c}}{K}, \mu_u \bar{z}, S(\bar{I})) \rightarrow \frac{\det(M(\lambda))}{(\lambda + \frac{r\bar{c}}{K})(\lambda + S(\bar{I}))} > 0 $$
which leads to the conclusion.

□

As a direct consequence, we can now simplify the problem.

**Proposition 9.** Every eigenvalue associated to the remission steady state, such that

$$Re(\lambda) > -\min(\frac{r\bar{c}}{K}, \mu_u\bar{z}, S(\bar{I}))$$

is a nonzero root of the equation

$$P(\lambda) + Q(\lambda) \int_{0}^{\infty} \frac{\bar{u}(x)}{\bar{I}} e^{-\lambda \frac{x}{g}} dx = 0, \quad (6)$$

where P, Q are the following polynomials:

$$P(\lambda) = \lambda^3 + (R + S)\lambda^2 + (RS + D)\lambda + DR,$$

$$Q(\lambda) = \left( \frac{\mu_c}{\mu_u} - 1 \right) D\lambda - DR$$

where we have introduced notations for the positive quantities

$$R = \frac{r\bar{c}}{K}, \quad S = S(\bar{I}), \quad D = -\mu_u s \frac{S'(\bar{I})\bar{I}}{S(\bar{I})}$$

*Proof.* Since 0 is not a root, we simply write from (5) the equation satisfied by $\lambda \det(M)$ when $\lambda \neq 0$. In particular, in this case

$$\lambda \int_{0}^{x/g} e^{-\lambda y} dy = (1 - e^{-\lambda x/g}).$$

With simple computations we derive:

$$\lambda \det(M(\lambda)) = \lambda(\lambda + R)(\lambda + S) - \bar{z}\mu_u S'(\bar{I})\bar{I}(\lambda + R) + \quad (7)$$

$$\bar{z}\mu_u S'(\bar{I})(\lambda + R)\left(\int_0^\infty u(x)e^{-\lambda x/g}dx\right) - \lambda\mu_u \bar{c}S'(\bar{I})\bar{c}\int_0^\infty \frac{1}{\bar{c}}u(x)e^{-\lambda x/g}.dx \quad (8)$$

Replacing $D = -\mu_u s \frac{S'(\bar{I})}{S(\bar{I})} \bar{I}$ and extending the polynomials, we obtain the result. □
Once we have established this expression of the characteristic equation, we can discuss the influence of the shape of the probability distribution $\bar{u}/\bar{I}$ on the stability. We can have a glimpse on the stability boundary in specific cases. In our equation there are two parameters that can be modified without changing the values of $\bar{c}, \bar{I}, \bar{z}$: that is the shape of the distribution $\bar{p} = \bar{u}/\bar{I}$ and the slope of the derivative $S'(\bar{I})$ which is encoded in the positive parameter D.

**The case of exponential distribution [2]**: The case of an exponential distribution, $\bar{p} = \gamma_0 e^{-\gamma_0 x}$ corresponds to the study done in [2]. In this case, the eigenvalues above $-\mu_u \bar{z}$ are the (nonzero) solutions of the equation:

$$ \lambda(\lambda + R)(\lambda + S) + D(\lambda + R) + D \left( \left( \frac{\mu_c}{\mu_u} - 1 \right) \lambda - R \right) \frac{\gamma_0}{\gamma_0 + \frac{\lambda}{g}} = 0. $$

This can be reduced (since 0 is not an acceptable root) to the third degree polynomial:

$$ (\lambda + \gamma_0 g)(\lambda + R)(\lambda + S) + D(\lambda + R + \frac{g\mu_c}{\mu_u}) = 0. $$

And then, we can check that any solution satisfies:

$$ Re(\lambda) \leq -\min(R, S, \gamma_0 g) < 0. $$

**Stability boundary for a Dirac distribution in $\tau$, when $\mu_c = \mu_u$ and $g = 1$**.

We are now investigating the extreme case of a Dirac distribution, for which explicit computations can be done.

To proceed it is useful to make the following changes on the characteristic equation:

$$ (\lambda + R)(\lambda + S) + D + DRE \left[ \frac{1 - e^{-\lambda X}}{\lambda} \right] = 0, $$

where E stands for the mean with respect to a probability distribution. Divide by $\frac{1}{R^2}$ and set $\bar{\lambda} = \frac{\lambda}{R}, \bar{S} = \frac{S}{R}, \bar{D} = \frac{D}{R^2}$. Hence,

$$ (\bar{\lambda} + 1)(\bar{\lambda} + \bar{S}) + \bar{D} + \bar{D}E \left[ \frac{1 - e^{-\bar{\lambda}Y}}{\bar{\lambda}} \right] = 0 $$

Where $Y = RX$. If $X$ is distributed with probability function $\bar{p}(x)dx$, then $Y$ has probability function $\bar{p}(\frac{y}{R})\frac{dy}{R}$.

In that, if $Y$ is destabilized with 1, $\bar{S}, \bar{D}$ (the characteristic equation has solutions with positive real parts), then $X$ is destabilized with $R, S, D$ and inversely.
To simplify the notations, we will omit using tilda notations. The equation that we will be studying is:

$$(\lambda + 1)(\lambda + S) + D + DE\left[\frac{1 - e^{\lambda Y}}{\lambda}\right] = 0$$

We are looking for purely imaginary roots $i\omega$ of the equations. This leads to $\omega > 0$ solution of:

$$-i\omega^3 - \omega^2(1 + S) + i\omega(S + D) + D = De^{-i\omega\tau} \quad (9)$$

It is very useful to notice that this equation implies in particular (comparing modulus of both sides)

$$(\omega^3 - (S + D)\omega)^2 + (\omega^2(1 + S) - D)^2 = D^2 \quad (10)$$

Which can be simplified by writing $y = \omega^2$

$$y^2 + y(S^2 + 1 - 2D) + S^2 + D^2 - 2D = 0$$

since 0 is not a solution. which is a quadratic polynomial of y. We can also rewrite it as a quadratic polynomial of D:

$$D^2 - 2D(y + 1) + y^2 + y(1 + S^2) + S^2 = 0.$$

In that we have two different approaches. We can either solve D as a function of $\omega^2$ or $\omega^2$ as a function of D.

**Lemma 1 (Solving $D(\omega)$).** Given a frequency $\omega > 0$, there exist (up to multiplicity)

* Either two positive solutions to the equation (10) if

$$(1 - S^2)(\omega^2 + 1) \geq 0 \quad (11)$$

which is a double root if this quantity is zero.

They are given by:

$$D_{\pm} = \omega^2 + 1 \pm \sqrt{(1 - S^2)(\omega^2 + 1)}$$

* or no nonnegative solutions otherwise.
When we start with S and D given, then a necessary condition for solutions to (10) to exist is:

$$S \le 1 \qquad (12)$$

and

$$D \ge 1 - \sqrt{1 - S^2}. \qquad (13)$$

*Proof.* If we consider D as the unknown, we are looking for positive roots of the quadratic polynomial:

$$P_{\omega^2}(D) = D^2 - 2D(\omega^2 + 1) + (\omega^4 + S^2 + \omega^2(1 + S^2)).$$

Noticing that this is decreasing on $\mathbb{R}_-$, and $P_{\omega^2}(0) > 0$, we see immediately that real roots are necessarily positive, so that the discriminant is non negative.

$$\Delta = 4(\omega^2 + 1)^2 - 4(\omega^4 + S^2 + \omega^2(1 + S^2)) = 4(1 - S^2)(\omega^2 + 1)$$

leading to condition (11). Calculation of the roots is then straightforward and we obtain the last condition by minimizing $D_-$.

**Lemma 2.** Assume (12) is satisfied, then a solution to (10) is given by

$$\omega_+^2 = \frac{1}{2} \left( 2D - 1 - S^2 + (1 - S^2) \sqrt{1 + \frac{4D}{1 - S^2}} \right). \qquad (14)$$

If additionally $D \ge (1 + \sqrt{1 - S^2})$, then we have a second solution that makes sense,

$$\omega_-^2 = \frac{1}{2} \left( 2D - 1 - S^2 - (1 - S^2) \sqrt{1 + \frac{4D}{1 - S^2}} \right). \qquad (15)$$

We also have the identities

$$\begin{cases} \forall D \ge (1 + \sqrt{1 - S^2}), \quad D_-(\omega_+(D)) = D, \\ \forall D \ge (1 + \sqrt{1 - S^2}) \quad D_+(\omega_-(D)) = D, \\ \forall \omega > 0, \quad \omega = \omega_+(D_-(\omega)) = \omega_-(D_+(\omega)) \end{cases} \qquad (16)$$

**Corollary 1.** For a general distribution (not a Dirac), if we have a solution of the form $i\omega$ to (6), then necessarily, we are in the conditions (12) and $\omega^2$ satisfies the inequalities

$$\begin{cases} 0 < \omega^2 \le \omega_+(D), & \text{if } D < (1 + \sqrt{1 - S^2}), \\ \omega_-^2(D) \le \omega^2 \le \omega_+(D), & \text{if } D \ge (1 + \sqrt{1 - S^2}) \end{cases} \qquad (17)$$
Proof. We simply solve the quadratic polynomial. If $D < (1 + \sqrt{1 - S^2})$ there is only one root. Otherwise there are two. The quadratic polynomial is negative in the regions of $\omega^2$ defined above. $\square$

**Lemma 3.** *The boundary is given by the graph:*

$$ \tau(D) = \tau(\omega_+(D)) = \frac{1}{\omega_+(D)} \left( \frac{3\pi}{2} + \arctan\left(\frac{D - \omega_+^2(D)}{\omega_+(D)S}\right) - \arctan(\omega_+(D)) \right). $$

*Proof.* From (9) we have:

$$ \begin{cases} -\omega^3 + \omega(S + D) = D\sin(-\omega\tau) \\ -\omega^2(1 + S) + D = D\cos(-\omega\tau) \end{cases} $$

So

$$ \frac{P(i\omega)}{D} = E[e^{iY}] $$

with respect to $\delta_{\tau}$. Where $P(\lambda) = -\lambda^3 + (1+S)\lambda^2 - (S+D)\lambda + D$. Evidently $\frac{P(i\omega)}{D}$ is a complex number of modulus 1.

Hence

$$ \left( \frac{\cos(\omega\tau)}{\sin(\omega\tau)} \right) = \left( \frac{\frac{-\omega^2(1+S)+D}{D}}{\frac{\omega^3-(S+D)\omega}{D}} \right) $$

Multiplying by

$$ \sqrt{\omega^2 + 1} (\sin(\arctan(\omega)) \cos(\arctan(\omega)) = (\omega \ 1), $$

We obtain

$$ D\sqrt{\omega^2 + 1} \sin(\omega\tau + \arctan(\omega)) = -\omega S(\omega^2 + 1) < 0 $$

With similar computations, we have :

$$ D\sqrt{\omega^2 + 1} \cos(\omega\tau + \arctan(\omega)) = (1 + \omega^2)(D - \omega^2). $$

Consequently

$$ \omega\tau + \arctan(\omega) \in [\pi, 2\pi] $$

or equivalently

$$ -\frac{\pi}{2} < \omega\tau + \arctan(\omega) - \frac{3\pi}{2} < \frac{\pi}{2} $$

So

$$ \tan(\omega\tau + \arctan(\omega) - \frac{3\pi}{2}) = \frac{\cos(\omega\tau + \arctan(\omega))}{-\sin(\omega\tau + \arctan(\omega))} = \frac{(1 + \omega^2)(D - \omega^2)}{\omega S(\omega^2 + 1)}. $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Stability boundary given by a single Dirac distribution</figcaption></figure>

Finally

$$ \tau = \frac{1}{\omega} \left( \frac{3\pi}{2} + \arctan\left(\frac{D - \omega^2}{\omega S}\right) - \arctan(\omega) \right). $$

Clearly $\tau$ is a decreasing function of $\omega$. So the boundary is given by $\omega_+$ and not $\omega_-$. □

We will be always looking for the minimal positive $\tau$.

**Stability boundary for two Dirac masses: Suboptimality of the Dirac solutions**

The characteristic equation for $\mu_u = \mu_c$ is the following:

$$ \frac{P(\lambda)}{D} = E(e^{-\lambda X}) \qquad (18) $$

where $E$ denotes the mean value with respect to the distribution $\bar{p}(x)$.

In the case of a first degree polynomial $P(\lambda) = \lambda + a$ it has been established in [1] that the Dirac is the less stable distribution among all distributions with the same mean. This means that if a distribution exists with mean $\bar{X}$ and for which we have instability, then the Dirac $\delta_{\bar{X}}$ is also unstable.
In what follows we will prove that this is not the case for our model. In particular, we will prove that for certain parameters we can find an unstable distribution, which is not a Dirac, with mean $\bar{X} < \tau_+$.

The question that we will answer is the following:

**Question 1.** Are there any parameters S, D such that:

$$ \inf \left\{ \tau_{\bar{p}} = \int_{0}^{\infty} x \bar{p}(x) dx; (18) \text{ admits solutions with } \mathrm{Re}(\lambda) \ge 0, \text{ for the distribution } \bar{p} \right\} < \tau_{+} $$

where $\tau_+$ is defined by lemma 3?

**Theorem 2.** Let $S_{max}$ be the point of intersection of the line $y = S$ and the function $g(S) = \frac{\pi^2}{32}(-S^3 - S^2 + S + 1)$. Then for all $S \in (0, S_{max})$ there exists a D and a distribution X such that $E(X) < \tau(S, D)$ and equation (18) is satisfied for some $i\omega$.

*Proof.* Let us write

$$ \begin{cases} x(\omega) = 1 - \frac{\omega^2(S+1)}{D} \\ y(\omega) = \frac{\omega^3 - \omega(S+D)}{D} \end{cases} $$

The point $(x, y)$ is inside the unit disk, so we can write it as

$$ x + iy = \sqrt{x^2 + y^2} e^{i\beta(\omega)} $$

with $\beta \in [0, 2\pi[$.

**Lemma 4.** Let $\omega \in (\omega_-, \omega_+)$ and consider the argument

$$ \theta'(\omega) = \pi - 2 \arctan \left( \frac{y}{1-x} \right), $$

the number

$$ q(\omega) = \frac{y^2 + (1-x)^2}{2(1-x)} \in [0, 1] $$

and the distribution :

$$ \bar{p}_{\omega} = q(\omega)\delta_{\frac{\theta'}{\omega}} + (1 - q(\omega))\delta_0. $$

Then, for every $\omega \in ]\omega_-, \omega_+[$ (or in $]0, \omega_+[$ if $\omega_-$ is not defined), we have

$$ \frac{P(i\omega)}{D} = \int_{0}^{\infty} e^{-i\omega x} \bar{p}_{\omega}(dx). $$

I.e (18) is satisfied with $\lambda = i\omega$ for the distribution $p_\omega$.
**Proof of the lemma**

We notice that for $\omega \in (\omega_-, \omega_+)$ we have $x^2 + y^2 \le 1$ by construction.
We easily check

$$
\begin{cases}
\cos(\theta') = \frac{y^2 - (1-x)^2}{y^2 + (1-x)^2} \\
\sin(\theta') = \frac{2y(1-x)}{y^2 + (1-x)^2}
\end{cases}
$$

Finally, if we denote $E_{\bar{p}_{\omega}}$ the mean with respect to the distriubtion $\bar{p}_{\omega}$, we
obtain then

$$
\begin{cases}
E_{\bar{p}_{\omega}}(\cos(\omega X)) = x(\omega) \\
E_{\bar{p}_{\omega}}(\sin(\omega X)) = y(\omega)
\end{cases}
$$

which is necessary and sufficient condition for (18) to be satisfied for $\lambda = i\omega$,
$\omega \neq 0$.
□

Obviously, we have $p(\omega_+) = 1$, $\theta'(\omega_+) = \beta(\omega_+)$, $\frac{q(\omega_+)\theta'(\omega_+)}{\omega_+} = \tau(\omega_+)$. The core of the proof lies in the following remark:

The core of the proof lies in the following remark:

**Lemma 5.** Let us denote $\tau_{\bar{p}}(\omega) = \int_{0}^{\infty} x \bar{p}_{\omega}(dx)$. Then, for every $S \in (0, S_{\max})$ there exists a $D$ such that the following holds:

$$
\frac{d}{d\omega} \Big|_{\omega=\omega_+} \tau_{\bar{p}}(\omega) > 0.
$$

Indeed, in this situation, $E_{\bar{p}_\omega}(X) < \tau_+$ for $\omega < \omega_+$ close enough to $\omega_+$,
and Question 1 has an affirmative answer.

*Proof.* Let us differentiate $\tau_{\bar{p}}(\omega)$.

$$
\frac{d}{d\omega} \tau_{\bar{p}}(\omega) = \theta' \frac{d}{d\omega} \frac{q}{\omega} + \frac{q}{\omega} \frac{d\theta'}{d\omega}(\omega).
$$

We compute the terms separately

$$
\theta'(\omega) = \pi - 2 \arctan \left( \frac{y}{1-x} \right) = \pi - 2 \arctan \left( \frac{\omega^2 - (S+D)}{\omega(S+1)} \right).
$$

So that,

$$
\begin{align*}
\frac{d\theta'}{d\omega}(\omega) &= -2 \frac{1}{1 + \frac{(\omega^2 - (S+D))^2}{\omega^2(S+1)^2}} \left( \frac{1}{S+1} + \frac{S+D}{\omega^2(S+1)} \right) \\
&= -2 \frac{(S+1)(\omega^2 + S + D)}{\omega^2(S+1)^2 + (\omega^2 - (S+D))^2}.
\end{align*}
$$
For the second term, we notice

$$ \frac{q}{\omega} = \frac{(\omega^3 - \omega(S+D))^2 + \omega^4(S+1)^2}{2\omega^2(S+1)D\omega} = \frac{(\omega^2 - (S+D))^2 + \omega^2(S+1)^2}{2(S+1)D\omega}. $$

Therefore

$$ \begin{aligned} \frac{d}{d\omega} \frac{q}{\omega} &= \frac{1}{2D(S+1)\omega^2} (2\omega^2 (2(\omega^2 - (S+D)) + (S+1)^2) - ((\omega^2 - (S+D))^2 + \omega^2(S+1)^2) \\ &= \frac{1}{2D(S+1)\omega^2} (3\omega^4 + \omega^2((S+1)^2 - 2S - 2D) - (S+D)^2). \end{aligned} $$

We end up with the following:

$$ \frac{d}{d\omega} \tau_{\bar{p}}(\omega) |_{\omega_+} = \theta' \left( \frac{1}{2D(S+1)\omega^2} (3\omega^4 + \omega^2((S+1)^2 - 2S - 2D) - (S+D)^2) \right) \\ - \frac{p}{\omega} 2 \frac{(S+1)(\omega^2 + S + D)}{\omega^2(S+1)^2 + (\omega^2 - (S+D))^2}. $$

Since we simply want an evaluation at $\omega = \omega_+$ and $p(\omega_+) = 1$ by construction, the formula is then a (little) simplified.

$$ \frac{d}{d\omega} \tau_{\bar{p}}(\omega) |_{\omega_+} = \theta' \left( \frac{1}{2D(S+1)\omega^2} (3\omega^4 + \omega^2((S+1)^2 - 2S - 2D) - (S+D)^2) \right) \\ - \frac{2}{\omega} \frac{(S+1)(\omega^2 + S + D)}{\omega^2(S+1)^2 + (\omega^2 - (S+D))^2}. $$

Finally, by definition of $\omega_+$, we notice that for $D = \frac{S(1+S)}{1-S}$, we have the equality $\omega_+^2 = S + D$ and thereby $y(\omega_+) = 0$, so that $\beta(\omega_+) = \theta'(\omega_+) = \pi$. If we are in this situation, the computation are simplified a lot more:

$$ \begin{aligned} \frac{d}{d\omega} \tau_{\bar{p}}(\omega) |_{\omega_+} = \pi & \left( \frac{1}{2D(S+1)(S+D)} (3(S+D)^2 + (S+D)((S+1)^2 - 2S - 2D) - (S+D)^2) \right) \\ & - \frac{2}{\sqrt{(S+D)}} \frac{(S+1)((S+D)+S+D)}{(S+D)(S+1)^2} \\ = \pi & \frac{(S+1)}{2D} - \frac{4}{(S+1)\sqrt{S+D}} \end{aligned} $$

Therefore the derivative has the sign of

$$ \pi - \frac{8D}{(S+1)^2 \sqrt{S+D}} = \pi - \frac{4\sqrt{2S}}{\sqrt{1-S}(S+1)}. $$
This quantity is positive iff $\frac{4\sqrt{2S}}{\sqrt{1-S}(S+1)} < \pi$,

iff $32\frac{S}{1-S} < \pi^2(S^2 + 2S + 1)$, iff $(S^2 + 2S + 1)(1 - S) > \frac{32}{\pi^2}S$.

As $g(S) = \frac{\pi^2}{32}(-S^3 - S^2 + S + 1)$ is decreasing on $(0, 1)$, there is unique point of intersection $S_{max}$ with the bisector. On $(0, S_{max})$, $g$ is above the bisector, so $S_{max}$ defines the interval we are looking for. □

### 3.4 Behavior of the system using realistic parameters

To test our results, we recall that the case of constant coefficients, leading to an exponential distribution for the steady state, has been estimated on patients. We refer to [2]. This gives us access to realistic parameters. Those parameters were estimated under treatment but we can modify (one of) them to have access to the non treated case. From the PDE system 2 we can obtain the ODE system 1 by setting:

$$\mu_c = \mu_u = \mu, a = a_1, g = 1, h(x) = -d_2,$$

We make the following choice of $S$:

$$S(y) = d \frac{(y - y_{min})(y - y_{max})}{y^2 + u_{min}y_{max}}.$$

and we take the same universal parameters as in [2]:

<table><caption>Table 1: Universal parameters for the constant coefficient case.</caption><tr><td>K</td><td>41.667</td></tr><tr><td>d<sub>2</sub></td><td>0.0375</td></tr><tr><td>r</td><td>0.00777</td></tr><tr><td>a<sub>1</sub> (without treatment)</td><td>1.350e+05</td></tr></table>

<table><caption>Table 2: Patient dependent parameters as estimated in table 3 of [2]</caption><thead><tr><th>Patient</th><th>d</th><th>μ</th><th>y<sub>min</sub></th><th>y<sub>max</sub></th></tr></thead><tbody><tr><td>1</td><td>0.051</td><td>3.647e-6</td><td>6.610e4</td><td>3.624e5</td></tr><tr><td>2</td><td>0.026</td><td>2.405e-8</td><td>3.831e4</td><td>3.055e5</td></tr><tr><td>3</td><td>0.054</td><td>4.224e-7</td><td>1.617e4</td><td>3.133e5</td></tr><tr><td>4</td><td>0.181</td><td>8.499e-6</td><td>1.206e3</td><td>1.090e4</td></tr><tr><td>5</td><td>0.038</td><td>5.723e-9</td><td>1.841e3</td><td>3.401e4</td></tr><tr><td>6</td><td>0.058</td><td>1.358e-9</td><td>7.143e3</td><td>7.576e4</td></tr></tbody></table>
With these parameters, there is always a positive steady state for (1) satisfying $\bar{y}_2 < y_{min}$. In the presence of a high steady state (which will satisfy $\bar{y}_2 > y_{max}$) this corresponds to a remission steady state.

Therefore, to apply our results, we will work with a value $\bar{I} = \bar{y}_2$, with $g = 1$ and $\mu_u = \mu_c$ to have insight into the effect of the distribution. We give here the relevant values for the patient dependent parameters.

Returning to what we had before normalisation, we obtain the following value for the critical $\tau$ in the case of a Dirac distribution:

$$ \tau_+(R, S, D) = \frac{1}{R} \tau_+ \left( 1, \frac{S}{R}, \frac{D}{R^2} \right). \qquad (19) $$

under the constraint (corresponding to (12) before normalisation):

$$ R^2 \geq S^2 \qquad (20) $$

and (corresponding to (13))

$$ D \geq R^2 \left(1 - \sqrt{1 - \frac{S^2}{R^2}}\right). \qquad (21) $$

Without treatment the remission equilibrium does not lose stability for any patient. We can see from the table 3 that condition 20 is never satisfied in the non treated case.

<table><caption>Table 3: Parameters for patients without treatment: condition 20 is never satisfied.</caption><thead><tr><th>Patient</th><th>I&#773;</th><th>R<sup>2</sup>/S<sup>2</sup></th></tr></thead><tbody><tr><td>1</td><td>6.0928e+04</td><td>1.7578e-06</td></tr><tr><td>2</td><td>3.8292e+04</td><td>0.0615</td></tr><tr><td>3</td><td>1.6053e+04</td><td>8.1206e-06</td></tr><tr><td>4</td><td>1.0176e+03</td><td>7.1767e-12</td></tr><tr><td>5</td><td>1.8408e+03</td><td>0.0012</td></tr><tr><td>6</td><td>7.1428e+03</td><td>0.1349</td></tr></tbody></table>

On the other hand, under treatment we observe destabilization for some patients. To account for treatment, we have to divide $\alpha_1$ by a factor $k_{inh}$ indicating the inhibitory effects. In table 4 we see the values of S, D and critical $\tau$ for the patients under treatment. Whenever $\tau = +\infty$ we have a stable equilibrium. In the case of treatment we do not know whether this destabilization corresponds to relapse or to another behavior (that could be the formation of a limit cycle for example). This is to be considered in some future work.
<table><caption>Table 4: Parameters for patients under treatment (fitted in [2]).</caption><thead><tr><th>Patient</th><th>a<sub>1</sub> (under treatment)</th><th>S</th><th>D</th><th>&tau; (in days)</th></tr></thead><tbody><tr><td>1</td><td>5.4e02</td><td>3.26e-03</td><td>2.4054e-04</td><td>+&infin;</td></tr><tr><td>2</td><td>1.2e02</td><td>1.44e-05</td><td>1.0526e-04</td><td>158.8137</td></tr><tr><td>3</td><td>3.2e02</td><td>3.723e-04</td><td>3.5612e-04</td><td>142.9637</td></tr><tr><td>4</td><td>0.2e02</td><td>2.489e-02</td><td>9.9090e-04</td><td>+&infin;</td></tr><tr><td>5</td><td>0.3e02</td><td>3.6e-06</td><td>2.4709e-04</td><td>100.8270</td></tr><tr><td>6</td><td>0.5e02</td><td>1.447e-06</td><td>3.1349e-04</td><td>89.6290</td></tr></tbody></table>

# 4 Conclusion

We have established the well posedness and the stability chart for a PDE system generalizing the system introduced in [2]. Although the generic behavior for high steady states and disease free steady states is not fundamentally modified (there might be more high steady states due to the maturity structure), the stability of the remission steady state might be strongly affected by the complexity of the maturity structure.

# Appendix A

Existence result:

**Theorem 3.** Let $u_0 \in L^\infty((0, \infty)) \cap L^1((0, \infty))$ and $(c^0, z^0) \in \mathbb{R}_+^2$. Assume $h \in L^\infty(\mathbb{R}^+)$ and $\limsup_{+\infty} h < 0$. Then, there is a weak solution $(c(t), u(x, t), z(t))$ of the system in the space $C([0, \infty), \mathbb{R}^+) \times C([0, \infty), L_+^1(\mathbb{R}^+)) \times C([0, \infty), \mathbb{R}^+)$.

The existence and uniqueness of the solution results from the Banach fixed point theorem. We consider the space

$$X = C([0, T], L_+^1(\mathbb{R}^+))$$

with the distance

$$d(u, v) = \sup_{t \in [0, T]} \|u - v\|_{L^1}$$

that makes the space complete. $T$ will be appropriately chosen to apply the fixed point theorem. We construct the following operator.

$$\begin{aligned} \Phi: &X \to X \\ \bar{u} &\mapsto u \end{aligned} \qquad (22)$$
where $u$ is a solution of the system:

$$ 
\left\{
\begin{aligned}
\bar{I}(t) &= \int_{0}^{\infty} \bar{u}(x, t)dx \\
\frac{dc}{dt} &= rc(1 - \frac{c}{K}) - \mu_c cz \\
\frac{\partial u}{\partial t}(x, t) + g\frac{\partial u}{\partial x}(x, t) &= h(x)u(x, t) - \mu_u u(x, t)z \\
gu(0, t) &= ac(t) \\
\frac{dz}{dt} &= s - S(\bar{I})z \\
u(x, 0) &= u_0(x) \\
c(0) &= c_0 \\
z(0) &= z_0
\end{aligned}
\right.
$$

For the proof of the theorem we will need the following lemma which proves that the operator $\Phi$ is well defined for all $t \in (0, T)$.

**Lemma 6.** Let $T$ be a fixed time. For $0 \le t \le T$ the components of the solution are well defined with bounds that satisfy:

(i) $0 \le c \le B_c$, where $B_c = \max\{K, c_0\}$.

(ii) $u(t, x) \le \begin{cases} u^0(x - gt)e^{\frac{H(x)-H(x-gt)}{g}}, & x - gt > 0, \\ aB_c e^{\frac{H(x)}{g}}, & x - gt \le 0 \end{cases}$

Integrating we obtain a uniform bound for $\int_0^\infty u(x, t)dx$:

$$ \int_0^\infty u(x, t)dx \le B_u < \infty, $$

with $B_u = B_u(I_0, T)$ a constant depending on $I_0 = \int_0^\infty u_0(x)dx$.

(iii) $0 \le z \le B_z(z_0, T)$, with $B_z(z_0, T) = (z^0 + \frac{s}{S^\infty})e^{S^\infty t} - \frac{s}{S^\infty}, S^\infty = \max(0, \sup(-S))$.

*Proof.* For bounds (i) and (iii), we firstly notice that the non negativity is obtained directly from the structure of the system. The upper bounds for (i) and (iii) are respectively consequence of the logistic growth and Gronwall lemma after noticing that $z$ satisfies the equation:

$$ \dot{z} \le s + mz, m = \sup_I S(I) $$

(ii) Non negativity of $u$ is a consequence of the non negativity of $ac$. For the upper bound on $I$, we can simply notice that $u$ satisfies

$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial t}(x, t) + g \frac{\partial u}{\partial x}(x, t) \le h(x)u(x, t) \\
\\
\displaystyle gu(0, t) \le aB_c.
\end{array}
\right.
$$
Therefore, using the characteristics of the transport equation, we have (by
noting $H(x) = \int_{0}^{x} h$),

$$
u(t,x) \le \begin{cases} u^0(x-gt)e^{\frac{H(x)-H(x-gt)}{g}}, & x-gt > 0, \\ aB_c e^{\frac{H(x)}{g}}, & x-gt \le 0 \end{cases}
$$

Therefore, we have the bound:

$$
\int_0^\infty u \le I_0 \sup_x e^{\frac{H(x)-H(x-gT)}{g}} + aB_c \int_0^{gT} e^{\frac{H(x)}{g}} dx.
$$

Remark: We have an upper bound, since the hypothesis on h ensures

$$
\sup_{0 \le x < y} e^{\frac{H(y)-H(x)}{g}} = K_h < +\infty
$$

From the previous lemma and the Cauchy-Lipschitz theorem we get the
existence and uniqueness of z and c for t ∈ (0, T). As z and c are Lipschitz,
they are uniformly continuous and they can be defined on the closed interval
[0, T]. Once z and c are defined, the boundary conditions are independent
of u so we have a solution of the transport equation. The operator Φ is
well defined. We show now that the operator has a unique fixed point. To
establish the contraction property, we consider two functions ¯u₁, ¯u₂ and let
cᵢ, uᵢ, zᵢ, i = 1, 2 be their images through Φ. In the following estimation we
use the Gronwall inequality.

**Lemma 7.** *Given two sources* $\bar{u}_1, \bar{u}_2$ *the following inequality holds on* $[0, T]$,
*where we have denoted* $\delta c = c_1 - c_2$, $\delta u = u_1 - u_2$, $\delta z = z_1 - z_2$,

$$
\frac{d}{dt} |\delta c| + |\delta z| + \|\delta u(t)\|_1 \le C(T) (|\delta c| + |\delta z| + \|\delta u(t)\|_1) + \bar{C}(T) \|\delta \bar{u}(t)\|_1
$$

Since the degradation terms have a negative contribution, we have im-
mediately

$$
\frac{d}{dt} |\delta c| \leq r |\delta c| + \mu B_c |\delta z|,
$$

Similarly, by the mean value theorem

$$
\begin{align*}
\frac{d}{dt} |\delta z| &\le \max(-S)|\delta z| + z(t)|S(\bar{I}_1) - S(\bar{I}_2)| \\
&\le \max(-S)|\delta z| + B_z(z^0, t) \|S'\|_{\infty} \|\bar{u}_1 - \bar{u}_2\|_1
\end{align*}
$$
Finally, multiplying the PDE on $\delta u$ by $\text{sign}(\delta u)$, we get

$$ \begin{cases} \partial_t |\delta u| + g \partial_x |\delta u| \le h(x) \delta u - \mu_u z(t) |\delta u| - \mu_u |u_1| |\delta z|, \\ g |\delta u(t, 0)| \le a |\delta c(t)| \end{cases} $$

And integrating, we obtain:

$$ \frac{d}{dt} \int_{0}^{\infty} |\delta u| \le a |\delta c| + (\sup h) \int_{0}^{\infty} |\delta u| + \mu_u B_u |\delta z|. $$

Putting everything together, we obtain the lemma with:

$$ C(T) = \max(a + r, \mu_c B_c + \max(-S) + \mu_u B_u, \sup h), \\ \bar{C}(T) = B_z(z^0, T) \|S'\|_{\infty} $$

And as a immediate corollary, we have:

**Corollary 2.** Given, $c^0, z^0, u^0$ for any $T \ge 0$, there exists $n$ such that $\Phi^n$ is a contraction. As a consequence the problem admits a unique solution.

From the Gronwall lemma estimate, we have immediately (using that $\delta u(0) = 0,$)

$$ \|\delta u\|_1(t) \le \bar{C}(T) \int_0^t e^{C(T)(t-s)} \|\delta \bar{u}\|_1(s) ds. $$

Then, we recall the classical iteration argument, denoting $\Phi^n$ the $n$ times composition of $\Phi$, and $\delta u^n = \Phi^n(\bar{u}_1) - \Phi^n(u_2)$ we have

$$ \|\delta u^n\|_1(t) \le \frac{(\bar{C}(T)e^{C(T)T}t)^n}{n!} \sup_{[0,T]} \|\delta \bar{u}(t)\|. $$

And finally, we have

$$ d(\Phi^n(\bar{u})\Phi^n(\bar{v})) \le \underbrace{\frac{(\bar{C}(T)e^{C(T)T}T)^n}{n!}}_{<1 \text{ for } n \text{ large enough}} d(\bar{u}, \bar{v}) $$

As a corollary of the Banach fixed-point theorem we deduce that $\Phi$ has a unique fixed point.

As $z$ and $c$ are defined for all $t$ in $[0, T]$, $u$ is also defined for all $t$ in $[0, T]$. Since $T$ can be chosen arbitrarily, this procedure defines a solution on every interval of time, and consequently, we have a global solution. $\square$
**Proof of theorem 1** We write the perturbed system around a steady state $(\bar{c}, \bar{u}, \bar{z})$. If $(c(t), u(x, t), z(t))$ is a solution, then the perturbation $(\delta c, \delta u, \delta z) = (c - \bar{c}, u - \bar{u}, z - \bar{z})$ satisfies:

$$ egin{cases} 
\frac{d\delta c}{dt} = (r - \frac{2r\bar{c}}{K} - \mu_c\bar{z})\delta z - \mu_c\bar{c}\delta c - \mu_c\delta c\delta z, \\ 
\frac{\partial \delta u}{\partial t}(x, t) + g\frac{\partial \delta u}{\partial x}(x, t) = h(x)\delta u(x, t) - \mu_u\bar{z}\delta u(x, t) - \mu_u\bar{u}(x)\delta z(t) - \mu_u\delta u\delta z, \\ 
g\delta u(0, t) = a\delta c(t) \\ 
\frac{dz}{dt} = -S(\bar{I})\delta z - \bar{z}S'(\bar{I})\delta I \\ 
(-S(\bar{I} + \delta I) + S(\bar{I}) + S'(\bar{I})\delta I)\bar{z} + (-S(\bar{I} + \delta I) + S(\bar{I}))\delta z 
\end{cases} $$

Therefore, if $L$ is the linear part of the system, we can write the equation in the form:

$$ \dot{Y} = LY + \omega(t) $$

where we have written

$$ Y(t) = \begin{pmatrix} c(t) \\ u(t, \cdot) \\ z(t) \end{pmatrix}, $$

$$ \omega(t, x) = \begin{pmatrix} -\mu_c \delta c(t) \delta z(t), \\ -\mu_u \delta u(t, x) \delta z(t), \\ (-S(\bar{I} + \delta I(t)) + S(\bar{I}) + S'(\bar{I}) \delta I(t)) \bar{z} + (-S(\bar{I} + \delta I(t)) + S(\bar{I})) \delta z(t) \end{pmatrix} $$

A few important remarks:

* We restrict ourselves to admissible perturbations (namely we impose initially $(\bar{c} + \delta c(0), \bar{u} + \delta u(0, \cdot), \bar{z} + \delta z(0)) \ge 0$ and in the appropriate space),

* we have already established (in the existence proof) the uniform bound of $I$ (and thereby of $\bar{I} + \delta I$) so that in a neighborhood of the steady state we can assume there exists a constant, independent of the perturbation, such that

$$ \bar{I} + \delta I(t) \le B $$

* this ensures, by Taylor expansion, the existence of a constant $K$ independent of the perturbation such that, if we start in the latter neighborhood of the steady state, we have

$$ |(S(\bar{I}+\delta I(t))-S(\bar{I})-S'(\bar{I})\delta I(t))\bar{z}-(S(\bar{I}+\delta I(t))-S(\bar{I}))\delta z(t)| \le K(|\delta I|^2+|\delta z|^2) $$
From all these remarks, we infer the existence of a constant still denoted by K, such that

$$ \|\omega(t)\|_{X} \leq K \|(\delta c, \delta u, \delta z)\|_{X}^{2} $$

Finally, if we denote $T_L$ the semigroup associated to the linear operator above, we can derive from Duhamel's formula

$$ (\delta c, \delta u, \delta z)(t) = T_L(t)(\delta c, \delta u, \delta z)(0) + \int_0^t T_L(t-s)\omega(s)ds. $$

In terms of norms, this can be written as:

$$ \|(\delta c, \delta u, \delta z)(t)\| = \|T_L(t)\| \|(\delta c, \delta u, \delta z)(0)\| + K \int_0^t \|T_L(t-s)\| \|(\delta c, \delta u, \delta z)(s)\|^2 ds. $$

Since the linear exponential stability of the steady state is characterized by the existence of two positive constants $M, \lambda_0$ such that $\|T_L(t)\| \leq Me^{-\lambda_0 t}$, we have, denoting $y = \|(\delta c, \delta u, \delta z)(t)\|$,

$$ y(t) \leq Me^{-\lambda_0 t}y_0 + KM \int_0^t e^{-\lambda_0(t-s)}y^2(s)ds, $$

which leads to

$$ y(t)e^{\lambda_0 t} \leq My_0 + KM \int_0^t e^{-\lambda_0 s}(y(s)e^{\lambda_0 s})^2 ds, $$

For $y_0 \leq \frac{\lambda_0}{4M^2K}$, we conclude by bootstrap argument that:

$$ y(t)e^{\lambda_0 t} \leq \frac{1}{2} \left( \frac{\lambda_0}{MK} + \sqrt{\left(\frac{\lambda_0}{MK}\right)^2 - 4\frac{\lambda_0 y_0}{K}} \right) $$

which leads to the conclusion.

# 5 Appendix B

## Spectral theory

We have used the following elements of spectral theory whose proofs can be found in [10]. Next to each result we give the page where one can find it in [10].

**Definition 2.** For a closed operator $A: D(A) \subseteq X \to X$ we define:
a) The point spectrum: (p.241)

$$P_{\sigma}(A) = \{\lambda \in \mathbb{C} : \lambda - A \text{ is not injective}\}$$

All $\lambda \in P_{\sigma}(A)$ are called the eigenvalues of A and every $x \neq 0 \in D(A)$ that satisfies $(\lambda - A)x = 0$ is an eigenvector of A associated to the eigenvalue $\lambda$.

b) The approximate point spectrum of A: (p.242)

$$A_{\sigma}(A) = \{\lambda \in \mathbb{C} : \lambda - A \text{ is not injective or } \mathrm{rg}(\lambda - A) \text{ is not closed in } X\}$$

All $\lambda \in A_{\sigma}(A)$ are called approximate eigenvalues of A.

c) The residual spectrum of A: (p.243)

$$R_{\sigma}(A) = \{\lambda \in \mathbb{C} : \mathrm{rg}(\lambda - A) \text{ is not dense in } X\}$$

**Lemma 8.** (p.242)

For a closed operator $A : D(A) \subset X \to X$ and a complex number $\lambda \in \mathbb{C}$ we have that: $\lambda \in A_{\sigma}(A)$ if and only if there is a sequence $(x_n)_{n \in \mathbb{N}} \subset D(A)$, called approximate eigenvector, such that $||x_n|| = 1$ and $\lim_{n \to \infty} ||Ax_n - \lambda x_n|| = 0$.

**Theorem 4.** *Spectral theorem for point and residual spectrum* (p.277)

For the generator $(A, D(A))$ of a strongly continuous semi group $(T(t))_{t \ge 0}$ defined on a Banach space, we have the following identities:

$$P_{\sigma}(T(t))\{0\} = e^{tP_{\sigma}(A)}$$

$$R_{\sigma}(T(t))\{0\} = e^{tR_{\sigma}(A)}$$

for every $t \ge 0$.

**Definition 3.** (p.250) Let $A : D(A) \subset X \to X$ be a closed operator. Then,

$$s(A) = \sup\{\Re(\lambda) : \lambda \in \sigma(A)\}$$

is the spectral bound of A.

Moreover,

$$\omega_0 = \inf\{\omega \in \mathbb{R} | \text{ such that there exists } M_{\omega} \ge 1 : ||T(t)|| \le M_{\omega}e^{\omega t} \text{ for every } t \ge 0\}$$

is the growth bound of A.
**Proposition 10.** (p.251) For the growth and the spectral bound of an operator A that is the infinitesimal generator of a semi group $(T(t))_{t \ge 0}$ we have:

$$ 
\begin{aligned}
-\infty \le s(A) \le \omega_0 &= \inf_{t>0} \frac{1}{t} \log ||T(t)|| = \lim_{t \to \infty} \frac{1}{t} \log ||T(t)|| \\
&= \frac{1}{t_0} \log r(T(t_0)) < \infty
\end{aligned}
$$

for every $t_0 > 0$. In particular, the spectral radius of $T(t)$ is given by

$$ r(T(t)) = e^{\omega_0 t}, $$

for every $t \ge 0$.

**Proposition 11.** (p.39)

For a strongly continuous semi group $(T(t))_{t \ge 0}$ there are constants $\omega \in \mathbb{R}$ and $M \ge 1$ such that

$$ ||T(t)|| \le Me^{\omega t} $$

for every $t \ge 0$.

# 6 Acknowledgements

We would like to thank Mosatfa Adimy and Samuel Bernard for their constructive discussions and suggestions.

# References

[1] S. Bernard and F. Crauste. Optimal linear stability condition for scalar differential equations with distributed delay. *Discrete & Continuous Dynamical Systems-B*, 20(7):1855, 2015.

[2] A. Besse, G. D. Clapp, S. Bernard, F. E. Nicolini, D. Levy, and T. Lepoutre. Stability Analysis of a Model of Interaction Between the Immune System and Cancer Cells in Chronic Myelogenous Leukemia. *Bulletin of Mathematical Biology*, pages 1–27, 2017.

[3] G. D. Clapp, T. Lepoutre, R. El Cheikh, S. Bernard, J. Ruby, H. Labussière-Wallet, F. E. Nicolini, and D. Levy. Implication of the autologous immune system in bcr-abl transcript variations in chronic myelogenous leukemia patients treated with imatinib. *Cancer Research*, 75(19):4053–4062, 2015.
[4] R. E. Clark, F. Polydoros, J. F. Apperley, D. Milojkovic, K. Rothwell, C. Pocock, J. Byrne, H. de Lavallade, W. Osborne, L. Robinson, S. G. O'Brien, L. Read, L. Foroni, and M. Copland. De-escalation of tyrosine kinase inhibitor therapy before complete treatment discontinuation in patients with chronic myeloid leukaemia (destiny): a non-randomised, phase 2 trial. *The Lancet Haematology*, 6(7):e375–e383, 2019.

[5] M. W. N. Deininger, J. M. Goldman, and J. V. Melo. The molecular biology of chronic myeloid leukemia. *Blood*, 96(10):3343-3356, 11 2000.

[6] A. C. Fassoni, I. Roeder, and I. Glauche. To cure or not to cure: consequences of immunological interactions in cml treatment. *Bulletin of mathematical biology*, 81(7):2345-2395, 2019.

[7] T. Hähnel, C. Baldow, A. C. Fassoni, J. Guilhot, F. Guilhot, S. Saussele, S. Mustjoki, S. Jilg, P. J. Jost, S. Dulucq, F.-X. Mahon, I. Roeder, and I. Glauche. Inferring immunological control mechanisms from tki dose alterations in cml patients. *bioRxiv*, 2019.

[8] M. Ilander, U. Olsson-Strömberg, H. Schlums, J. Guilhot, O. Brück, H. Lähteenmäki, T. Kasanen, P. Koskenvesa, S. Söderlund, M. Höglund, et al. Increased proportion of mature nk cells is associated with successful imatinib discontinuation in chronic myeloid leukemia. *Leukemia*, 31(5):1108-1116, 2017.

[9] F.-X. Mahon, D. Réa, J. Guilhot, F. Guilhot, F. Huguet, F. Nicolini, L. Legros, A. Charbonnier, A. Guerci, B. Varet, G. Etienne, J. Reiffers, and P. Rousselot. Discontinuation of imatinib in patients with chronic myeloid leukaemia who have maintained complete molecular remission for at least 2 years: the prospective, multicentre stop imatinib (stim) trial. *The Lancet Oncology*, 11(11):1029–1035, 2010.

[10] K.-J. E. R. Nagel. *One-parameter semigroups for linear evolution equations*, volume 194. Springer-Verlag New York, 2000.

[11] D. Rea, G. Henry, Z. Khaznadar, G. Etienne, F. Guilhot, F. Nicolini, J. Guilhot, P. Rousselot, F. Huguet, L. Legros, et al. Natural killer-cell counts are associated with molecular relapse-free survival after imatinib discontinuation in chronic myeloid leukemia: the immunostim study. *Haematologica*, 102(8):1368-1377, 2017.
[12] S. Saussele, J. Richter, A. Hochhaus, and F. Mahon. The concept of treatment-free remission in chronic myeloid leukemia. *Leukemia*, 30(8):1638-1647, 2016.

[13] A. Yong and A. Hughes. Immune effector recovery in chronic myeloid leukemia and treatment-free remission. *Frontiers in Immunology*, 8, 04 2017.