Apollos Besse

Geoffrey D Clapp

Samuel Bernard

Franck E Nicolini

Doron Levy

Thomas Lepoutre

Stability Analysis of a Model of Interaction Between the Immune System and Cancer Cells in Chronic Myelogenous Leukemia

We describe here a simple model for the interaction between leukemic cells and the autologous immune response in chronic phase chronic myelogenous leukemia (CML). This model is a simplified version of the model we proposed in [Clapp et al., Cancer Research, 75:4053-4062, 2015]. Our simplification is based on the observation that certain key characteristics of the dynamics of CML can be captured with a three compartments model: two for the leukemic cells (stem cells and mature cells) and one for the immune response. We characterize the existence of steady states and their stability for generic forms of immunosuppressive effects of leukemic cells. We provide a complete co-dimension one bifurcation analysis. Our results show how clinical response to tyrosine kinase inhibitors treatment is compatible with the existence of a stable low-disease, treatment-free steady state.

Introduction

Chronic myeloid leukemia is a clonal disease that arises from a discrete population of hematopoietic stem cells residing in the bone marrow harmed by a single fusion gene, BCR-ABL, encoding for a deregulated tyrosine kinase activity responsible for the disease [START_REF] Deininger | The molecular biology of chronic myeloid leukemia[END_REF]. Unlike most cancers, the measurement of the BCR-ABL transcripts can quantify precisely the tumor load and its evolution throughout time under targeted therapies (i.e. tyrosine kinase inhibitors (TKI)). The first-in-class TKI Imatinib (Glivec R /Gleevec R , Novartis) is a targeted therapy that inhibits the deregulated tyrosine kinase activity of ABL in leukemic cells, and induces within 18 hours cell death by apoptosis [START_REF] Graham | Primitive, quiescent, Philadelphia-positive stem cells from patients with chronic myeloid leukemia are insensitive to STI571 in vitro[END_REF]. This TKI and the other TKIs derived since 2005 from IM revolutionized the paradigms of treatment of CML and transformed this disease from a mostly fatal illness into a life-long indolent pathology [START_REF] O'brien | Imatinib compared with interferon and low-dose cytarabine for newly diagnosed chronic-phase chronic myeloid leukemia The New England[END_REF]. However, in most cases, under TKI, the disease remains detectable at low levels. It is rare when CML becomes undetectable, possibly cured. After stopping targeted therapies, most patients relapse and require a new TKI treatment. Yet, recent clinical trials [START_REF] Mahon | 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[END_REF][START_REF] Preudhomme | Imatinib plus peginterferon alfa-2a in chronic myeloid leukemia The New England[END_REF][START_REF] Rousselot | Imatinib mesylate discontinuation in patients with chronic myelogenous leukemia in complete molecular remission for more than 2 years[END_REF] have shown that a long TKI treatment duration may result with remission after treatment withdrawal, and, in some cases, although detectable levels of transcripts might reappear, patients may end up in an indolent detectable disease not requiring further treatments.

There is increased evidence that such remissions (so-called treatment-free remissions) do not require disease eradication [START_REF] Ross | Patients with chronic myeloid leukemia who maintain a complete molecular response after stopping imatinib treatment have evidence of persistent leukemia by DNA PCR[END_REF]. Such coexistence of low-level leukemic cell load and non-active relapse suggest the occurrence of leukemic cell control mechanisms. In CML, it is well-known that the autologous or allogeneic immune responses against the tumor have a critical importance in disease cure or long-term control. In an autologous setting, complete hematologic and cytogenetic remission, even in accelerated phase CML, has been observed after generating ex-vivo autologous leukemia-reactive cytotoxic T-lymphocyte infusions [START_REF] Falkenburg | Complete remission of accelerated phase chronic myeloid leukemia by treatment with leukemia-reactive cytotoxic T lymphocytes[END_REF]. Nowadays, in successful treatment-free remission setting, elevated counts of natural killer cells [START_REF] Ilander | IFNα induces prolonged remissions modeling curative immunologic responses in chronic myeloid leukemia Oncoimmunology[END_REF] and plasmacytoid dendritic cells [START_REF] Burchert | Frequency of CTLA-4 Receptor Ligand (CD86, B7.2)-Positive Plasmacytoid Dendritic Cells Predicts Risk of Disease Recurrence after Tyrosine-Kinase Inhibitor Discontinuation in Chronic Myeloid Leukemia: Results from a Prospective Substudy of the Euroski Trial[END_REF] at time of TKI withdrawal have been observed, supporting the hypothesis that the autologous immune system might play a critical role in controlling the disease after TKI withdrawal.

Several mathematical models of the dynamics of CML during treatment and after treatment cessation have provided insight into the mechanism of remission and relapse. We briefly comment on some of these models. For a more comprehensive review of mathematical models for CML we refer to [START_REF] Clapp | A Review of Mathematical Models for Lymphoma and Leukemia[END_REF].

Michor et al. used a mathematical model to explain the failure of complete eradication of the disease under TKI treatments by the development of resistance [START_REF] Michor | Dynamics of chronic myeloid leukaemia[END_REF]. Komarova and Wodarz proposed a model that incorporates the stem cell ability to become quiescent and to develop resistance to treatment [START_REF] Komarova | Effect of cellular quiescence on the success of targeted CML therapy[END_REF]. An alternative model was proposed by Roeder et al. [START_REF] Roeder | Dynamic modeling of imatinib-treated chronic myeloid leukemia: functional insights and clinical implications[END_REF]. There, an agent-based model in which the fate of each cell is stochastically determined, according to its characteristics and its environment. The model describes the competition of leukemic stem cells versus disease-free stem cells in the bone marrow, and the effect of TKIs on competition [START_REF] Roeder | Dynamic modeling of imatinib-treated chronic myeloid leukemia: functional insights and clinical implications[END_REF]. The potential long-term effect of immune response was studied by Kim et al. in [START_REF] Kim | Dynamics and potential impact of the immune response to chronic myelogenous leukemia[END_REF] by adding an immune compartment to the model of [START_REF] Michor | Dynamics of chronic myeloid leukaemia[END_REF].

Recently, we proposed a model describing the interactions between leukemic cells and the autologous immune system in IM-treated patients [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF] (see also [START_REF] Clapp | BCR-ABL transcript variations in chronic myelogenous leukemia patiens on imatinib first-line: Possible role of the autologous immune system[END_REF]). The model was motivated by an observation that many patients who otherwise responded well to IM therapy still showed variations in their BCR-ABL transcripts. Our modeling results suggested that, at diagnosis, a patient's leukemic load is able to partially or fully suppress immune response developed in a ma-jority of patients towards the CML clone(s). IM therapy drives the residual leukemic cell population into the "immune window", allowing the patient's autologous immune cells to expand and eventually mount an efficient recognition of the residual leukemic burden. This drives the leukemic load below this immune window, allowing the residual leukemic population to partially recover until another weaker immune response is initiated. Thus, we suspect that the autologous immune response may explain the oscillations in the BCR-ABL transcripts regularly observed in patients on IM.

