A class of constrained optimal control problems arising in an immunotherapy cancer remission process

1.   Introduction

Myeloproliferative Neoplasms (MPNs) is a group of hematopoietic stem cell disorders, including essential thrombocytosis (ET), polycythemia vera (PV) and primary myelofibrosis (PMF) ^[1,2]. The pathogenesis of these neoplasms is yet to be fully discovered. For patients with MPNs, the mutation $JAK2$V617F is found present in the most cases of ET (50%) and in 95% of the cases with PV and PMF ultimately leading to acute myeloid leukemia (AML) ^[3,4]. This suggests a biological continuum where the diseases evolve from early cancers (ET and PV) into the advanced myelofibrosis stage, with an increasing load of $JAK2$V617F mutations from a low burden at ET and PV to a high load ^[2,5]. MPNs imply an increased risk for the development of other cancers ^[1,4].

Recent research supports that MPNs can be regarded as chronic inflammatory diseases and MPNs has been described as a "human inflammation model", which leads to premature atherosclerosis, clonal evolution and an increased risk of second cancers. ^[2,3,6]. This is based on evidence from clinical observations, experiments and molecular studies ^[3].

Several insightful theoretical studies have been published on control dynamics of biological networks. Mathematical models have been proposed ^[7,8,9] describing the control networks for regulation of stem cell lineage.Mathematical modelling of cancer is useful for understanding of cancer initiation, progression ^[10,11], to confirm or dismiss biological/medical hypotheses, and to study effects of single or multi modality treatments in silico. The mathematical model presented in ^[12] shows that successful therapy may eliminate tumour stem cells. A five-dimensional model given in ^[13] includes active and quiescent stem cells, progenitor cells, mature cells and one immune compartment describing chronic myelogenous leukemia.In ^[14] a mathematical model of cancer stem cell dynamics is proposed and the different scenarios of cancer initiation and possible treatments strategies have been discussed. The mathematical model given in ^[15] is useful for investigating the impact of cytokine-dependence of acute myeloid leukemic cells. In addition, the model allows distinguishing between cytokine-dependent and cytokine-independent acute myeloid leukemia (AML) and both scenarios are supported by patient data.

However, only a few mathematical models of MPNs exist. Some work includes investigation of the origin of myeloid malignancies with MPNs as a particular example ^[16]. In ^[17], a two dimensional model of MPNs is investigated without including the immune response dynamically. The Cancitis model including chronic inflammation as the trigger and driver of MPNs was proposed in ^[5]. In ^[5], T-cells are not explicitly considered whereas, in ^[18], the effect of these cells has been included specifically. The analysis of a two dimensional mathematical model ^[18] is used to discuss in silico effect of existing and novel treatments. The model presented here is identical to the model presented in ^[5] except for the simpler functional form of the stem cells niche interaction used here and in ^[18].

In the present paper we conduct a thorough mathematical investigation of the Cancitis model and explore the intricate coupling between inflammation and MPNs. We address the following questions which have not been systematically investigated previously: Which steady states of the system are feasible and which trajectories are attracted to the steady states? How do the number and stability of the steady states change when varying the parameters, in particular, the exogenous inflammatory stimuli, self-renewal and death rates of stem cells, and inhibitory strength of the stem cell niche interaction? Which set of clustered parameters control the dynamics of the system? Does the analysis suggest correlated parameters? The bio-medical applications of the model analysis are discussed, e.g. how the inflammation influences the transition between healthy and diseased states. In addition, the analysis predicts effects of ongoing and potential combination therapies.

2.   The Cancitis model

The Cancitis model stated in ^[5] is illustrated in Figure 1, with the system of differential equations shown in system (2.1). In this section the details of the model and the reasoning behind it is presented.

The model describes the proliferation of hematopoietic stem cells (HSC) into hematopoietic mature cells (HMC) and likewise malignant stem cells (MSC) into malignant mature cells (MMC). Additionally, the model considers the number of dead cells and the level of inflammation. The debris from the dead cells stimulates the immune response, which in turn affects the self-renewal rate of both HSCs and MSCs.

The model consists of six ordinary differential equations one for each compartment; the number of HSC ($x_0$), the number of HMC ($x_1$), the number of MSC ($y_0$), the number of MMC ($y_1$), the number of dead cells ($a$), and the level of inflammation ($s$).

The equations are of the general form,

and read specifically,

with

The expressions for the inhibitory niche feedback are chosen as Michaelis-Menten-like functions in contrast to ^[5],

A stem cell can proliferate in three ways; symmetric self-renewal (resulting in two new stem cells), asymmetric self-renewal (resulting in one stem cell and one progenitor cell) and symmetric differentiation (resulting in two progenitor cells). The rate of self-renewal is denoted as $r_x$ and $r_y$ for HSC and MSC respectively. The self-renewal of stem cells is known to be inhibited by self-regulating niche feedback ^[19], resulting in decreased self-renewal when the level of stem cells in the bone marrow is high. Adopting the approach taken in ^[12], ^[20] and ^[21], this is implemented by Michaelis-Menten-like functions $\phi_x(x_0, y_0)$ and $\phi_y(x_0, y_0)$, shown in Eq (2.3b). Allowing the feedback to be different for HSC and MSC, the constants $c_{xx}$ and $c_{xy}$ capture the effects of HSC and MSC on the self-renewal of HSC, while $c_{yx}$ and $c_{yy}$ capture the corresponding effects on the self-renewal of MSC. Additionally, the inflammatory level also affects the self-renewal ^[22,23]. This leads a to self-renewal term per cell of $r_x \phi_x s$ and $r_y \phi_y s$ for HSC and MSC respectively. The parameter $c_{ij}$ describes the inhibitory strength of the signalling bone marrow niche feedback from cell type $j$ onto cell type $i$. It is generally assumed that $c_{yy}\leq c_{yx} \leq c_{xy} \leq c_{xx}$, since malignant cells are less sensitive to inhibitive niche feedback than hematopoietic cells ^[22,24].

In ^[25], a multi compartmental model is proposed relying on a single external feedback mechanism. It is shown that the equilibrium level of mature cells depends only on the self-renewal parameters for the HSC and it is independent of the other compartments. Therefore, the progenitor cells are considered as intermediate steps between stem cells and mature cells, and are implicitly accounted for by multiplication factors $A_x$ and $A_y$ for HSC and MSC respectively. The rate at which the HSC reduces in transforming to HMC is denoted by $a_x$ while the similar rate for MSC transforming to MMC is denoted by $a_y$. As such, the HMC and MMC accordingly increase with rates $a_x A_x$ and $a_y A_y$ respectively. To account for the cell apoptosis, the four types of cells are removed with rates $d_{x_0}$, $d_{x_1}$, $d_{y_0}$ and $d_{y_1}$, for the corresponding cell types.

Genetic mutations are by nature to be described as Poisson processes ^{[26,27,28,29]}. However, not all mutations are malignant; only mutation which happens on a particular location of the DNA, i.e. at specific amino acids causes a specific mutation, e.g. the JAK2V617F mutation. The probability for hitting a specific location is about 1/30000. In ^[30] the average mutation probability is estimated to 0.0035 per year, which corresponds to a specific mutation probability of $0.0035/30000 = 1.2 10^{-7}$ per year. Thus, the probability for one specific malignant mutation is about $10^{-7}$ per cell per year. Moreover, the mutation is affected by the inflammation, $s$ ^[31,32], which is explicitly stated, and resulting in the effective mutation rate $r_m s$. Assuming three sequential mutations are needed to generate a specific malignant stem cell the resulting probability becomes much higher ($10^{-25}$ per year per cell if the mutations are assumed independent). This could be implemented in the otherwise deterministic model but it would increase the computational cost, since it depends on both the probability of a single cell mutation and the number of potential mutating cells at a given time, which itself is determined by the preceding mutational history. To avoid such complications we initialize by having a single MSC and none MMC, and put the mutation rate to zero. This is justified by the fact that the probability of a single cell mutating is small compared to the self-renewal of the MSCs. Thus, the first mutation drives the development leaving a later identical mutation insignificant to the dynamics, which is confirmed by numerical simulations.

The number of dead cells has an up-regulatory effect on the immune response denoted $r_s$. External environmental factors also influence the inflammatory level. This is captured in the model by the term $I$. Throughout we take $I > 0$, as a perfect sterile environment is an utopic idealization. This term may vary over time due to environmental changes, but in our mathematical analysis we will consider $I$ as piecewise constant. The inflammation, $s$, is down-regulated naturally by the eliminating rate $e_s$.

Additionally, the change in the amount of dead cells per time is given by the death rate times the number of cells minus the clearance by the immune system. As given in ^[33] clearance is described by a second order equation $-e_a a s$ since the engulfed immune cells have to meet the dead cells debris to mediate endocytosis. Thus, clearance is bilinear in both $a$ and $s$ representing the activity of the immune system, eliminating the dead cells with an elimination rate $e_a$.

Initial conditions for the Cancitis model in equations are needed for the given system of differential equations (2.1–2.3b). Here, we mainly focus on the model after the first mutation, i.e. with $y_0(0) = 1$, $y_1(0) = 0$, $r_m = 0$, and the remaining variables as those in the healthy steady state (see below). All other parameter values are assumed to be positive.

2.1.   Steady states of the model

The stable steady states are attractors in the six dimensional phase space. This motivates systematic study of the existence and location of steady states and how this is affected by perturbing parameter values.

Motivated by the biology where the number of cells and concentrations are required to be non-negative numbers, we will use the terminology that a steady state is admissible if and only if all the components are non-negative, i.e. if and only if a steady state is in the non-negative octahedron.

Consider the system of Eqs (2.1–2.3b). For steady state solutions, $\dot{x}_0 = \dot{x}_1 = \dot{y}_0 = \dot{y}_1 = \dot{a} = \dot{s} = 0$,

The values of $x_0$, $x_1$, $y_0$, $y_1$, $a$ and $s$ at steady state are denoted as $\bar{x}_0$, $\bar{x}_1$, $\bar{y}_0$, $\bar{y}_1$, $\bar{a}$ and $\bar{s}$ respectively. The admissible steady states can be classified as,

● a trivial steady state if and only if $\bar{x}_0 = \bar{y}_0 = 0$,

● a (purely) hematopoietic steady state if and only if $\bar{y}_0 = 0$ and $\bar{x}_0 > 0$,

● a (purely) malignant steady state if and only if $\bar{x}_0 = 0$ and $\bar{y}_0 > 0$, or

● a co-existing steady state if and only if $\bar{x}_0 > 0$ and $\bar{y}_0 > 0$.