In this paper we propose a new model for describing the interaction between leukemic cells and the autologous immune system in chronic phase CML. This model is a simplified version of our original model from [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF], in which we reduce the number of equations and simplify some of the interaction terms. Our simplification is based on the observation that certain key characteristics of the dynamics of CML can be captured with a three compartments model: two for the leukemic cells (proliferating stem cells and mature cells) and one for the immune response.

The new model is sufficiently accessible to allow us to conduct mathematical analysis. Indeed, we characterize the existence of steady states and their stability for generic forms of immunosuppressive effects of leukemic cells, and provide a complete co-dimension one bifurcation analysis. Our results show how clinical response to TKI treatment is compatible with the existence of a stable lowdisease, treatment-free steady state.

The structure of this manuscript is as follows. In Section 2 we briefly overview the mathematical model of [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF] after which we introduce our new model [START_REF] Clapp | A Review of Mathematical Models for Lymphoma and Leukemia[END_REF]. The positive steady states of the system (4) are characterized in Section 2. A stability analysis of these steady states in conducted in Section 3. This section includes our main Theorems, Theorem 1 and 2, which describe the bifurcation diagrams of the steady states. Applications of the stability analysis are given in Section 4, and a concluding discussion is provided in Section 6.

The Model

The Clapp et al. model

The mathematical model of Clapp et al. described the dynamics of the interaction between cancer cells and the immune response in CML [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF]. This model is given as a system of five ODEs, describing the dynamics of four states for the leukemic cells (quiescent and cycling stem cells, progenitor cells, and mature cells), and the dynamics of the immune response.

Let y 0 and y 1 represent the concentration of quiescent and cycling leukemic stem cells, respectively, y 2 the concentration of progenitor cells, y 3 the concentration of mature leukemic cells, and z the concentration of active immune cells.

The model of [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF] is written as the following system of ODEs:

                                         y 0 = b 1 y 1 -a 0 y 0 - µy 0 z 1 + εy 2 3 , y 1 = a 0 y 0 -b 1 y 1 + ry 1 1 - y 1 K -d 1 y 1 - µy 1 z 1 + εy 2 3 , y 2 = a 1 inh 1 y 1 -d 2 y 2 - µy 2 z 1 + εy 2 3 , y 3 = a 2 inh 2 y 2 -d 3 y 3 - µy 3 z 1 + εy 2 3 , z = s -dz + αy 3 z 1 + εy 2 3 .
(1)

The model equations [START_REF] Burchert | Frequency of CTLA-4 Receptor Ligand (CD86, B7.2)-Positive Plasmacytoid Dendritic Cells Predicts Risk of Disease Recurrence after Tyrosine-Kinase Inhibitor Discontinuation in Chronic Myeloid Leukemia: Results from a Prospective Substudy of the Euroski Trial[END_REF] describe transitions between leukemic populations (at rates b 1 , a 0 , a 1 and a 2 ). The different populations have associated death rates (d 1 , d 2 , d 3 and d). The cycling stem cells are assumed to grow logistically with growth rate r and carrying capacity K. The last terms in the first four equations in (1) represent the death of leukemic cells caused by an immune response. The mass action term µy i z represents the killing of leukemic cells by the immune system, where µ is the maximal rate (per immune cell) at which an immune cell will engage and kill a leukemic cell.

The last equation in [START_REF] Burchert | Frequency of CTLA-4 Receptor Ligand (CD86, B7.2)-Positive Plasmacytoid Dendritic Cells Predicts Risk of Disease Recurrence after Tyrosine-Kinase Inhibitor Discontinuation in Chronic Myeloid Leukemia: Results from a Prospective Substudy of the Euroski Trial[END_REF] represents the concentration of autologous immune cells. The first term, s, is a constant source term for the immune cells, who die at rate d. The mass action term αy 3 z represents the expansion (proliferation) of the immune cell pool in response to its leukemia stimulus, which occurs with maximal rate per leukemic cell. Only the contributions of the mature leukemic cells y 3 to immune stimulation are included since they are a much larger population than the immature leukemic cells.

This model is based on assuming that immunosuppression acts in two ways. First, mature leukemic cells inhibit the expansion of immune cells. In the last equation, the immune cell expansion term αy 3 z is divided by 1 + εy 2 3 , where the constant ε determines the strength of the immunosuppression. Second, mature leukemic cells are assumed to decrease the killing capacity µ of activated immune cells, also by a factor of 1 + εy 2 3 . By implementing immunosuppression in this way, Clapp et al. encoded an autologous immune response that is effective only with intermediate levels of leukemic cells. When the leukemic load is small, only a small number of immune cells is stimulated to respond. On the other hand, although large leukemic loads provide a stronger stimulus, the leukemic cells are able of suppressing the efficacy of the immune system. Thus, the immune response will be negligible when the leukemic load is either very small, at levels undetectable by the immune system, or very large, at levels that overwhelm and suppress the immune system. A strong immune response can occur only when the leukemic load y 3 is at an intermediate level, within a range [y min , y max ], the so-called immune window. The immune window is defined as the range of y 3 for which the rate of immune stimulation exceeds the death rate. The parameters inh 1 and inh 2 represent the inhibition of leukemic cell amplification by TKI treatment. TKI therapy may be used to drive the leukemic load into this immune window, allowing the autologous immune system to assist the drug in the elimination of the residual leukemic cells.

Introducing a simplified model

In [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF] we provided no theoretical analysis of the model [START_REF] Burchert | Frequency of CTLA-4 Receptor Ligand (CD86, B7.2)-Positive Plasmacytoid Dendritic Cells Predicts Risk of Disease Recurrence after Tyrosine-Kinase Inhibitor Discontinuation in Chronic Myeloid Leukemia: Results from a Prospective Substudy of the Euroski Trial[END_REF]. Analysis of the model in its full generality is beyond the reach of this paper. Instead, we introduce in this subsection a reduced model which captures some of the essential features of the original model.

Our main simplification is to reduce the original system to three ODEs, keeping only two states for the leukemic cells. In addition, we assume that immunosuppression does not directly suppress the killing capacity µ of immune cells. These assumptions make the theoretical analysis accessible.