The admissibility of steady states necessitates certain inequalities to be fulfilled, which leads to restriction on the parameters, e.g. from Eq (2.4f), non-negativity of $\bar{a}$ requires $\bar{s} \geq \frac{I}{e_s}$. $I$ is assumed to be positive thus $\frac{I}{e_s} > 0$. Hence

for any admissible steady state. Note that, $\bar{x}_1$ and $\bar{y}_1$ are non-negative if and only if $\bar{x}_0$ and $\bar{y}_0$ are non-negative, respectively. Substituting expressions of $\bar{x}_1$ and $\bar{y}_1$ in Eq (2.4e) gives,

where $\beta_x = a_x A_x+d_{x0}$ and $\beta_y = a_y A_y+ d_{y0}$. Thus, Eqs (2.4e–2.4f) and (2.6) result in a second order equation in $\bar{s}$ having the general solution,

where $\zeta = \frac{4 r_s e_s}{e_a I^2} > 0$. $\bar{s}_-$ is negative for positive $\bar{x}_0$ or $\bar{y}_0$. In case, $(\bar{x}_0, \bar{y}_0) = (0, 0)$, Eq (2.7) leads to $\bar{s}_{-} = 0$ contradicting $\bar{s} > 0$. Thus, $s = \bar{s}_-$ is not biologically realizable and we therefore put $\bar{s} = \bar{s}_+$ in the further analysis. Note that, for non-trivial steady states, $\bar{s} > \frac{I}{e_s}$. In addition, non-negativity of $\bar{x}_0$ and $\bar{y}_0$ implies non-negativity of $\bar{x}_1$, $\bar{y}_1$, $\bar{s}$ and $\bar{a}$.

Hence, the existence of admissible steady states $\bar{x}_0$, $\bar{y}_0$, $\bar{x}_1$, $\bar{y}_1$, $\bar{s}$ and $\bar{a}$ follows from Eqs (2.4a), (2.4b) and (2.7).

Below we make a complete analysis of the existence of various steady states depending on how $I$ relates to the remaining parameters. This choice is due to the fact that the external inflammatory stimuli $I$ is of great interest in health care and to elucidate consequences of using inflammation inhibitors as part of treatment.

Proposition 1. A trivial steady state $E_{0}$ always exists,

Proof. Follow directly from Eqs (2.4e) and (2.5).

Hematopoietic steady states may exist depending on the rest of the parameter values. As above we chose the inflammatory stimuli $I$ as the leading parameter and make a complete analysis of possible hematopoietic steady states. The analysis of the existence of the hematopoietic steady states depends crucially on the following lumped parameters,

the last two always being positive.

Proposition 2. Two hematopoietic steady states $E_{H\pm}$ may exist in the following cases,

● If $\alpha_x^2 \leq \frac{\zeta_{H3}}{4}$ then $E_{H+}$ exists if and only if $I > 0$.

● If $\frac{\zeta_{H3}}{4} < \alpha_x^2 \leq \zeta_{H3}$ then $E_{H+}$ exists if and only if $I \geq I_{H}$.

● If $\zeta_{H3} < \alpha_x^2$ then $E_{H+}$ exists if and only if $I > \zeta_{H2}$.

● If $\alpha_x^2 \leq \frac{\zeta_{H3}}{4}$ then $E_{H-}$ exists if and only if $I \leq \zeta_{H2}$.

● If $\frac{\zeta_{H3}}{4} < \alpha_x^2 \leq \zeta_{H3}$ then $E_{H-}$ exists if and only if $I_{H} \leq I \leq \zeta_{H2}$.

● If $\zeta_{H3} < \alpha_x^2$ then $E_{H-}$ does not exist.

In case of existence, $E_{H+} = (\bar{x}_{0H+}, \bar{x}_{1H+}, 0, 0, \bar{a}_{H+}, \bar{s}_{H+})$ is given by $\bar{x}_{0H+} = \frac{1}{2 e_s c_{xx}\alpha_x} \left(I-\zeta_{H1} + \sqrt{(\zeta_{H1}-I)^2-4 e_s \alpha_x(\zeta_{H2}-I)}\right)$, $\bar{s}_{H+} = \alpha_x (1+c_{xx} \bar{x}_{0H+})$, $\bar{a}_{H+} = \frac{\beta_x \bar{x}_{0H+}}{e_a \bar{s}_{H+}}$, and $\bar{x}_{1H+} = \frac{a_x A_x \bar{x}_{0H+}}{d_{x1}}$ whereas $E_{H-} = (\bar{x}_{0H-}, \bar{x}_{0H-}, 0, 0, \bar{a}_{H-}, \bar{s}_{H-})$ is given by $\bar{x}_{0H-} = \frac{1}{2 e_s c_{xx}\alpha_x} \left(I-\zeta_{H1} - \sqrt{(\zeta_{H1}-I)^2-4 e_s \alpha_x(\zeta_{H2}-I)}\right)$, $\bar{s}_{H-} = \alpha_x (1+c_{xx} \bar{x}_{0H-})$, $\bar{a}_{H-} = \frac{\beta_x \bar{x}_{0H-}}{e_a \bar{s}_{H-}}$, and $\bar{x}_{1H-} = \frac{a_x A_x \bar{x}_{0H-}}{d_{x1}}$.

Proof. A hematopoietic steady state $E_{H}$ follows from Eqs (2.4a) and (2.7) with $y_0 = y_1 = 0$ as possible positive solutions to

For the solutions to (2.13) to be real,