Let y 1 represent the concentration of cycling stem cells, y 2 the concentration of mature leukemic cells, and z the concentration of immune cells. We consider the following system of ODEs:

               y 1 = ry 1 1 - y 1 K -d 1 y 1 -µy 1 z, y 2 = a 1 k inh y 1 -d 2 y 2 -µy 2 z, z = s -f (y 2 )z. (2)
We note that the system (2) is a simplified version of system [START_REF] Burchert | Frequency of CTLA-4 Receptor Ligand (CD86, B7.2)-Positive Plasmacytoid Dendritic Cells Predicts Risk of Disease Recurrence after Tyrosine-Kinase Inhibitor Discontinuation in Chronic Myeloid Leukemia: Results from a Prospective Substudy of the Euroski Trial[END_REF]. We also note that y 2 in the system (2) plays the role of y 3 in the system (1). The parameter a 1 is the product of the stem cell differentiation rate and the amplification factor between stem and mature cells. The effect of the treatment is represented by dividing a 1 by an inhibition coefficient k inh . The parameters r and K are the cycling stem cells growth rate and carrying capacity, respectively. The natural death rates of each population are given by d, d 1 , d 2 , and µ is the probability that an immune cell will interact with a leukemic cell and kill it.

The function f represents the net suppression effect of the leukemic cells on the immune system. It mixes the natural death rate of immune cells d = f (0), the competitive (immune response) and the immunosuppressive effect of cancer cells. The suppression rate f can take negative values; in this case, immune response outbalances immunosuppression.

In all what follows, we assume

f ∈ C 2 (R + ),
and

∃Y > 0, f < 0, on [0, Y ), f > 0, on (Y, +∞) (3)
In other words, we assume f to be strictly decreasing for low amount of cancer cells and strictly increasing for high amount of cancer cells (its monotonicity changes exactly once). An immediate important consequence of this hypothesis is the following result:

Proposition 1 (Zeros of f ) Let f ∈ C 2 (R + ) satisfy (3), then either 1. f has no zero on R + , 2. f has one zero y min ∈ (0, Y ],
3. f has two zeros y min ∈ (0, Y ) and y max ∈ (Y, +∞)

The notations used in this proposition will be kept throughout the manuscript.

In order to eliminate the parameter d 1 from (2), we make a change of variables r new ← r old -d 1 and K new ← (r new /r old )K old . We also set a 1 ← a 1 /k inh . This leads to the following system:

             y 1 = ry 1 1 - y 1 K -µy 1 z, y 2 = a 1 y 1 -d 2 y 2 -µy 2 z, z = s -f (y 2 )z. ( 4
)
3 Steady states for the simplified model

We want to determine the positive steady states of the system (4). We start with a simple characterization.

Proposition 2

The system (4) always has a unique, disease-free non-negative steady state (0, 0, s f (0) ). Other positive steady states are given by the solutions to the equations

ȳ1 = (r + d 2 )ȳ 2 a 1 + r K ȳ2 , (5) z
= s f (ȳ 2 ) . (6)
Here, ȳ2 is a positive zero of the function P given by

P (X) = f (X)(X -M ) + µs d 2 M d 2 r + X , (7)
with M = a1K d2 , such that f (ȳ 2 ) > 0. For any positive steady state, we have

0 < ȳ2 ≤ M. (8)
Proof. The steady states of (4) are the solutions of the nonlinear system:

             0 = ry 1 1 - y 1 K -µy 1 z, 0 = a 1 y 1 -d 2 y 2 -µy 2 z, 0 = s -f (y 2 )z. (9)
The point (0, 0, s f (0) ) is always a steady state. It corresponds to the disease-free steady state, where the leukemic cell populations vanish and the immune system is maintained at a basal, surveillance level s f (0) . The positivity of all coefficients provides bounds for ȳ1 and ȳ2 :

ȳ1 ≤ K, ȳ2 ≤ a 1 K d 2 := M.
Other steady states can be found by:

             µz = r 1 - y 1 K , µy 2 z = a 1 y 1 -d 2 y 2 , f (y 2 )z = s, (10)
which leads to:

ȳ1 = (r + d 2 ) ȳ2 a 1 + r K ȳ2 , (11) and z
= s f (y 2 ) . ( 12
)
As we are looking for non-negative steady states, equation [START_REF] Kim | Dynamics and potential impact of the immune response to chronic myelogenous leukemia[END_REF] implies that f (ȳ 2 ) must be non-negative. Since the steady state values ȳ1 and z are uniquely determined by ȳ2 , and we reduce the problem to an equation on ȳ2 alone. For ȳ1 , we have

1 - ȳ1 K = M -ȳ2 M + r d2 ȳ2
.

Then, equations ( 10) and ( 12) imply:

µs f (ȳ 2 ) = r M -ȳ2 M + r d2 ȳ2
, which means that,

µs r M + r d 2 ȳ2 = f (ȳ 2 )(M -ȳ2 ).
This leads to the desired equation for ȳ2 :

f (ȳ 2 )(ȳ 2 -M ) + µs d 2 M d 2 r + ȳ2 = 0.
The following lemma provides the location of the steady states with respect to the zeros of f . Lemma 1. Keeping the notations of Proposition 1, positive steady states of the system satisfy ȳ2 ∈ (0, y min ], or ȳ2 ∈ [y max , +∞).

Let (C 1 ) denote the condition

(C 1 ) : d > µs r .
Proposition 3 Assume (C 1 ) is satisfied, there exists a disease steady state for (4) such that ȳ2 ∈ (0, M ). Moreover, if f vanishes and (with notation of Proposition 1) y min < M , there is a unique steady state for (4) that satisfies

0 < ȳ2 < y min < M.
All other steady states satisfy y 2 ≥ y max and y 2 < M .

Proof. Assuming (C 1 ), we have

P (0) = -d - µs r M < 0, P (M ) = µs d 2 d 2 r + 1 M > 0.
Thus, P changes sign in the interval (0, M ), which implies, from Proposition 2 that there is a steady state with ȳ2 ∈ (0, M ). Moreover, f (ȳ 2 ) must be positive, so the solution for ȳ2 corresponds to a positive steady state, according to Proposition 2.

Keeping the notations of proposition 1, assume y min ∈ (0, M ). Because f (0) = d is positive, f (y 2 ) stays positive for y 2 < y min . From the definition of P in (7), we have P (y min ) > 0. Since P (0) < 0 and since y min < M , P admits at least one zero ȳ2 ∈ (0, y min ).

To prove uniqueness, we highlight the fact that P is monotonous on [0; y min ]:

P (X) = f (X) + f (X)(X -M ) + θ
is positive on [0; y min ], since function f stays positive and strictly decreasing. Moreover, at the unique zero ȳ2 < y min , f (ȳ 2 ) is still positive, and Proposition 2 ensures that ȳ2 corresponds to a positive steady state x = (ȳ 1 , ȳ2 , z).

When (C 1 ) is satisfied, the smallest disease steady state (in terms of y 2 ) is denoted as the low disease steady state. When d = µs r , P (0) = 0, and the low disease steady state becomes disease-free. For d ≤ µs r , the steady state becomes negative, and does not exists anymore.

In order to draw a complete picture of the dynamics of the system (4), it is important to characterize the stability of the steady states. Because the system is bounded (positive cell concentrations cannot grow unbounded), we expect the system to converge either to a stable steady state or perhaps to a nonsteady-state bounded attractor, such as a limit cycle. We will show that under biologically relevant assumptions, there always exists an asymptotically stable steady state, and that none of the steady states can be destabilized by a Hopf bifurcation. We are able to provide a complete local co-dimension 1 bifurcation diagram. Although more complex bifurcation diagrams are possible, they are non generic as they occur only for specific combinations of two or more parameters.

Stability of the disease-free steady state

In the disease-free steady state (0, 0, s d ), where d = f (0), there are no leukemic cells, but the concentration of immune cells is positive, due to a basal cell production level s.

Proposition 4

The disease-free steady state (0, 0, s d ) is asymptotically stable if and only if r < µs d . If r = µs d , then this disease-free steady state is locally attractive for positive solutions.

Proof. The linearized system around the disease-free steady state is

                 y 1 = ry 1 - µs d y 1 , y 2 = a 1 y 1 -d 2 y 2 - µs d y 2 , z = -dz + s d f (0)y 2 . ( 13
)
The disease-free Jacobian matrix is

J 0 =     r -µs d 0 0 a 1 -d 2 -µs d 0 0 s d f (0) -d     . (14)
The matrix ( 14) is triangular, and its eigenvalues lie on the diagonal, which are all real coefficients. Two eigenvalues are always strictly negative. The third one, given by r -µs d , can be either negative or positive. Therefore, the disease-free steady state is asymptotically stable if rd < µs, and it is unstable if rd > µs (condition (C 1 )).

In addition, if rd = µs, it is easy to verify that for ε > 0 small enough,

D ε = [0; ε]×[0; a1 d2 ε]×[z -r Kµ ;
+∞) is a stable set for the equations. For solutions starting in D ε , since f is non-increasing on [0; a1 d2 ε] and since z ≥ s -f (0)z, we have

lim inf +∞ z ≥ z.
Finally, since y 1 = -ry 2 1 /K + µy 1 (z -z), we have lim sup

+∞ y 1 ≤ Kµ r lim sup +∞ (z -z) ≤ 0.
Hence, due to the non-negativity of y 1 , y 1 (+∞) = 0, and using the equations ( 4), we conclude that (y 1 (t), y 2 (t), z(t)) → (0, 0, z).

From now on, we assume that the disease-free steady state is unstable, and we want to know if there are other stable steady states.

Stability of the disease steady states

We denote by (ȳ 1 , ȳ2 , z) a non-negative steady state of system (4), with ȳ1 and/or ȳ2 positive. The Jacobian matrix is

J =     r -2 r K ȳ1 -µz 0 -µȳ 1 a 1 -d 2 -µz -µȳ 2 0 -f (ȳ 2 )z -f (ȳ 2 )     . (15)
Using [START_REF] Ilander | IFNα induces prolonged remissions modeling curative immunologic responses in chronic myeloid leukemia Oncoimmunology[END_REF] the matrix (15) becomes

J =     -r K ȳ1 0 -µȳ 1 a 1 -a 1 ȳ1 ȳ2 -µȳ 2 0 -f (ȳ 2 )z -f (ȳ 2 )     . (16)
In order to determine the stability of the steady state, we compute the characteristic polynomial χ

J χ J (X) = X + r K ȳ1 (X + a 1 ȳ1 ȳ2 )(X + f (ȳ 2 )) -µȳ 2 zf (ȳ 2 ) -a 1 µȳ 1 zf (ȳ 2 ) = X 3 + r K ȳ1 + a 1 ȳ1 ȳ2 + f (ȳ 2 ) X 2 + r K ȳ1 a 1 ȳ1 ȳ2 + f (ȳ 2 ) + a 1 ȳ1 ȳ2 f (ȳ 2 ) -µȳ 2 zf (ȳ 2 ) X + r K ȳ1 a 1 ȳ1 ȳ2 f (ȳ 2 ) -µȳ 2 zf (ȳ 2 ) -a 1 µȳ 1 zf (ȳ 2 ) ,
and using the relation r

K ȳ1 + a 1 ȳ1 ȳ2 = r + d 2 , χ J (X) =X 3 + X 2 (r + d 2 + f (ȳ 2 )) + (r + d 2 )f (ȳ 2 ) -µz ȳ2 f (ȳ 2 ) + r K ȳ1 a 1 ȳ1 ȳ2 X + r K ȳ1 a 1 ȳ1 ȳ2 f (ȳ 2 ) -(r + d 2 )µz ȳ2 f (ȳ 2 ) . (17)
In what follows we will again keep notation of Proposition 1. In case 3 of proposition 1, we will distinguish the low steady state (ȳ 2 ≤ y min ) and the others (ȳ 2 ≥ y max ).

Stability of the low disease steady state

We are interested in proving the stability of the disease steady state given by Proposition 3 under condition (C 1 ). It is the only disease steady state ( ȳ1 , ȳ2 , z) such that ȳ2 ≤ y min .

Proposition 5 Assume the low disease steady state of (4) exists. It is asymptotically stable if and only if the disease-free steady state is unstable, i.e., if and only if condition (C 1 ) is satisfied. For rd = µs, there is a transcritical bifurcation between the disease-free steady state and the low disease steady state.

Proof. We can assume that the disease-free steady state is unstable, which implies that rd > µs and the existence of the low disease steady state. Either f has no zero in (0, M ), or f admits a smallest zero y min ∈ (0, M ). In either case, at the low disease steady state, we have f (ȳ 2 ) < 0. This means that the coefficients of χ J are all positive. We need to prove that χ J verifies the condition of Lemma 4 (in Appendix).

a χ b χ -c χ = (r + d 2 + f (ȳ 2 )) (r + d 2 )f (ȳ 2 ) -µz ȳ2 f (ȳ 2 ) + r K ȳ1 a 1 ȳ1 ȳ2 - r K ȳ1 a 1 ȳ1 ȳ2 f (ȳ 2 ) -(r + d 2 )µz ȳ2 f (ȳ 2 ) = (r + d 2 ) (r + d 2 )f (ȳ 2 ) + r K ȳ1 a 1 ȳ1 ȳ2 +f (ȳ 2 ) ((r + d 2 )f (ȳ 2 ) -µz ȳ2 f (ȳ 2 )) > 0.
Here we used the fact that f (ȳ 2 ) < 0.

By Lemma 4, we know that the dominating root of χ J has a negative real part, which means that the dominating eigenvalue of J has a negative real part. Hence the low steady state is asymptotically stable.

For rd = µs, the low disease steady state crosses the disease-free steady state. At this point, the fusion steady state is locally stable for positive solutions (as seen in Subsection 4.1). For rd < µs, it becomes negative, and is necessarily unstable, i.e., it is a transcritical bifurcation.

The asymptotic stability of the low disease steady state corresponds to a state where the immune system is able to keep the leukemic cell population at a low level, even without treatment. However, we demonstrate numerically that the basin of attraction of low disease steady state is small, and solutions typically converge to the large disease steady state.

Stability of the other steady states

The existence of steady states other than the disease-free and the low disease steady states depends on the exact shape of the leukemia-induced immunosuppression function f . We assume that we are in the case 3 of Proposition 1 and moreover y max < M . We denote this condition as (C 2 ).

(C 2 ): f has two zeros y min and y max such that 0 < y min < y max < M .

The assumption about the two zeros follows from the biological assumption that there exists a range of leukemic cell concentrations for which the immune system is able to mount an immune response (when f is negative, the growth rate of z is strictly positive). If M ≤ y max , the only possible disease steady state is the low disease steady state, which is asymptotically stable under (C 1 ). From now on, we assume that M > y max , opening the door for other disease steady states to exist. All these disease steady states verify f ( ȳ2 ) > 0.

We give a function that does not depend on s and µ, that gives the steady states and their stability. We denote

A(X) = f (X)(M -X), B(X) = X + d2 r M , θ = µs/d 2 and θ max = sup x∈(ymax,M ) { A B (x)}.
All disease steady states are entirely determined by their component ȳ2 defined as a zero of the function P defined in [START_REF] Deininger | The molecular biology of chronic myeloid leukemia[END_REF]. We express P as a parametric linear combination of two functions A and B, P θ (X) = -A(X) + θB(X).

(

) 18
Theorem 1. If the function A B admits only one critical point, then the following holds:

1. For all θ in [0, θ max ) there exist two disease steady states x(1) and x(2) with

y max ≤ ȳ(1) 2 < ȳ(2)
2 ≤ M . The disease steady state x(1) is unstable and x(2) is asymptotically stable. 2. For all θ > θ max there is no steady state with ȳ2 in (y max , M ). 3. There is an saddle-node bifurcation for θ = θ max , the corresponding double steady states is generically unstable.

Proof. The proof is constructed as follows:

-We give a necessary and sufficient condition for stability, that is based on condition (C 2 ) (Lemma 1). -We connect this condition to the functions A and B, which are independent of the parameter θ (Lemma 2).

branch could co-exist with the ones we highlight. Another consequence is that, at ȳ2 fixed, the stability is known regardless of the value of parameters s and µ. We conducted the whole analysis by preferring the variable y 2 over y 1 and z because it is the only practical biologically measurable quantity. Our results are especially convenient since they only depend on this quantity.

Application of the stability analysis

In this section we will apply the previous stability analysis with two different functions modeling the immune system.

The original model

As in [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF], we set f (x) = d -α x 1+εx 2 . Its profile is represented in Figure 2. In this case, the derivative f (x) = -α 1-εx 2 (1+εx 2 ) 2 , vanishes only for x = ε -1/2 . We assume that α 2 ≥ 4εd 2 , and denote:

y min = α - √ α 2 -4εd 2 2εd , y max = α + √ α 2 -4εd 2 2εd .
As soon as M > y max , we can apply Corollary 1. The function P θ is then the third-order polynomial,

P θ = (X -y min )(X -y max )(X -M ) + θ εd X + M d 2 r (1 + εX 2 ).
With this particular example, we can give an upper bound for A B on R + :

A B (X ≥ 0) = d (X 2 -(y min + y max )X + 1 ε )(M -X) (X + d2 r M )(X 2 + 1 ε ) ≤ d M -X X + d2 r M ≤ rd d 2 .
Hence, θ max is bounded by rd d2 . This value corresponds to the transcritical bifurcation between the low disease and the disease-free steady states. This means that the saddle-node bifurcation occurs for lower values of θ than the transcritical bifurcation.

The bifurcation diagram and stability chart given by Corollary 1 are reproduced in Figures 3 and4. Fig. 3: Bifurcation diagram depending on parameter s variation (immune cells source). Solid red curves are for asymptotically stable steady states with real dominant eigenvalue, dashed red curves are for asymptotically stable steady states with complex dominant eigenvalue, and dotted curves are for unstable steady states. On the left graph the disease-free steady state corresponds to ȳ2 = 0, so it does not appear with the log-scale. For s big enough, the two highest steady states disappear, then low and trivial steady states cross each other.

We applied the simplified model ( 4) to the patient data from [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF] in order to compare the relevance of two models in Figure 5.

For many patients, the fits in Figure 5 are comparable between the complete model ( 1) and the simplified model (4). Intriguingly, for some patients (such as Patients 1, 3 and 5), the fits seem more satisfactory with the simplified model. Typically when the initial decay is slow, it looks as if the oscillations are better controlled by the simplified model ( 4). It might be due to the fact that in order to be consistent with [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF], we decided to assume a priori that y min and y max exist, which then constrains the parameters. Compared with the full model, (1), the model (4) cannot capture rapid oscillations and sharp decays.

A mechanistic model

We apply here our procedure to a second example, a mechanistic model for describing the dynamics of CML [START_REF] Clapp | Applying Mathematical Models to Study the Role of the Immune System in Chronic Myelogenous Leukemia[END_REF]. In this model, there are no terms representing immunosuppression of the immune-leukemia interactions. Instead, mature Fig. 5: Fits of the original (dashed line) and reduced (solid line) models to six patients treated with first-line Imatinib. The base-10 log of the BCR-ABL ratio is plotted against time, in months. Dots represent patient data. Dotted lines approximate the minimum leukemic level that is detectable by RT-PCR. Dots along this line represent zero measurements, meaning CML cells were not detected. The parameters used for the fits to the original and simplified model are respectively in Tables 2 and3, universal parameters are in Table 1. For the immune dynamics, f is given by f (x) = d + νx -α x 1+εx . Its profile is represented in Figure 6. In this case, the derivative f (x) = ν -α(1 + εx) -2 , vanishes only for x = 1 ε ( α ν -1). We assume (α -εd -ν) 2 -4ενd > 0, and denote:

y min = (α -εd -ν) -(α -εd -ν) 2 -4ενd 2εν , y max = (α -εd -ν) + (α -εd -ν) 2 -4ενd 2εν .
As soon as M > y max , we can apply Corollary 1. The bifurcation diagram and stability chart given by Corollary 1 are reproduced in Figures 7 and8.

Discussion and conclusion

In this manuscript we presented a simplified, analytically tractable version of a model for tumor-immune interaction recently proposed [START_REF] Clapp | Implication of the autologous immune system in BCR-ABL transcript variations in chronic myelogenous leukemia patients treated with imatinib[END_REF]. This simplification allows us to conduct a complete analytic study of the model. We showed the existence of a disease-free steady state and established conditions for the existence of disease steady states. We characterized the stability of each of these steady states, and identified the transitions in the number of steady states and in their stability. Fig. 7: Bifurcation diagram depending on parameter s variation (immune cells source). Solid red curves are for asymptotically stable steady states with real dominant eigenvalue, dashed red curves are for asymptotically stable steady states with complex dominant eigenvalue, and dotted curves are for unstable steady states. On the left graph the disease-free steady state corresponds to ȳ2 = 0, so it does not appear with the log-scale. For s big enough, low and trivial steady states cross each other, then the two highest steady states disappear. Bistable: two disease s.s.

Monostable: one disease s.s.

Monostable: one non-disease s.s.

Bistable: one disease-free s.s., one disease s.s. £ (immune system e cacy) Fig. 8: Stability chart in the plan (θ, M ), for Subsection 5.2, where θ = µs d2 . We still have the possibility of having two stable steady states (one with y 2 < y min one with y 2 > y max ). One of the main difference is that a stable disease free steady state may coexist with a stable disease steady state. In the first model, this is mainly due to the fact that the immunosuppression is not saturated and high levels of y 2 leads to very small values of z, below s f (0) .

It was anticipated that the simplified model would possess a disease-free steady state. This state corresponds to the asymptotic solution when treatment is effective. Whether the leukemic cells completely vanishes under treatment is a biological open question [START_REF] Chu | Persistence of leukemia stem cells in chronic myelogenous leukemia patients in prolonged remission with imatinib treatment[END_REF][START_REF] Ross | Patients with chronic myeloid leukemia who maintain a complete molecular response after stopping imatinib treatment have evidence of persistent leukemia by DNA PCR[END_REF]. The model we analyze leads to a low disease steady state (with leukemic load below the immune window) that corresponds to the asymptotic solution when treatment is not sufficiently efficient to eradicate the disease, but is sufficient to allow the immune system to keep it under control at a low level.

We obtained complete generic one-parameter bifurcation diagrams (Figures 3 and7), and showed that a low-disease steady state exists for a wide range of parameters. When such a state exists, it is stable. The low-disease steady state becomes non-positive and disappear from the biologically relevant phase-space at a transcritical bifurcation with an unstable disease-free steady state, at which they exchange their stability. Large disease steady state can coexist; they are generically created and destroyed through saddle-node bifurcations. No steady state can ever be destabilized through a Hopf bifurcation. Although we have not shown that no limit cycles can occur, we could not find any non-local bifurcations leading to limit cycles, such homoclinic bifurcations, or saddle-node bifurcation on an invariant cycle (SNIC). The low-disease steady state is the only steady state that can be a focus, and it is always a stable focus.

Two qualitatively different bifurcation diagrams were identified: (i) a stable disease-free steady state and existence of a large disease steady state and mutually exclusive; and (ii) co-existence of a stable disease-free steady state and a large disease steady state possible.

The asymptotic stability we highlighted is a local stability of steady states, which does not guarantee the convergence of any solution to the steady state. We have not yet derived any theoretical result on the basin of attraction of asymptotically stable steady states. Numerically, the low steady state basin of attraction seems much smaller than the basin of attraction of the high steady state. That could be interpreted by the scarcity of solutions converging to the low disease steady state. For example, this model allows a state with low amount of leukemic cells and high amount of immune cells to converge to the highest disease steady state with high amount of leukemic cells and low amount of immune cells.

The stability results obtained in this paper could help in interpreting the clinical state of remission after treatment cessation. Indeed, remission could correspond to attraction by the low disease steady state, while relapse could correspond to attraction by a higher stable disease steady state. It is useful to see the weight of each parameter in the dynamics of solutions. For example, Figures 4 and8 highlight the relevance controlling the autologous immune system efficacy in order to constrain solutions to converge to a low disease steady state or to a disease-free steady state.

Fitting either the original or the simplified model to individual patients yields excellent results in most cases (Figure 5). The simplified model has a tendency to predict oscillatory levels of BCR/ABL ratio, due to the focus nature of the low disease steady state. Whether the fluctuations in clinical BCR/ABL ratio correspond to deterministic oscillations or to stochastic fluctuations is not clear, but our steady state analysis suggests that oscillations are a signature of the immune control of the CML.

The bifurcation diagrams suggest three possible scenarios for the outcome of a treatment cessation. In any of these scenarios, relapse is possible. First, if there is co-existence of a stable large disease steady state and a stable low steady state, no cure is possible and the success of treatment cessation is determined by how attractive the low-disease steady state is. Second, if there is co-existence of a large disease steady state and a stable disease-free steady state, it is possible that treatment cessation (or continuation) could actually lead to eradication of the tumor. Third, no large steady state exists, and the low-disease steady state is the clinical disease state. In this case, the disease-free steady state is always unstable, and treatment cessation is bound to fail.

Our stability analysis is adaptable to a large class of functions modeling the immune system degradation. The double-monotony of this function is a biologically acceptable condition, corresponding to immune activation and immunosuppression.

Fig. 2 :

 2 Fig. 2: Functions f and P profiles. The part where P is positive is represented in red.

×10 6 MFig. 4 :

 64 Fig. 4: Stability chart in the plan (θ, M ), for Subsection 5.1.

Fig. 6 :

 6 Fig.6: Functions f and P profiles. The part where P is positive is represented in red.

Table 1 :

 1 Universal parameters values used in the reduced model.

	Parameter Value
	r	0.007775
	K	41.667
	a1	1.350e5
	d2	0.0375
	s	z(0) * d
	y1(0)	K
	y2(0)	1.5e8
	z(0)	120