In case $I\geq\zeta_{H2}$, (2.14) is always fulfilled. In case $I < \zeta_{H2}$, (2.14) is equivalent to

Solving for $I$ we get,

with

From Eqs (2.14) and (2.16) it follows that the solutions to Eq (2.13) are real for $I\geq\zeta_{H2}$ or $I_H\leq I < \zeta_{H2}$ in case $I_H < \zeta_{H2}$.

Whenever the solutions to (2.13) are real, they are given by

which depends on the sign of the following five quantities,

where the last one immediately implies that the criteria for real solutions of Eq (2.13) is $I\geq I_H$. If the solutions, $\bar{x}_{0H\pm}$, are positive and real, then the formulas for the remaining variables easily follow from Eqs (2.4a–2.4f).

To continue we first consider $\bar{x}_{0H+}$ and afterwards $\bar{x}_{0H-}$.

For $\alpha_x^2 \leq \frac{\zeta_{H3}}{4}$, it follows from Eqs (2.19)–(2.22) that $\zeta_{H1} < 0$, $\zeta_{H1} < \zeta_{H2}$, $I_H\leq 0$, and $\zeta_{H1} < I_H$. Thus, $\bar{x}_{0H+} > 0$ if and only if $I > 0$.

For $\frac{\zeta_{H3}}{4} < \alpha_x^2 \leq \frac{\zeta_{H3}}{2}$, it follows from Eqs (2.19)–(2.22) that $\zeta_{H1}\leq 0$, $\zeta_{H1} < \zeta_{H2}$, $I_H > 0$, and $\zeta_{H1} < I_H$. Thus, $\bar{x}_{0H+} > 0$ if and only if $I > I_H$.

For $\frac{\zeta_{H3}}{2} < \alpha_x^2 \leq \zeta_{H3}$, it follows from Eqs (2.19)–(2.22) that $\zeta_{H1} > 0$, $\zeta_{H1}\leq \zeta_{H2}$, $I_H > 0$, $\zeta_{H1} < I_H$, and $I_H\leq\zeta_{H2}$. Thus, $\bar{x}_{0H+} > 0$ if and only if $I > I_H$.

For $\zeta_{H3} < \alpha_x^2$ are $\zeta_{H1} > 0$, it follows from Eqs (2.19)–(2.23) that $\zeta_{H1} > \zeta_{H2}$, $I_H > 0$, $\zeta_{H1} > I_H$, and $I_H\leq\zeta_{H2}$. Thus, $\bar{x}_{0H+} > 0$ if and only if $I > \zeta_{H2}$.

Similar, $\bar{x}_{0H-}$ is real if and only if $I \geq I_H$ and $\zeta_{H1} < I < \zeta_{H2}$.

For $\alpha_x^2 \leq \frac{\zeta_{H3}}{4}$, it follows from Eqs (2.19) and (2.20) that $\zeta_{H1} < 0 < \zeta_{H2}$, and $I_H\leq 0$. Thus, $\bar{x}_{0H-} > 0$ if and only if $I < \zeta_{H2}$.

For $\frac{\zeta_{H3}}{4} < \alpha_x^2 \leq \frac{\zeta_{H3}}{2}$, it follows from Eqs (2.19), (2.20) and (2.23) that $\zeta_{H1}\leq 0 < \zeta_{H2}$, $I_H > 0$, and $I_H < \zeta_{H2}$. Thus, $\bar{x}_{0H-} > 0$ if and only if $I_H\leq I < \zeta_{H2}$.

For $\frac{\zeta_{H3}}{2} < \alpha_x^2 \leq \zeta_{H3}$, it follows from Eqs (2.19)- (2.21) and (2.23) that $0 < \zeta_{H1} < \zeta_{H2}$, $I_H > 0$, and $I_H < \zeta_{H2}$. Thus, $\bar{x}_{0H-} > 0$ if and only if $I_H\leq I < \zeta_{H2}$.

For $\zeta_{H3} < \alpha_x^2$, it follows from Eqs (2.20) and (2.21) that $\zeta_{H1} > \zeta_{H2}$ and $I_H > 0$. Thus, $\bar{x}_{0H-} > 0$ if and only if $\zeta_{H1} < I < \zeta_{H2}$, which is a contradiction.

The conditions for the existence of the hematopoietic steady states are summarized in Table 1.

Malignant steady states may exist depending on therange of the parameters. As above we chose the inflammatory stimuli $I$ as our leading parameter and make a complete analysis of possible malignant steady states. The analysis of the existence of the malignant steady states depends crucially on the following lumped parameters,

the last two being positive.

Proposition 3. Two malignant steady states $E_{L\pm}$ may exist in the following cases,

● If $\alpha_y^2 \leq \frac{\zeta_{L3}}{4}$ then $E_{L+}$ exists if and only if $I > 0$.

● If $\frac{\zeta_{L3}}{4} < \alpha_y^2 \leq \zeta_{L3}$ then $E_{L+}$ exists if and only if $I \geq I_{L}$.

● If $\zeta_{L3} < \alpha_y^2$ then $E_{L+}$ exists if and only if $I > \zeta_{L2}$.

● If $\alpha_y^2 \leq \frac{\zeta_{L3}}{4}$ then $E_{L-}$ exists if and only if $I \leq \zeta_{L2}$.

● If $\frac{\zeta_{L3}}{4} < \alpha_y^2 \leq \zeta_{L3}$ then $E_{L-}$ exists if and only if $I_{L} \leq I \leq \zeta_{L2}$.