Table 2 :

 2 Estimated parameters with the original model[START_REF] Burchert | Frequency of CTLA-4 Receptor Ligand (CD86, B7.2)-Positive Plasmacytoid Dendritic Cells Predicts Risk of Disease Recurrence after Tyrosine-Kinase Inhibitor Discontinuation in Chronic Myeloid Leukemia: Results from a Prospective Substudy of the Euroski Trial[END_REF]

	Patient inh1	inh2	d	µ	ymin	ymax
	1 4.612 92.3215 0.031 9.964e-7 4.994e4 5.598e5
	2 1.456 545.150 0.099 1.504e-8 3.765e4 2.759e5
	3 591.591 14.568 0.040 2.371e-7 3.132e3 2.228e4
	4 486.315 226.000 0.075 2.879e-8 3.536e2 1.684e3
	5 50.988 79.645 0.005 1.271e-6 1.182e3 5.482e4
	6 30.208 359.979 0.371 2.263e-7 4.959e3 1.353e4

Table 3 :

 3 Estimated parameters with the simplified model (4) Menten term, αy 2 z/(1 + εy 2 ), that represents the stimulation of immune cells by leukemic cells, with maximum rate αε -1 and Michaelis constant ε -1 ; and (ii) a mass action term, νy 2 z, representing the killing of immune cells by leukemic cells, at a rate ν.

	Patient k inh	d	µ	ymin	ymax
	1 2.521e2 0.051 3.647e-6 6.610e4 3.624e5
	2 1.133e3 0.026 2.405e-8 3.831e4 3.055e5
	3 4.205e2 0.054 4.224e-7 1.617e4 3.133e5
	4 5.691e3 0.181 8.499e-6 1.206e3 1.090e4
	5 4.594e3 0.038 5.723e-9 1.841e3 3.401e4
	6 2.853e3 0.058 1.358e-9 7.143e3 7.576e4

Acknowledgements The work of GC was supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE1322106. The work of DL was supported in part by the John Simon Guggenheim Memorial Foundation. The work was supported by the Inria Partnerships Program grant "Modelling Leukemia".

The parameter θ will be used as the bifurcation parameter. It has no biological meaning, but it allows the polynomial B to be unitary, and thus simplifies the calculations. Otherwise it could be feasible to use parameters µ or s that represent the immune system aggressiveness and the immune cells source, as we do in the Section 5. The parameter θ may vary in [0; rd/d 2 ] to respect condition (C 1 ). The limit value θ = 0 is allowed to be reached when s = 0, but not when µ = 0. Indeed, in this last case the link between steady state and roots of polynomial P does not hold anymore.

The following lemma highlights the link between the stability of a steady state x and the sign of the jacobian matrix at x, under (C 2 ). Its proof is given in the Appendix.

Lemma 1 Let x = (ȳ 1 , ȳ2 , z) be a steady state of (4) such that f (ȳ 2 ) > 0. The steady state x is asymptotically stable if and only if, det(J(x)) < 0.

Lemma 1 affirms that the key to stability is the Jacobian matrix determinant. The following two lemmas highlight the link between this determinant and the polynomials A and B. Their proofs are given in the Appendix.

Lemma 2 Let x = (ȳ 1 , ȳ2 , z) be a steady state of (4). The determinant of the corresponding Jacobian matrix is:

Lemma 3 Let x ≥ 0. The following conditions are equivalent:

x is a double root of P θ , where θ := A(x) B(x) .

We first note that the unique zero of f must lie between the zeros y min and y max of f . This implies that f is decreasing in [0; y min ] and increasing in [y max ; +∞). Therefore, f stays positive on [y max ; M ], steady states x(1) and x(2) verify the condition (C 2 ).

Hence, Lemma 1 ensures that the stability of the steady state x = (ȳ 1 , ȳ2 , z) is given by the sign of det(J(x)), which is equal to the sign of

The steady states are zeros of P θ , so they satisfy A(x) B(x) = θ. Since A and B do not depend on θ, there exists a maximal value for θ such that this equation admits a solution in (y max ; M ), θ max = sup x∈(ymax,M ) { A B (x)}. This bound stands for a real x 0 ∈ (y max ; M ) because A B (y max ) = A B (M ) = 0. The function A B -AB vanishes in x 0 because it is an extremum for A B .Then x 0 is a double zero of P θmax .

We set θ < θ max . A double zero x 1 for P θ satisfies A B (x 1 ) = θ, and (A + θB )(x 1 ) = 0. Then (A B -AB )(x 1 ) = 0, so x 1 is a critical point of A B . By our hypothesis A B admits a unique critical point, so x 1 = x 0 . This is a contradiction since θ < θ max . This means that P θ admits only simple zeros in (y max , M ). We have

Since P θ (y max ) > 0 and P θ (M ) > 0, P θ admits exactly two zero in (y max , M ), distributed on both sides of x 0 . A third zeros would induce to a second critical point for A B , which is forbidden. The function A B changes its direction only at its critical point x 0 . It is easy to verify that (A B -AB )(M ) < 0 and (A B -AB )(y max ) > 0, hence A B is strictly increasing on [y max ; x 0 ] and strictly decreasing on [x 0 ; M ]. We can now explore the three cases of the theorem:

1. If θ < θ max , there exist two distinct steady states x(1) and x(2) such that ȳ(1)

2 . In addition x(1) is unstable and x(2) is asymptotically stable. 2. If θ > θ max , P θ cannot vanish on [y max , M ], so only the disease-free and maybe the low disease steady state remains. We have seen in Proposition 5 that, as soon as the disease-free steady state is unstable, the low one is asymptotically stable. 3. If θ = θ max , then there exists a unique steady state in [y max , M ], resulting from the fusion of two steady states with different stabilities: it is a saddlenode bifurcation. Let x = (y 1 , y 2 , z) be a perturbation of the steady state x such that y 2 stays in [y max , M ]. We have

Hence z is always increasing near the steady state which, therefore, cannot be stable.

The next theorem is a stronger form of Theorem 1.

Theorem 2. Denote by

For all θ in [0; θ max )\{θ 1 , . . . , θ n }, there exist an positive even number of steady states with ȳ2 in (y max ; M ). Their stability is alternated, the largest in term of value of ȳ2 is asymptotically stable. 2. For all θ > θ max there is no steady state with ȳ2 in (y max ; M ). 3. Generically, there is a saddle-node bifurcation for each θ ∈ {θ 1 , . . . , θ n }.

Corresponding double steady states are generically unstable.

Proof. Using arguments identical to those used in proving Theorem 1, we conclude that x 0 is a double zero of P θ . We set θ ∈ [0; θ max )\{θ 1 , . . . , θ n }. The

Since P θ (y max ) > 0 and P θ (M ) > 0, P θ admits a positive even number of zeros in (y max ; M ). The function A B stays monotonic between two consecutive extrema. In addition, the sign of its derivative alternates between extrema. Calculations give A B (M ) = 0 and (A B -AB )(M ) < 0, so A B is locally positive and decreasing on the left side of M .

We fix θ < θ max , and denote

B must be locally decreasing near x θ . We can now explore the three cases of the theorem:

1. If θ ∈ [0, θ max )\{θ 1 , . . . , θ n }, the steady states set is

ȳ1 and z verify 5 and 6 .

Since A B is decreasing in max{ ȳ2 / x ∈ S θ ), that steady state is asymptotically stable by Lemmas 2 and 1. The stability of the other steady states alternates in terms of ȳ2 value. Since #S θ is even, the lower steady state is unstable. 2. If θ > θ max , P θ cannot vanish on [y max , M ], so only the disease-free and maybe the low disease steady state remains. We have seen in Proposition 5 that as soon as the disease-free steady state is unstable, the low one is asymptotically stable.

3. If θ = θ max then there exists a unique steady state, resulting from the fusion of two steady states with different stabilities: it is a saddle-node bifurcation. 4. If θ ∈ {θ 1 , . . . , θ n }, we let x = (y 1 , y 2 , z) be a perturbation of the steady state x such that y 2 stays in [

In both cases, z is monotonic next to the steady state which, therefore, cannot be stable.

Corollary 1 Let f be a rational fraction of the form f = N D , where D is a positive polynomial with degree two or less, and N is a quadratic polynomial that admits two roots y min and y max in [0, M ], such that f vanishes only once in [0, M ]. Then there exists a value θ max such that:

1. For any θ in [0, θ max ) there exist two disease steady states with ȳ2 in (y max , M ).

The lower one is unstable and the higher one is asymptotically stable. 2. For any θ in (θ max , +∞) there is no steady state with ȳ2 in (y max , M ). 3. There is an unstable saddle-node bifurcation at θ = θ max .

Proof. The function f admits two zeros in (0, M ), which we denote as y min and y max . We write f as f (x) = (x-ymin)(x-ymax)

, and adapt Theorem 2: The first steady state is (0, 0, s d ). It is the disease-free steady state. Other steady states are given by:

, and z = sD(ȳ 2 ) (ȳ 2 -y min )(ȳ 2 -y max )

,

where ȳ2 is a positive zero of the cubic polynomial:

such that ȳ2 / ∈ [y min , y max ]. We keep the previous notations: A(X) = (X -y min )(X -y max )(X -M ) and B(X) = (M d2 r + X)D(X). To prove the corollary, we only have to verify that A B has only one critical point in (y max , M ). We assume that A B admits a second critical point x 1 in (y max , M ), and denote θ 1 = A B (x 1 ). Then x 1 is a double root for P θ1 . For all θ lower than θ max , P θ admits one root in [0, y min ] and two distinct roots in [y max , M ]. Since it is a cubic polynomial, there are no more roots. Hence P θ cannot have any double root in [y max , M ].

This means that θ 1 = θ max . Since P θmax is a cubic polynomial that admits x 0 and x 1 as double roots, necessarily x 1 = x 0 . We showed that A B admits one and only one critical point. Applying Theorem 1 concludes the proof.

We have seen in this section that the bifurcation diagram is uniquely determined by transposing the graph of the function A B . This means that no other

APPENDIX

The following equations are satisfied by positive steady states. They will be used in technical lemmas.

.

Lemma 1 Let x = (ȳ 1 , ȳ2 , z) be a steady state of ( 4) such that f (ȳ 2 ) > 0. The steady state x is asymptotically stable if and only if, det(J(x)) < 0.

Proof. The determinant of J(x) is the product of all its eigenvalues, so it is equal to -χ J (0). We note that the polynomial χ J is convex on R + , and that

Here we used f (ȳ 2 ) > 0 given by Proposition 2. We consider two cases:

Case 1: If χ J (0) < 0, then there exists a root in R * + . In this case the steady state x is unstable.

Case 2: If χ J (0) > 0, then χ J (0) > 0 and, by convexity, χ J stays non-negative on R * + . If χ J admits three real roots, they are all negative, and the steady state is asymptotically stable. Otherwise, there exists two conjugate complex roots. In this case, we denote by x the negative real root, and by z one of the two complex roots. First, note that since f (ȳ 2 ) > 0, we have:

As χ J has only one sign change on R, at x, we deduce that -r K ȳ1 < x. Yet,

< 0.

We showed that complex roots have negative real parts in the case where f (ȳ 2 ) is positive. Therefore, the steady state is asymptotically stable.

Lemma 2 Let x = (ȳ 1 , ȳ2 , z) be a steady state of (4). The determinant of the corresponding Jacobian matrix det(J(x)) has a sign opposite to that of A B -AB .

Proof. We compute

Let x = (ȳ 1 , ȳ2 , z) be a steady state for the system (4). The determinant of the Jacobian matrix at x is det

Hence sign(det(J(x))) = -sign((A B -AB )(ȳ 2 )).

Lemma 3 Let x ≥ 0. The following conditions are equivalent:

x is a double root of P θ , where θ := A(x) B(x) .

Proof. Let x be a positive zero of A B -AB . We fix θ := A(x) B(x) . Then,

Hence x is a double zero of P θ . Reciprocally, let (θ, x) be such that x is a double zero of P θ . Since P θ (x) = 0, θ = A B (x). Also, since P θ (x) = 0, (A B -AB )(x) = 0, which means that x is a zero of A B -AB . Proof. First, positivity of all coefficients ensures that the real roots of P cannot be positive. Second, as P (0) = c > 0, P necessarily admits a real, negative root. It remains to characterize the two other roots. If real and negative, then the dominating root is negative and the lemma is proven. Therefore we may assume that the two remaining roots are complex, and we need to determine the sign of their real part.

Let λ be the first negative root of P . We can factor P = X -λ X 2 + (a + λ)X + (b + aλ + λ 2 ) .

By our assumption, P admits a pair of complex roots. We want to compare In order to find the sign of Re(z), we need to compare a and λ. We have P (-a) = -ab + c. Since we assume that, ab > c, we obtain P (-a) < 0. As P has only one real root λ, we obtain -a < λ. This leads to Re(z) < 0, and the dominating root of P has its real part negative.