● If $\zeta_{L3} < \alpha_y^2$ then $E_{L-}$ does not exist.

In case of existence, $E_{L+} = (0, 0, \bar{y}_{0L+}, \bar{y}_{1L+}, \bar{a}_{L+}, \bar{s}_{L+})$ is given by $\bar{y}_{0L+} = \frac{1}{2 e_s c_{yy}\alpha_y} \left(I-\zeta_{L1} + \sqrt{(\zeta_{L1}-I)^2-4 e_s \alpha_y(\zeta_{L2}-I)}\right)$, $\bar{s}_{L+} = \alpha_y (1+c_{yy} \bar{y}_{0L+})$, $\bar{a}_{L+} = \frac{\beta_y \bar{y}_{0L+}}{e_a \bar{s}_{L+}}$, and $\bar{y}_{1L+} = \frac{a_x A_x \bar{y}_{0L+}}{d_{y1}}$ whereas $E_{L-} = (0, 0, \bar{y}_{0L-}, \bar{y}_{0L-}, \bar{a}_{L-}, \bar{s}_{L-})$ is given by $\bar{y}_{0L-} = \frac{1}{2 e_s c_{yy}\alpha_y} \left(I-\zeta_{L1} - \sqrt{(\zeta_{L1}-I)^2-4 e_s \alpha_y(\zeta_{L2}-I)}\right)$, $\bar{s}_{L-} = \alpha_y (1+c_{yy} \bar{y}_{0L-})$, $\bar{a}_{L-} = \frac{\beta_y \bar{y}_{0L-}}{e_a \bar{s}_{L-}}$, and $\bar{y}_{1L-} = \frac{a_x A_x \bar{y}_{0L-}}{d_{y1}}$.

Proof. Due to symmetry in indices $x$ and $y$, the proof for the malignant case is equivalent to that for the hematopoietic case except index $H$ has to be substituted by $L$.

The result is summarized in Table 2.

The existence of a co-existing steady state

is far more cumbersome to deal with, since a wealth of sub-cases may arrive depending on various inequality-relations between the parameters. To avoid many tedious but straight forward calculations we limit ourself to the non-degenerate cases where $\zeta_{C1} = \alpha_yc_{yx}- \alpha_xc_{xx} \neq 0$ and $\zeta_{C2} = \alpha_yc_{yy}- \alpha_xc_{xy} \neq 0$.

From Eqs (2.4a and 2.4b), a linear relation between $\bar{x}_0$ and $\bar{y}_0$ directly follows,

where $\zeta_{C3} = \alpha_x-\alpha_y$. Thus, for the non-degenerate cases,

which geometrically corresponds to a straight line through $\left(0, \frac{\zeta_{C3}}{\zeta_{C2}}\right)$ and $\left(\frac{\zeta_{C3}}{\zeta_{C1}}, 0\right)$. Hence, two generic cases arrive, for $\left(0, \frac{\zeta_{C1}}{\zeta_{C2}}\right)$ corresponding to positive slope, $\frac{\zeta_{C1}}{\zeta_{C2}} < 0$ corresponding to negative slope, $\frac{\zeta_{C1}}{\zeta_{C2}} > 0$. The first case defines a half line in the positive octahedron and in this case $\bar{x}_{0C}\in(\max\{0, \frac{\zeta_{C3}}{\zeta_{C1}}\}; \infty)$ and $\bar{y}_{0C}\in(\max\{0, \frac{\zeta_{C3}}{\zeta_{C2}}\}; \infty)$. The second case corresponds to either no admissible solution (if and only if $\frac{\zeta_{C3}}{\zeta_{C2}} < 0$ and $\frac{\zeta_{C3}}{\zeta_{C1}} < 0$) or a line segment in the positive octahedron which requires that $\frac{\zeta_{C3}}{\zeta_{C2}} > 0$ and $\frac{\zeta_{C3}}{\zeta_{C1}} > 0$ and in that case are $\bar{x}_{0C}\in(0, \frac{\zeta_{C3}}{\zeta_{C1}})$ and $\bar{y}_{0C}\in(0, \frac{\zeta_{C3}}{\zeta_{C2}})$. From Eq (2.4a) and (2.29),

with $m_0 = \alpha_x (c_{xy}\frac{\zeta_{C3}}{\zeta_{C2}}+1)$ and $m_1 = \alpha_x (c_{xx}-c_{xy}\frac{\zeta_{C1}}{\zeta_{C2}})$. Before continuing, it is emphasized that $\zeta_1$, $\zeta_2$, $\zeta_3$, $m_0$, and $m_1$ all are independent of $I$ but may be positive, negative or in case of $m_0$ and $m_1$ zero. From Eq (2.7) it follows that a real and positive $\bar{s}$ exist for $(\bar{x}_{0C}, \bar{y}_{0C})\in \mathbb{R}_+ \times \mathbb{R}_+$,

where $\zeta_0 = \frac{4 r_s e_s}{e_a} > 0$. Similarly, a negative real root exists. Substituting (2.29) into (2.31) give,

Combining Eq (2.30) and (2.32) results in,

where $n_0 = 2e_sm_0$, $n_1 = 2e_sm_1$, $n_2 = \frac{\zeta_0 \beta_y \zeta_{C3}}{\zeta_{C2}}$ and $n_3 = \zeta_0 \left(\beta_x-\beta_y\frac{\zeta_{C1}}{\zeta_{C2}}\right)$.

Note that Eq (2.33) has no real solution if either of $f(x) = n_0-I+n_1 x$ and $g(x) = I^2+n_2+n_3 x$ are negative. Thus, if both $f(x)$ and $g(x)$ are positive, Eq (2.33) is equivalent to,

which may have up to two real positive solutions. Hence, there can be at most two coexistence steady states. More specifically,

is positive if and only if $f(x_{0C+}) > 0$, $g(x_{0C+}) > 0$, and

Similar,

is positive if and only if $f(x_{0C-}) > 0$, $g(x_{0C-}) > 0$, and

Note, some possibilities of equality signs in the inequalities are left out for simplification reasons. Equality may occur on a set of measure zero which is unlikely for a noisy biological system and including these possibilities makes the analysis much more messy. For practical purposes one may first calculate the two (possibly complex) roots $x$ of Eq (2.34) and afterwards examine whether these are real and positive, whether $f(x) > 0$ and $g(x) > 0$, and whether the corresponding $\bar{y}_{0C}$ calculated from Eq (2.29) is positive, thus the remaining component of $E_{C}$ will be positive too and the steady state admissible.

Continuing analytically is possible but becomes somehow cumbersome and instead we point out that for any choice of parameter values, there can be at most two coexistence steady states, their existence and value depending on the admissibility of $x_{0C+}$ (Eq (2.35)) and $x_{0C-}$ (Eq (2.37)).

3.   Stability and bifurcation analysis

In this section we analytically and numerically examine the stability properties of the various admissible steady states of Eq (2.1) in terms of selected parameters.

3.1.   Stability properties of the trivial steady state

The Jacobian of the trivial steady states $E_0$ is a triangular matrix and four of the six eigenvalues, $-d_{x1}, -d_{y1}, -e_s$ and $\frac{-I e_a}{e_s}$, are negative but the two, $\frac{r_x}{e_s} \left(I-e_s \alpha_x\right)$ and $\frac{r_y}{e_s} \left(I-e_s \alpha_y\right)$, may be positive, negative, zero. Thus, by the Hartman-Grobman Theorem ^[34]

Lemma 1. $E_0$ is asymptotically stable if $I < e_s \min \{\alpha_x, \alpha_y\}$, whereas it is unstable if $I > e_s \min \{\alpha_x, \alpha_y\}$.

At $E_{H\pm}$ the Jacobian for the hematopoietic states can be calculated (see Supplementary) and the resulting sixth order characteristic equation shows that $E_{H\pm}$ are stable for

However, this is not the generic case, since $\alpha_y < \alpha_x$ (and $c_{yx} \leq c_{xx}$), which contradicts $\bar{s}_H = \alpha_x (1+c_{xx}\bar{x}_{0H})$. Intensive numerical investigations shows that $E_{H\pm}$ are unstable.

The stability of $E_L$ is similar to that for the hematopoietic steady state except that it is stable if

which is fulfilled in the generic case, since $\alpha_x < \alpha_y$ (and $c_{xy} \leq c_{yy}$). This follows from $\bar{s}_L = \alpha_y (1+c_{yy}\bar{y}_{0L})$. The Jacobian may be found in supplementary.

Lastly, consider the co-existing steady state. The Jacobian at $E_{C}$ may be found in supplementary. However, it is hard to prove any result analytically and we therefore do the stability investigation numerically the in next section.

3.1.1.   Numerical Simulations and treatment scenarios

In this section, we focus on numerical results. The default values of parameters used in Figure 2 are given in Table 3. The values are the same as given in ^[18].

The model has been investigated for various choices of parameter values. In Figure 2, clusters of five important parameters, $\mathcal{C} = \frac{c_{xx}}{c_{yy}}$, $\mathcal{R} = \frac{\zeta_{H2}}{\zeta_{L2}} = \frac{\alpha_x}{\alpha_y}$ and $I$ are considered to investigate the number of steady states and their stability. In the default case $\mathcal{R} > 1$ (Figure 2a), a trivial steady state always exists, and for low inflammation, i.e., $I < \zeta_{L2}$ it is stable otherwise it is unstable. For $I > \zeta_{L2}$, a purely malignant steady state becomes admissible. For values of $I$ where the trivial and the malignant steady states are admissible, the malignant steady state is stable whereas the trivial steady state is unstable. An unstable hematopoietic steady state becomes admissible as $I$ becomes larger than the threshold value $\zeta_{H2}$, and increasing $I$ further causes emergence of a stable co-existing steady state while the malignant steady state becomes unstable. Thus, for $I > \zeta_{H2}$ and $\mathcal{C}$ sufficiently small, four steady states appear namely the trivial, the hematopoietic, the malignant and the co-existing steady states where the co-existing steady state is stable and the rest are unstable. This illustrates that the co-existing steady state depends on $I$, $\mathcal{C}$ and $\mathcal{R}$. Increasing $\mathcal{C}$ from a small, initial value makes the co-existing steady state vanish and the malignant steady state becomes stable whereas the trivial and the hematopoietic steady states remain unstable. Secondly, consider the second case where $\mathcal{R} = 1$ implying that $\zeta_{H2} = \zeta_{L2}$ (Figure 2b). Increasing $I$ across this value generates an unstable hematopoietic steady state and a malignant steady state simultaneously. For $\mathcal{C} < 1$ the malignant steady state is unstable, and a stable coexistence steady state is created as $I$ increase past $\zeta_{H2}$. For $\mathcal{C} > 1$ no coexistence steady state is created, instead the malignant steady state is stable. Hence, for $\mathcal{R} = 1$, decreasing $\mathcal{C}$ may change the topology from a stable malignant steady state to a stable coexistence steady state i.e. improving the prognosis from disease escape to disease equilibrium. The stable co-existing steady state bifurcates from the trivial steady state and remains stable until $\mathcal{C} = 1$. As $\mathcal{C}$ exceeds $1$, the co-existing steady state disappears, the malignant steady state becomes stable and the trivial and the hematopoietic steady state become unstable.

In the remaining panels, $\mathcal{R} < 1$, which implies that a stable hematopoietic steady state is created as the first transition to appear when increasing $I$ from low values past the threshold value $\zeta_{H2}$. Simultaneously, the trivial steady state becomes unstable. In Figure 2c where $\mathcal{R} = 0.97$ the hematopoietic steady state remains stable for low values of $\mathcal{C}$ until $I$ passes a threshold value where a stable coexistence steady state is created leaving the hematopoietic steady state unstable.

For larger values of $\mathcal{C}$ there is no coexistence steady state. Instead, as $I$ is increased, a region of bistability appears with a stable hematopoietic steady state and a stable malignant steady state. Increasing $I$ further the hematopoietic steady state becomes unstable. Hence, to reduce disease load, in the case of $\mathcal{R} < 1$, and large values of $\mathcal{C}$ and $I$, it may be optimal treatment to reduce the $\mathcal{C}$ value prior to reducing the inflammatory level to avoid being stuck in the bassin of attraction of the malignant steady state.

In Figure 2d where $\mathcal{R} = 0.93$, the coexistence steady state no longer appears, the region of bistability has shrunk and a hematopoietic stable steady state is more dominant.

In Figure 2e and f, $\mathcal{R}$ is decreased to $0.77$ and $0.5$ respectively, and the bistability region is no longer visible. For $I > \zeta_{H2}$ a hematopoietic steady state is the only stable steady state. Figure 2 indicates that reducing $\mathcal{C}$ and $\mathcal{R}$ should be targets of intervention. A reduction of $I$ may improve prognosis as well, for example for parameter values as in 2c.

4.   Discussion and conclusion

A mechanism-based model published in ^[5] - the Cancitis model - describing the interaction of the hematopoietic cells, malignant cells and inflammation is analysed here. A thorough mathematical investigation of the model is presented in this paper which did not appear previously. We conducted an analytical analysis of the steady states and showed that four kinds of steady states may exist i.e. trivial, hematopoietic, malignant and co-existing steady states. We characterized the stability of each of these steady states and identified transitions conditions in the number of steady states and in their stability. Trivial, hematopoietic, malignant and coexistence steady states all appear for some parameter values. The steady states are highly relevant as all trajectories appear to approach a steady state after some time - see Figure 3. The case of bistability is visualized in the bottom right panel of Figure 3, with the basin of attraction shown in the $(x_0, y_0)$-plane using initial condition $(x_1, y_1, a, s) = (4\times 10^{11}, 4\times 10^{11}, 600, 2)$. The initial conditions for $x_0$ and $y_0$ are varied in a range $1-10^{5}$. The malignant steady state has a large bassin of attraction (region (i)), while region (ii) marks the bassin of attraction for the hematopoietic steady state.

The intuitive interpretation in most bio-medical literature attributes the main cause for cancer development to the frequency of stem cell division. Another main cause is the regulatory feedback that allows stem cells residing in niche to further divide into blood cell required in blood stream. Our investigation is in agreement with this perception and quantifies this intuitive concept. Furthermore, it shows that stem cell population is important to target in treatment to prevent disease progression.

In ^[14] and ^[15] a model without immune interaction is presented. The authors discuss a fraction similar to $\mathcal{R}$ and show that it is important for the dynamics of the system. It has been shown ^[15] that the leukemic cell load can be temporarily reduced if the growth of HSC is larger than that of leukemic cells for cytokine-dependent AML.

It is generally assumed that $c_{yy} \leq c_{xx}$ since malignant cells might be less sensitive to environmental crowding ^[22] and ^[24]. The ratio $\mathcal{C}$ of inhibition of the hematopoietic relative to malignant cells is one of several important prognostic markers. For large values of $I$, bi-stable and mono-stable regions depend upon $\mathcal{C}$. It can be observed in Figure 2 that for small values of $\mathcal{C}$, i.e., $c_{yy} \geq c_{xx}$, either the hematopoietic steady state is stable or the co-existing steady state is stable which can be interpreted as a good prognosis. However, large values of $\mathcal{C}$ may lead to a worse situation, e.g. in one case, the malignant steady state is stable or there exists bi-stability of the hematopoietic and the malignant steady states (see Figure 2c). In addition to the ratio of inhibitive niche feedback, the ratio $\mathcal{R}$ is also important to consider, since it determines how robust the hematopoietic condition may be and how disastrously a potential blood cancer disease will develop. Thus for $\mathcal{R} > 1$ we have a more serious situation than for $\mathcal{R} < 1$ showing that if this reproduction ratio exceeds the threshold $\mathcal{R}_0 = 1$, it is more disastrous than if it is below $\mathcal{R}_0$.

The $JAK2$V617F allele burden is expected to increase due to the expansion of malignant cells. The $JAK2$V617F allele burden is interpreted as the ratio of malignant cells to the total number of mature cells. The model predicted $JAK2$V617F allele burden is shown in Figure 4 for the region where $E_C$ is stable. Perturbation of a parameter may improve or impair prognosis when the coexistence point is the stable attractor. The top panel of Figure 4, shows that decreasing $\mathcal{C}$ and $\mathcal{R}$ improve prognosis by lowering the allele burden. Contrarily, increasing $I$, causes a decay in allele burden. This suggests that inflammatory inhibitors could counteract treatments in this case. In other cases, increasing $I$ typically leads to a worse prognosis, considering Figure 2.

The model presented here may inform clinical practice to make group specific treatment protocols with particular focus on the inflammatory components which may accelerate or dampen the disease progression. Interventions should address decreasing $\mathcal{C}$ and $\mathcal{R}$ and potentially $I$ but the latter depends on the remaining parameter values as adverse effects may be observed.

Conflict of interest

No conflict of interest.

Supplementary

Stability analysis of Steady states:

At $E_{H\pm}$ the Jacobian of the purely hematopoietic steady state becomes,

where

$a_{11} = r_x \left(\frac{\bar{s}_{H}}{(1+c_{xx} \bar{x}_{0H})^2}-\alpha_x\right)$,

$a_{13} = -\frac{r_x c_{xy} \bar{s}_{H} \bar{x}_{0H}}{(1+c_{xx} \bar{x}_{0H})^2}$,

$a_{16} = \frac{r_x \bar{x}_{0H}}{(1+c_{xx} \bar{x}_{0H})}$,

$a_{21} = a_x A_x$,

$a_{22} = -d_{x1}$,

$a_{33} = r_y \left(\frac{\bar{s}_{H}}{1+c_{yx} \bar{x}_{0H}}-\alpha_y\right)$,

$a_{43} = a_y A_y$,

$a_{44} = -d_{y1}$,

$a_{51} = d_{x0}$,

$a_{52} = d_{x1}$,

$a_{53} = d_{y0}$,

$a_{54} = d_{y1}$,

$a_{55} = -e_a \bar{s}_{H}$,

$a_{56} = -e_a \bar{a}_{H}$,

$a_{65} = r_s$,

$a_{66} = -e_s$,

and rest of the elements of $J_{E_H}$ are zero.

At $E_{L\pm}$ the Jacobian of the purely malignant steady state.

where

$a_{11} = r_x\left(\frac{ \bar{s}_{L}}{1+c_{xy}\bar{y}_{0L}}-\alpha_x\right)$,

$a_{21} = a_x A_x$,

$a_{22} = -d_{x1}$,

$a_{31} = -\frac{r_y c_{yx}\bar{s}_{L} \bar{y}_{0L}}{(1+c_{yy} \bar{y}_{0L})^2}$,

$a_{33} = r_y \left(\frac{\bar{s}_{L}}{(1+c_{yy}\bar{y}_{0L})^2}-\alpha_y\right)$,

$a_{36} = \frac{r_y \bar{y}_{0L}}{1+c_{yy}\bar{y}_{0L}}$.

$a_{43} = a_y A_y$.

$a_{44} = -d_{y1}$.

$a_{51} = d_{x0}$.

$a_{52} = d_{x1}$.

$a_{53} = d_{y0}$.

$a_{54} = d_{y1}$.

$a_{55} = -e_a \bar{s}_{L}$.

$a_{56} = -e_a \bar{a}_{L}$.

$a_{65} = r_s$.

$a_{66} = -e_s$

and rest of the elements of array are zero.

At $E_{C\pm}$ the Jacobian of the co-existing steady state becomes,

where

$a_{11} = r_x \left(\frac{\bar{s}_C}{(1+c_{xx}\bar{x}_{0C}+c_{xy} \bar{y}_{0C})^2}-\alpha_x \right)$,

$a_{13} = -\frac{r_x c_{xy} \bar{s}_C \bar{x}_{0C}}{(1+c_{xx}\bar{x}_{0C}+c_{xy}\bar{y}_{0C})^2}$,

$a_{16} = \frac{r_x \bar{x}_{0C}}{(1+c_{xx}\bar{x}_{0C}+c_{xy} \bar{y}_{0C})}$,

$a_{21} = a_x A_x$,

$a_{22} = -d_{x1}$,

$a_{31} = -\frac{r_y c_{yx} \bar{s}_C \bar{y}_{0C}}{(1+c_{yx}\bar{x}_{0C}+c_{yy}\bar{y}_{0C})^2}$,

$a_{33} = r_y \left(\frac{\bar{s}_C}{(1+c_{yx}\bar{x}_{0C}+c_{yy} \bar{y}_{0C})^2}-\alpha_y \right)$,

$a_{36} = \frac{r_y \bar{y}_{0C}}{(1+c_{yx}\bar{x}_{0C}+c_{yy} \bar{y}_{0C})}$,

$a_{43} = a_y A_y$,

$a_{44} = -d_{y1}$,

$a_{51} = d_{x0}$,

$a_{52} = d_{x1}$,

$a_{53} = d_{y0}$,

$a_{54} = d_{y1}$,

$a_{55} = -e_a \bar{s}_C$,

$a_{56} = -e_a \bar{a}_C$,

$a_{65} = r_s$,

$a_{66} = -e_s$

and rest of the elements of the $J_{E_C}$ are zero.