A structured population model of clonal selection in acute leukemias with multiple maturation stages

TL;DR

This paper develops a structured population model to understand clonal selection in acute leukemias, highlighting the role of self-renewal fractions in driving clonal dominance, and providing analytical and numerical insights based on patient data.

Contribution

It introduces a multi-compartmental, analytically tractable model of leukemia clonal dynamics that incorporates multiple maturation stages and is calibrated with real patient data.

Findings

Clonal selection is primarily driven by self-renewal fractions of leukemic stem cells.

Proliferation rates of non-stem cells have limited impact on clonal selection.

The model offers a formal understanding of how clonal dominance emerges in leukemia.

Abstract

Recent progress in genetic techniques has shed light on the complex co-evolution of malignant cell clones in leukemias. However, several aspects of clonal selection still remain unclear. In this paper, we present a multi-compartmental continuously structured population model of selection dynamics in acute leukemias, which consists of a system of coupled integro-differential equations. Our model can be analysed in a more efficient way than classical models formulated in terms of ordinary differential equations. Exploiting the analytical tractability of this model, we investigate how clonal selection is shaped by the self-renewal fraction and the proliferation rate of leukemic cells at different maturation stages. We integrate analytical results with numerical solutions of a calibrated version of the model based on real patient data. In summary, our mathematical results formalise the biological notion that clonal selection is driven by the self-renewal fraction of leukemic stem cells and the clones that possess the highest value of this parameter are ultimately selected. Moreover, we demonstrate that the self-renewal fraction and the proliferation rate of non-stem cells do not have a substantial impact on clonal selection. Taken together, our results indicate that interclonal variability in the self-renewal fraction of leukemic stem cells provides the necessary substrate for clonal selection to act upon.

Equations

$${\rho }_i(t):={\int }_0^1{n}_i(t,x)\mathrm{d}x\text{mbox}withi=1,\dots ,M.$$

$$a_i\in {\mathrm{W}}^{1,\infty }\left([0,1]\right)\text{mbox}with{a}_i:[0,1]\to \left(\frac{1}{2},1\right)\text{mbox}foralli=1,\dots ,M-1,$$

$$p_i\in {\mathrm{W}}^{1,\infty }\left([0,1]\right)\text{mbox}with{p}_i:[0,1]\to \left(0,1\right)\text{mbox}foralli=1,\dots ,M-1,$$

$$\frac{1}{2}<a_i(x)<a_1(x)<1\text{mbox}forallx\in [0,1],\text{mbox}foralli=2,\dots ,M-1.$$

$$\left\{\begin{array}[]{ll}\displaystyle{\frac{\partial}{\partial t}n_{1}(t,x)=\left(\frac{2\,a_{1}(x)}{1+K\rho_{M}(t)}-1\right)p_{1}(x)\,n_{1}(t,x),}\\ \\ \displaystyle{\frac{\partial}{\partial t}n_{i}(t,x)=2\left(1-\frac{a_{i-1}(x)}{1+K\rho_{M}(t)}\right)p_{i-1}(x)\,n_{i-1}(t,x)}\\ \\ \displaystyle{\phantom{\frac{\partial n_{i}}{\partial t}(t,x)=}\;+\left(\frac{2\,a_{i}(x)}{1+K\rho_{M}(t)}-1\right)p_{i}(x)\,n_{i}(t,x),\quad i=2,\ldots,M-1,}\\ \\ \displaystyle{\frac{\partial}{\partial t}n_{M}(t,x)=2\left(1-\frac{a_{M-1}(x)}{1+K\rho_{M}(t)}\right)p_{M-1}(x)\,n_{M-1}(t,x)-d\,n_{M}(t,x),}\\ \\ \displaystyle{\rho_{M}(t)=\int_{0}^{1}n_{M}(t,x)\,{\rm d}x,}\end{array}\right.$$

$$n_i(0,x)=n_i^0(x)\in \mathrm{C}([0,1]),0<n_i^0<\infty \text{mbox}on[0,1]\text{mbox}fori=1,\dots ,M,$$

$$\left\{\begin{array}[]{ll}\displaystyle{\varepsilon\,\frac{\partial}{\partial t}n_{1\varepsilon}(t,x)=P_{1}(\rho_{M\varepsilon}(t),x)\,n_{1\varepsilon}(t,x),}\\ \\ \displaystyle{\varepsilon\,\frac{\partial}{\partial t}n_{i\varepsilon}(t,x)=Q_{i-1}(\rho_{M\varepsilon}(t),x)\,n_{i-1\varepsilon}(t,x)}\\ \displaystyle{\phantom{\varepsilon\,\frac{\partial n_{i\varepsilon}}{\partial t}(t,x)=}\;+P_{i}(\rho_{M\varepsilon}(t),x)\,n_{i\varepsilon}(t,x),\quad i=2,\ldots,M-1,}\\ \displaystyle{\varepsilon\,\frac{\partial}{\partial t}n_{M\varepsilon}(t,x)=Q_{M-1}(\rho_{M\varepsilon}(t),x)\,n_{M-1\varepsilon}(t,x)-d\,n_{M\varepsilon}(t,x),}\\ \\ \displaystyle{\rho_{M\varepsilon}(t)=\int_{0}^{1}n_{M\varepsilon}(t,x)\,{\rm d}x,}\\ \\ n_{i\varepsilon}(0,x)=n_{i}^{0}(x),\quad i=1,\ldots,M,\end{array}\right.$$

$$P_i({\rho }_{M\epsilon },x):=\left(\frac{2a_i(x)}{1+K{\rho }_{M\epsilon }}-1\right)p_i(x)\text{mbox}fori=1,\dots ,M-1$$

$$Q_i({\rho }_{M\epsilon },x):=2\left(1-\frac{a_i(x)}{1+K{\rho }_{M\epsilon }}\right)p_i(x)\text{mbox}fori=1,\dots ,M-1.$$

$${\rho }_{i\epsilon }(t):={\int }_0^1{n}_{i\epsilon }(t,x)\mathrm{d}x\text{mbox}withi=1,\dots ,M$$

$$R_{i\epsilon }(t,x):={\int }_0^t{P}_i({\rho }_{M\epsilon }(s),x)\mathrm{d}s\text{mbox}fori=1,\dots ,M-1.$$

$$n_{i\epsilon }\underset{\epsilon \to 0}{\rightharpoonup }{n}_i\text{mbox}on{w}^{\ast }-{\mathrm{L}}^{\infty }\left((0,T),{\mathrm{M}}^1([0,1])\right)\text{mbox}fori=1,\dots ,M$$

$$R_{i\epsilon }\underset{\epsilon \to 0}{\to }{R}_i\text{mbox}uniformlyin[0,T]\times [0,1]\text{mbox}fori=1,\dots ,M-1,$$

$$R_i(t,x):={\int }_0^t{P}_i({\rho }_M(s),x)\mathrm{d}s\text{and}{R}_i\in {\mathrm{W}}^{1,\infty }\left((0,T)\times [0,1]\right).$$

$$\underset{x\in [0,1]}{\max }{R}_1(t,x)=0\text{mbox}foranyt\in [0,T]$$

$$\underset{x\in [0,1]}{\max }{R}_i(t,x)<0\text{mbox}foranyt\in [0,T],\text{mbox}fori=2,\dots ,M-1.$$

$$\text{supp}(n_i(t,\cdot ))=\underset{x\in [0,1]}{\mathrm{argmax}}{a}_1(x)\text{mbox}fori=1,\dots ,M.$$

$${\rho }_{i\epsilon }(t)<{\overline{\rho }}_i\text{mbox}fori=1,\dots ,M,\text{mbox}forallt\ge 0\text{mbox}andforany\epsilon >0$$

$$\text{mbox}forall(t,x)\in [0,\infty )\times [0,1]\text{we have}{R}_{i\epsilon }(t,x)\underset{\epsilon \to 0}{\to }{R}_i(t,x),$$

$$R_1(t,x)\le 0\text{mbox}forall(t,x)\in (0,T)\times [0,1]$$

$$R_i(t,x)<0\text{mbox}forall(t,x)\in (0,T)\times [0,1],\text{mbox}foralli=2,\dots ,M-1.$$

$${\int }_0^1{n}_{1\epsilon }(t,x)\mathrm{d}x={\int }_0^1{n}_1^0(x)e^{\frac{R_{1\epsilon }(t,x)}{\epsilon }}\mathrm{d}x\ge e^{\frac{\sigma }{\epsilon }}{\int }_{\hat{x}-\sigma }^{\hat{x}+\sigma }{n}_1^0(x)\mathrm{d}x\underset{\epsilon \to 0}{\to }\infty$$

$${\int }_0^t{P}_1({\rho }_M(s),x)\mathrm{d}s=p_1(x){\int }_0^t\left(\frac{2a_1(x)}{1+K{\rho }_M(s)}-1\right)\mathrm{d}s\le 0$$

$${\int }_0^t\left(\frac{2a_1(x)}{1+K{\rho }_M(s)}-1\right)\mathrm{d}s\le 0\text{mbox}forall(t,x)\in (0,T)\times [0,1].$$

$${\int }_0^t\left(\frac{2a_i(x)}{1+K{\rho }_M(s)}-1\right)\mathrm{d}s<0\text{mbox}forall(t,x)\in (0,T)\times [0,1]$$

$$p_i(x){\int }_0^t\left(\frac{2a_i(x)}{1+K{\rho }_M(s)}-1\right)\mathrm{d}s<0\text{mbox}forall(t,x)\in (0,T)\times [0,1]$$

$$\text{mbox}if{R}_1(t,\cdot )<0\text{mbox}on[0,1]\text{mbox}then{n}_i(t,\cdot )=0\text{mbox}a.e.on[0,1]$$

$$\underset{\epsilon \to 0}{\lim }{\int }_{\hat{t}-\sigma }^{\hat{t}+\sigma }{\int }_{\hat{x}-\sigma }^{\hat{x}+\sigma }{n}_{1\epsilon }(t,x)\mathrm{d}x\mathrm{d}t$$

$$\underset{\epsilon \to 0}{\lim }{\int }_{\hat{t}-\sigma }^{\hat{t}+\sigma }{\int }_{\hat{x}-\sigma }^{\hat{x}+\sigma }{n}_{1\epsilon }(t,x)\mathrm{d}x\mathrm{d}t=0.$$

$${\int }_{\hat{t}-\sigma }^{\hat{t}+\sigma }{\int }_0^1\phi (x)n_1(t,x)\mathrm{d}x\mathrm{d}t=\underset{\epsilon \to 0}{\lim }{\int }_{\hat{t}-\sigma }^{\hat{t}+\sigma }{\int }_0^1\phi (x)n_{1\epsilon }(t,x)\mathrm{d}x\mathrm{d}t$$

Full text

1 Introduction

Statement of the biological problem. Leukemias are malignant diseases of the hematopoietic system. Similarly to the healthy hematopoietic system, the leukemic cell bulk is organised as a hierarchy of multiple maturation stages (i.e. maturation compartments) – from stem cells through a number of increasingly mature progenitor cells to the most mature cells (Bonnet and Dick, 1997, Hope et al., 2004). Red blood cells, white blood cells and platelets constitute the most mature compartment of healthy cells, whereas non-functional leukemic blasts are the most mature leukemic cells. The extensive growth of non-functional leukemic blasts leads to impaired hematopoiesis and causes a shortage of healthy blood cells.

In contrast to the most mature cells, which do not divide and die at a constant rate, stem and progenitor cells can proliferate and give rise to progeny cells which are either at the same maturation stage as their parent (self-renewal) or at a subsequent maturation stage (differentiation). These processes can be quantitatively characterised in terms of two parameters: the cell proliferation rate, i.e. the number of cell divisions per unit of time, and the cell self-renewal fraction, i.e. the probability that a progeny adopts the same cell fate as its parent (Marciniak-Czochra et al., 2009, Stiehl and Marciniak-Czochra, 2012). There is both theoretical (Stiehl et al., 2015) and experimental (Jung et al., 2015, Metzeler et al., 2013, Wang et al., 2017) evidence suggesting that the proliferation rate and the self-renewal fraction of leukemic cells have a significant impact on the disease dynamics and on patient prognosis (Stiehl and Marciniak-Czochra, 2017).

The collection of all stem, progenitor and most mature cells which carry the same set of genetic alterations defines a leukemic clone. Recent experimental evidence indicates that the leukemic cell bulk of an individual patient is composed of multiple clones carrying different mutations (Anderson et al., 2011, Ding et al., 2012, Ley et al., 2008) and having different functional properties, among which different proliferation rates and self-renewal fractions (Eppert et al., 2011, Heuser et al., 2009). Such clonal heterogeneity poses a major obstacle to successful therapy and management of disease relapse (Choi et al., 2007, van Delft et al., 2011, Lutz et al., 2013a, b). In fact, it has been reported that in many cases of acute lymphoblastic leukemias (ALL) and acute myeloid leukemias (AML) relapse is triggered by the selection of clones that have been present as minor (small) clones at the time of diagnosis rather than by newly acquired mutations (Choi et al., 2007, van Delft et al., 2011, Ding et al., 2012, Jan and Majeti, 2013).

It is becoming increasingly apparent that clonal heterogeneity results from a selection process at the cellular level which triggers the expansion of some clones and the out-competition of others (Belderbos et al., 2017, Choi et al., 2007, Hirsch et al., 2016, Noetzli et al., 2017, Stiehl et al., 2014a, Wu et al., 2017). However, several aspects of clonal selection are so far not well understood and a number of important questions remain open (Stiehl and Marciniak-Czochra, 2019). In this paper, we focus on clonal selection taking place before diagnosis or relapse of acute leukemias. We use a mathematical modelling approach to address the following questions:

• Q1

What is the role of the proliferation rate and the self-renewal fraction of leukemic cells in clonal selection?

• Q2

Can we observe clonal selection among clones that have identical stem-cell properties (in terms of proliferation rate and self-renewal fraction) and differ in their progenitor-cell properties?

• Q3

What are the necessary conditions (in terms of proliferation rate and self-renewal fraction) for the long-term coexistence of different clones?

Mathematical framework. A well-established mathematical approach to describing the dynamics of multiple leukemic clones is to use models with multiple compartments formulated in terms of ordinary differential equations (ODEs) (Stiehl et al., 2014a, Werner et al., 2013). In these models, every cell is characterised by a pair of indices $i=1,\ldots,M$ and $j=0,\ldots,J$. The index $i$ corresponds to the maturation stage of the cell (i.e. the compartment to which the cell belongs). The index $j=1,\ldots,J$ indicates what leukemic clone the cell belongs to, whilst the index $j=0$ is conventionally associated to healthy cells. In this modelling framework, a collection of parameters $p^{j}_{i}$ and $a^{j}_{i}$ is introduced to model, respectively, the proliferation rate and the self-renewal fraction of cells of clone $j=1,\ldots,J$ at the maturation stage $i=1,\ldots,M-1$. Clonal heterogeneity is incorporated into the model by allowing the values of these parameters to change from one clone to the other. The dynamic of every clone is described by a system of $M$ coupled ODEs that track the time evolution of the density of cells at different maturation stages. An additional system of $M$ coupled ODEs is introduced to model the dynamic of the healthy cells. Coherently with biological findings (Kondo et al., 1991, Layton et al., 1989, Shinjo et al., 1997), all ODEs of the model are coupled through a feedback signal that regulates the cell dynamics and depends on the total density of cells at the maturation stage $M$ (i.e. all the most mature healthy and leukemic cells). A prototypical version of such ODE models can be found in Appendix A of this paper.

These models consist of systems of $(J+1)\,M$ coupled ODEs. Due to the high degree of clonal heterogeneity usually observed in leukemia patients, the biologically realistic values of $J$ can be very high. As a result, the analysis of these models becomes very hard (if not impossible) in scenarios that are clinically relevant. This poses limitations to the robustness of the biological conclusions that can be obtained using these models. To overcome these limitations, here we present a modelling framework whereby the index $j$ is replaced by a continuous structuring variable $x$. This variable can be seen as a parameterisation of the self-renewal fraction and the proliferation rate of the different clones. Hence, in our modelling framework the parameters $p^{j}_{i}$ and $a^{j}_{i}$ of the ODE model are replaced by some functions $p_{i}(x)$ and $a_{i}(x)$ which represent, respectively, the proliferation rate and the self-renewal fraction of cells that are at the maturation stage $i$ and belong to the clone identified by the variable $x$. This multi-compartmental continuously structured population model consists of a system of $M$ coupled integro-differential equations (IDEs) that can be analysed in a more efficient way than its ODE counterpart.

Main results of this paper. Exploiting the analytical tractability of the model, we address questions Q1-Q3 listed above by elucidating how clonal selection is shaped by the self-renewal fraction and the proliferation rate of leukemic cells at different maturation stages. Analytical results are integrated with numerical solutions of a calibrated version of the model based on patient data from the existing literature. In summary, our results formalise the idea that clonal selection is controlled by the self-renewal fraction of leukemic stem cells and the clones with the highest value of this parameter are ultimately selected. This implies that only the clones whose stem cells are characterised by the highest self-renewal fraction can stably coexist. Finally, our results indicate that interclonal variability in the self-renewal fraction of leukemic stem cells provides the necessary substrate for clonal selection to act upon.

Our analytical work follows earlier papers on the asymptotic analysis of IDEs that arise in the mathematical modelling of selection dynamics in populations structured by physiological traits  (Barles et al., 2009, Busse et al., 2016, Calsina et al., 2013, Chisholm et al., 2016, Delitala and Lorenzi, 2012, Desvillettes et al., 2008, Diekmann et al., 2005, Lorenzi et al., 2014, Lorz et al., 2011, Perthame and Barles, 2008, Raoul, 2011). In particular, Busse et al. (2016) have studied a basic version of the model, in which only two maturation stages are considered (i.e. $M=2$) and all the clones have the same cell proliferation rate, i.e. the functions $p_{1}(x),\ldots,p_{M-1}(x)$ are constant. The main novelty of our work is that we let the number of maturation stages $M$ be arbitrarily large and we allow the cell proliferation rate to vary from clone to clone, i.e. the functions $p_{1}(x),\ldots,p_{M-1}(x)$ are not necessarily constant. This makes the application domain of our results significantly wider and strengthens the robustness of our biological conclusions. Due to these additional layers of complexity, our analysis builds on a method of proof which is different from that proposed by Busse et al. and is based on asymptotic arguments related to those proposed by Desvillettes et al. in (Desvillettes et al., 2008). We remark that our analytical results rely on general assumptions and, therefore, are applicable to different subtypes of leukemia.

2 Description of the model

We present a multi-compartmental continuously structured population model for the dynamics of cells of multiple leukemic clones and healthy cells at different maturation stages. The model is given by the system of IDEs (5) and is defined as a continuous version of the multi-compartmental ODE model presented in (Stiehl et al., 2014a) – for the sake of completeness, we provide a short description of such an ODE model in Appendix A – and a generalisation of the continuously-structured population model with two maturation compartments considered in (Busse et al., 2016). The key ideas underlying our model are summarised by the schematic diagrams in Fig. 1.

As illustrated by the scheme in Fig. 1(B), we capture the high degree of clonal heterogeneity usually observed in leukemia patients by introducing a continuous structuring variable $x\in[0,1]$, and we assume that $x=0$ corresponds to healthy cells whereas different leukemic clones are characterised by different values of $x\in(0,1]$.

At time $t\in[0,T]$, where $T>0$ is an arbitrary final time, the population densities of stem cells (compartment $i=1$), increasingly mature progenitor cells (compartments $i=2,\ldots,M-1$), and mature cells and leukemic blasts (compartment $i=M$) are represented by the functions $n_{i}(t,x)\geq 0$. At every time instant $t$, the total density of cells at the $i$-th maturation stage is defined as

[TABLE]

As illustrated by the scheme in Fig. 1(A), mature cells and leukemic blasts do not divide and die at a constant rate $d>0$. Moreover, we use the functions $p_{i}(x)$ and $a_{i}(x)$ to model, respectively, the proliferation rate and the self-renewal fraction of cells of clone $x$ at the maturation stage $i=1,\ldots,M-1$. The flux to mitosis of cells of clone $x$ at the maturation stage $i=1,\ldots,M-1$ at time $t$ is given by $n_{i}(t,x)\,p_{i}(x)$. These $n_{i}(t,x)\,p_{i}(x)$ cells divide into $2\,n_{i}(t,x)\,p_{i}(x)$ offspring cells. The fraction $a_{i}(x)\,s(t)$ of the offspring cells is at the maturation stage $i$ (self-renewal), while the fraction $1-a_{i}(x)\,s(t)$ is at the maturation stage $i+1$ (differentiation). The factor $s(t)$ models the concentration of feedback signal that promotes the self-renewal of dividing cells and is absorbed by mature cells and leukemic blasts. In agreement with the biological findings presented in (Kondo et al., 1991, Layton et al., 1989, Shinjo et al., 1997), we let $s(t)$ be a monotonically decreasing function of $\rho_{M}(t)$. In particular, we use the definition $s(t)=\frac{1}{1+K\rho_{M}(t)}$, where the parameter $K>0$ is related to the degradation rate of the feedback signal (Marciniak-Czochra et al., 2009, Stiehl et al., 2014c, Stiehl and Marciniak-Czochra, 2011, 2012). Such a definition can be derived from an ODE for $s(t)$ using a quasi-stationary approximation (Getto et al., 2013, Marciniak-Czochra et al., 2018).

Throughout the paper, we assume

[TABLE]

[TABLE]

[TABLE]

Assumptions (2)-(4) rely on the following biological considerations. It has been shown that the self-renewal fraction of stem cells has to be larger than $\frac{1}{2}$ to allow for cell expansion (Stiehl and Marciniak-Czochra, 2011, 2012). This justifies assumption (2) for $i=1$. Moreover, we justify assumptions (2) for $i=2,\ldots,M-1$ on the basis of biological evidence indicating that progenitor cells at different maturation stages are able to expand with little or no influx from the stem-cell compartment (Roelofs et al., 2003, Yamamoto et al., 2013). Previous experimental and theoretical studies (Hope et al., 2004, Reya et al., 2001, Stiehl and Marciniak-Czochra, 2011) have shown that stem cells have a higher self-renewal fraction than the progenitor cells to which they give rise. This justifies assumption (4). Furthermore, the cellular proliferation rates are bounded from above due to the time required by genome replication. A reasonable upper bound is between 1 and 2 divisions per day (Morgan et al., 2007), from which one can deduce assumptions (3).

The evolution of the cell population density functions $n_{i}(t,x)$ is governed by the following system of coupled IDEs

[TABLE]

subject to the initial conditions below

[TABLE]

which correspond to a biologically consistent scenario where numerous leukemic clones are present at the time of diagnosis.

3 Analysis of clonal selection

We prove a general asymptotic result (vid. Section 3.1) that sheds light on the way in which the self-renewal fraction and the proliferation rate of cells at different maturation stages impact on clonal selection (vid. Section 3.2).

3.1 A general asymptotic result

Building on previous studies on the long-time behaviour of continuously structured population models (Barles et al., 2009, Chisholm et al., 2016, Desvillettes et al., 2008, Diekmann et al., 2005, Lorz et al., 2011, Perthame and Barles, 2008), we introduce a small parameter $\varepsilon>0$ and use the time scaling $t\mapsto\frac{t}{\varepsilon}$, so that considering the asymptotic regime $\varepsilon\rightarrow 0$ is equivalent to studying the behaviour of the solutions to the IDEs of the model over many cell generations. With this scaling, the Cauchy problem for the cell population density functions $n_{i}\left(\frac{t}{\varepsilon},x\right)=n_{i\varepsilon}(t,x)$ reads as

[TABLE]

where

[TABLE]

and

[TABLE]

Throughout this section we use the notation

[TABLE]

and

[TABLE]

A general result on the asymptotic behaviour of the cell population density functions $n_{i\varepsilon}(t,x)$ for $\varepsilon\rightarrow 0$ (i.e. in the limit of many cell generations) is established by Theorem 3.1.

Theorem 3.1 (A general asymptotic result)

Under assumptions (2)-(4) and (6), the solutions to the Cauchy problem (7) are such that, up to extraction of subsequences,

[TABLE]

and

[TABLE]

where

[TABLE]

Moreover, the limits $R_{i}$ are such that

[TABLE]

and

[TABLE]

Finally, the limits $n_{i}$ are such that for a.e. $t\in[0,T]$

[TABLE]

Proof

We divide the proof of Theorem 3.1 into five parts.

Part 1: Non-negativity of $n_{i\varepsilon}(t,x)$ and continuity of $\rho_{i\varepsilon}(t)$. For all $\varepsilon>0$, standard arguments based on the Banach fixed point theorem allow one to prove that, under assumptions (2)-(4) and (6), the Cauchy problem (7) admits a solution of non-negative components $n_{i\varepsilon}\in{\rm C}([0,\infty),{\rm L}^{1}([0,1]))$.

Part 2: Uniform upper bounds for $\rho_{i\varepsilon}(t)$. The following uniform upper bounds hold

[TABLE]

with $0<\overline{\rho}_{i}<\infty$, the proof of which is provided in Appendix B.

Part 3: Proof of (9). The upper bounds (15) allow us to use the Banach-Alaoglu theorem to conclude that, up to extraction of subsequences, the asymptotic result (9) is verified.

Part 4: Proof of (10) and (11). The asymptotic result (9) on $n_{M\,\varepsilon}$ allows us to conclude that there exists a subsequence of $R_{i\varepsilon}$, that we denote again as $R_{i\varepsilon}$, such that

[TABLE]

for all $i=1,\ldots,M-1$, where $R_{i}(t,x)$ is defined according to (11). Furthermore, assumptions (2) and (3) ensure that, for any $\varepsilon>0$, the function $R_{i\varepsilon}$ and its first derivatives with respect to $t$ and $x$ are bounded in ${\rm L}^{\infty}\left((0,T)\times[0,1]\right)$. Therefore, $R_{i\varepsilon}$ belongs to ${\rm W}^{1,\infty}\left((0,T)\times[0,1]\right)$ for any $\varepsilon>0$ and, using the fact that ${\rm W}^{1,\infty}((0,T)\times[0,1])$ is compactly embedded in ${\rm C}([0,T]\times[0,1])$, we conclude that the uniform convergence result (10) is verified.

Part 5: Proof of (12)-(14). We prove the results (12)-(14) in four steps.

Part 5 – Step 1. We prove that

[TABLE]

and

[TABLE]

By contradiction, assume that there exists $(\hat{t},\hat{x})\in[0,T]\times(0,1)$ such that $R_{1}(\hat{t},\hat{x})>0$. The convergence result (10) for $i=1$ implies that $R_{1\varepsilon}(t,x)\geq\sigma$ for some $\sigma>0$, as long as $|t-\hat{t}|\leq\sigma$, $|x-\hat{x}|\leq\sigma$ and $\sigma\geq\varepsilon>0$. Since $n_{1}^{0}>0$ on $[0,1]$ we conclude that

[TABLE]

for all $t\in[\hat{t}-\sigma,\hat{t}+\sigma]$. This contradicts the upper bound (15) for $\rho_{1\varepsilon}(t)$, thus proving (17). Moreover, such a result on $R_{1}$ ensures that

[TABLE]

for all $(t,x)\in(0,T)\times[0,1]$ and, since $p_{1}>0$ on $[0,1]$ [cf. assumption (3)], the above inequality implies that

[TABLE]

Since $a_{1}(x)>a_{i}(x)$ for all $x\in[0,1]$ and for all $i=2,\ldots,M-1$ [cf. assumption (4)], the latter inequality allows one to conclude that for all $i=2,\ldots,M-1$

[TABLE]

and using the fact that $p_{i}>0$ on $[0,1]$ [cf. assumption (3)] one obtains

[TABLE]

for all $i=2,\ldots,M-1$. This concludes the proof of (18), which implies that (13) is verified.

Part 5 – Step 2. We prove that

[TABLE]

for all $i=1,\ldots,M$. Throughout this step we consider $(\hat{t},\hat{x})\in(0,T)\times(0,1)$ such that $R_{1}(\hat{t},\hat{x})<0$.

Proof of (19) for $i=1$. The uniform convergence result (10) for $i=1$ ensures that there exists some $\sigma>0$ such that $R_{1\varepsilon}(t,x)\leq-\sigma$ for $|t-\hat{t}|\leq\sigma$, $|x-\hat{x}|\leq\sigma$ and $\sigma\geq\varepsilon>0$. This allows us to conclude that

[TABLE]

that is,

[TABLE]

The weak convergence result (9) for $n_{1\varepsilon}(t,x)$ ensures that

[TABLE]

for every test function $\varphi:[0,1]\rightarrow\mathbb{R}$. Therefore, choosing a smooth test function $\varphi$ that satisfies the following conditions

[TABLE]

where $\mathbf{1}$ denotes the indicator function, and using the fact that $n_{1\varepsilon}$ is non-negative we find

[TABLE]

Substituting (20) into the latter integral inequality we conclude that

[TABLE]

for every smooth test function that satisfies (22). Hence, the result  (19) for $i=1$ is verified.

Proof of (19) for $i=2$. Multiplying by a test function $\varphi:[0,1]\rightarrow\mathbb{R}$ both sides of the IDE (7) for $n_{2\varepsilon}$ and integrating over the set $[\hat{t}-\sigma,\hat{t}+\sigma]\times[0,1]$ we find

[TABLE]

Since the uniform upper bound (15) for $\rho_{2\varepsilon}$ ensures that

[TABLE]

we conclude that

[TABLE]

Choosing a smooth test function that satisfies (22) and using (21) and (24) along with the fact that $Q_{1}(\rho_{M\,\varepsilon}(t),x)>0$ on $[0,T]\times[0,1]$ yields

[TABLE]

Since $P_{2}(\rho_{M\,\varepsilon}(t),x)<0$ on $[0,1]$ for a.e. $t\in(0,T)$ when $\varepsilon\to 0$, the result given by (25) allows us to conclude that

[TABLE]

for every smooth test function that satisfies (22). The weak convergence result (9) for $n_{2\varepsilon}(t,x)$ ensures that

[TABLE]

Hence, using (26) and (27) we conclude that

[TABLE]

for every smooth test function that satisfies (22). Therefore, the result (19) for $i=2$ is verified.

Proof of (19) for all $i=3,\ldots,M-1$. Using a bootstrap argument based on the method of proof that we have used for the case $i=2$, one can prove that for all $i=3,\ldots,M-1$

[TABLE]

for every smooth test function that satisfies (22). Hence, the results (19) for $i~{}=~{}3,\ldots,M-1$ are verified.

Proof of (19) for $i=M$. Multiplying by a test function $\varphi:[0,1]\rightarrow\mathbb{R}$ both sides of the IDE (7) for $n_{M\varepsilon}$ and integrating over the set $[\hat{t}-\sigma,\hat{t}+\sigma]\times[0,1]$ we obtain

[TABLE]

Since the uniform upper bound (15) for $\rho_{M\varepsilon}$ ensures that

[TABLE]

we conclude that

[TABLE]

Choosing a smooth test function that satisfies (22) and using the result (28) for $n_{M-1\varepsilon}(t,x)$ along with the fact that $Q_{M-1}(\rho_{M\,\varepsilon}(t),x)>0$ on $[0,T]\times[0,1]$ gives

[TABLE]

The weak convergence result (9) for $n_{M\varepsilon}(t,x)$ ensures that

[TABLE]

Hence, using (29) and (30) we conclude that

[TABLE]

for every smooth test function that satisfies (22). Therefore, the result (19) for $i=M$ is verified.

Part 5 – Step 3. We prove that

[TABLE]

[TABLE]

and

[TABLE]

We begin by noting that, since $p_{1}>0$ on $[0,1]$, if

[TABLE]

for some $t\in[0,T]$ then [vid. the definition (11) of the function $R_{1}$]

[TABLE]

and, therefore,

[TABLE]

Hence, if (31) holds true then (33) is verified.

To prove (31) and (32) we proceed as follows. Assume by contradiction that there exist $\hat{t}\in[0,T)$ and $\sigma>0$ with $\hat{t}+\sigma\leq T$ such that

[TABLE]

whereas

[TABLE]

Under assumptions (35) and (36), the result (19) on $n_{M}(t,x)$ allows one to conclude that

[TABLE]

We take $\displaystyle{\hat{x}\in\underset{x\in[0,1]}{\operatorname{arg}\,\operatorname{max}}\;R_{1}(\hat{t},x)}$. Using (37) we find that under assumptions (35) and (36)

[TABLE]

that is,

[TABLE]

Due to assumption (2) this yields $R_{1}(\hat{t}+\sigma,\hat{x})>0$, which contradicts the result (17) on $R_{1}$. Hence,

[TABLE]

as well. Therefore, the results (31) and (32) are verified. Hence, both the result (12) and, under the assumption (6) on $n_{1}^{0}$, the result (14) on the limit $n_{1}$ hold.

Part 5 – Step 4. To prove the result (14) on the limit $n_{2}$ we proceed as follows. The weak convergence result (9) for $n_{1\varepsilon}$ along with the result (14) on the limit $n_{1}$ and the fact that $Q_{1}(\rho_{M\,\varepsilon}(t),x)>0$ on $[0,T]\times[0,1]$ ensures that for all $\tau\in[0,T)$ and for all $\sigma>0$ with $\tau+\sigma\leq T$ we have

[TABLE]

Moreover, integrating over the set $[\tau,\tau+\sigma]\times[0,1]$ both sides of the IDE (7) for $n_{2\varepsilon}$ we find

[TABLE]

and, since the upper bound (15) for $\rho_{2\varepsilon}$ ensures that

[TABLE]

we conclude that

[TABLE]

This along with (38) implies that

[TABLE]

Coherently with the fact that $P_{2}(\rho_{M\,\varepsilon}(t),x)<0$ on $[0,1]$ for a.e. $t\in(0,T)$ when $\varepsilon\to 0$, the result given by (39) yields

[TABLE]

This along with the calculations carried out in Part 5 – Step 2 allow us to conclude that the result (14) on the limit $n_{2}$ is verified. The results (14) for the limits $n_{i}$ with $i=3,\ldots,M$ can be proved in a similar way through a bootstrap argument.

$\Box$

The general asymptotic result established by Theorem 3.1 put on a rigorous mathematical basis the idea that clonal selection is driven by the self-renewal fraction of leukemic stem cells, as exemplified by Corollaries 3.2-3.4 given in Section 3.2.

3.2 Biological implications of Theorem 3.1

Building on the ideas presented in previous studies (Stiehl et al., 2015, Stiehl and Marciniak-Czochra, 2011, 2012), here we focus on the case where there are $M=3$ possible maturation stages, that is, stem cells ($i=1$), progenitor cells ($i=2$) and mature cells/leukemic blasts ($i=3$).

In this case, Corollary 3.2 of Theorem 3.1 shows that if

[TABLE]

then in the limit $\varepsilon\rightarrow 0$ the population density functions of stem cells $n_{1\varepsilon}(t,x)$, progenitor cells $n_{2\varepsilon}(t,x)$ and mature cells/leukemic blasts $n_{3\varepsilon}(t,x)$ become concentrated as Dirac masses centred at the point $\overline{x}$. Analogously, Corollary 3.3 of Theorem 3.1 shows that if

[TABLE]

then when $\varepsilon\rightarrow 0$ the cell population density functions $n_{1\varepsilon}(t,x)$, $n_{2\varepsilon}(t,x)$ and $n_{3\varepsilon}(t,x)$ become concentrated as weighted sums of Dirac masses centred at the points $\{\overline{x}_{1},\dots,\overline{x}_{N}\}$. In both cases, the centres of the Dirac masses do not depend on the functions $p_{1}(x)$, $a_{2}(x)$ and $p_{2}(x)$.

From a biological point of view, the centres of the Dirac masses can be understood as the leukemic clones that are selected for in the limit of many cell generations. Therefore, the asymptotic results of Corollary 3.2 and Corollary 3.3 formalise the idea that clonal selection is controlled by the self-renewal fraction of leukemic stem cells, as the leukemic clone(s) with the highest stem cell self-renewal fraction are ultimately selected.

The results of Corollary 3.2 and Corollary 3.3 are complemented by Corollary 3.4 of Theorem 3.1, which shows that if the function $a_{1}(x)$ is constant, i.e. if

[TABLE]

then in the asymptotic regime $\varepsilon\rightarrow 0$ the cell population density functions $n_{1\varepsilon}(t,x)$, $n_{2\varepsilon}(t,x)$ and $n_{3\varepsilon}(t,x)$ do not become concentrated as Dirac masses. Biologically, this indicates that if the stem cell self-renewal fraction is the same for all leukemic clones, then clonal selection will not occur and no specific clones will be selected.

Corollary 3.2** (Selection of one single clone)**

Assume $M=3$ and let the assumptions of Theorem 3.1 along with the additional assumption (41) hold. Then, the measures $n_{1}$, $n_{2}$ and $n_{3}$ given by Theorem 3.1 are such that

[TABLE]

Proof

For $M=3$, under the additional assumption (41), the result given by (44) is a straightforward consequence of the general asymptotic result (14).

$\Box$

Corollary 3.3** (Selection of multiple clones)**

Assume $M=3$ and let the assumptions of Theorem 3.1 along with the additional assumption (42) hold. Then, the measures $n_{1}$, $n_{2}$ and $n_{3}$ given by Theorem 3.1 are such that

[TABLE]

Proof

For $M=3$, under the additional assumption (42), the result given by (45) is a straightforward consequence of the general asymptotic result (14).

$\Box$

Corollary 3.4** (Absence of clonal selection)**

Assume $M=3$ and let the assumptions of Theorem 3.1 along with the additional assumption (43) hold. Moreover, let $M=3$. Then, the measures $n_{1}$, $n_{2}$ and $n_{3}$ given by Theorem 3.1 are such that

[TABLE]

Proof

For $M=3$, under the additional assumption (43), the result given by (46) is a straightforward consequence of the general asymptotic result (14).

$\Box$

4 Numerical solutions

We integrate the asymptotic results established by Corollaries 3.2 - 3.4 of Theorem 3.1 with numerical solutions of the Cauchy problem defined by the system of IDEs (5) with $M=3$ complemented with biologically relevant initial conditions (vid. Section 4.2). We parameterise the model based on patient data from the existing literature (vid. Section 4.1).

4.1 Setup of numerical simulations and model calibration

To construct numerical solutions, we choose a uniform discretisation of the interval $[0,1]$ that consists of 1000 points. We assume $t\in[0,T]$ and we approximate the IDE system for the population density functions $n_{1}(t,x)$, $n_{2}(t,x)$ and $n_{3}(t,x)$ using the forward Euler method with step size $10^{-4}$. We choose $T=10^{4}$. Under the parameter settings considered here, such a value of $T$ corresponds to a time frame of approximatively 30 years after the appearance of the first leukemic clones. This is biologically reasonable since acute leukemia has a pronounced peak of incidence in late adulthood.

Numerical computations are performed in Matlab for two main parameter settings: one under which the total cell density functions $\rho_{1}(t)$, $\rho_{2}(t)$ and $\rho_{3}(t)$ converge to some stable values (vid. Section 4.1.1), and the other such that the total cell density functions undergo oscillations (vid. Section 4.1.2).

4.1.1 Model calibration 1

To model an initial scenario whereby, due to clonal heterogeneity, all possible leukemic clones are present in small numbers in every maturation compartment, and the total densities of healthy cells are close to the cell counts at physiological equilibrium, we use the following initial data

[TABLE]

with

[TABLE]

Such initial conditions satisfy assumptions (6).

Multi-compartment models of hematopoiesis after bone marrow transplantation allow to estimate the proliferation rates and the self-renewal fractions of healthy stem and progenitor cells, as well as the parameters of the feedback signal (Stiehl et al., 2015, 2014c). In particular, coherently with the estimations performed in (Stiehl et al., 2014c), we assume that

[TABLE]

and

[TABLE]

Clearance rates of mature cells and leukemic blasts can be estimated based on patient data and they appear to be between $0.1/day$ and $2.3/day$ (Cartwright et al., 1964, Malinowska et al., 2002, Savitskiy et al., 2003). For this reason, we choose

[TABLE]

In agreement with the estimations performed in (Stiehl et al., 2015, 2014c), which are based on clinical data of blast dynamics in relapsing patients, we impose the following conditions

[TABLE]

The values of the function $a_{2}(x)$ are constrained by assumption (4). Moreover, since it is well accepted that stem cells divide at lower rates compared to progenitor cells (Adams et al., 2015, Cronkite, 1979, Shepherd et al., 2004), we impose the following additional conditions

[TABLE]

We let the continuous structuring variable $x$ parameterise the cell proliferation rate and we use of the following definitions

[TABLE]

In order to fulfil conditions (49), we choose

[TABLE]

Moreover, we choose

[TABLE]

so that conditions (52) and (53) are satisfied. To take into account the complex relationship between proliferation and self-renewal observed in leukemic cells (Doulatov et al., 2009, Ghosh et al., 2016, Kikushige et al., 2015, Lagunas-Rangel and Chávez-Valencia, 2017, Wang et al., 2010, Yassin et al., 2010), we assume the correspondence between the self-renewal fraction and the cell proliferation rate to be non-bijective and among all possible definitions which satisfy conditions (4), (49) and (52) we choose

[TABLE]

and

[TABLE]

or

[TABLE]

with

[TABLE]

or

[TABLE]

Definition (57) models a biological scenario where the leukemic clone identified by $x=0.4$ has the highest self-renewal fraction among progenitor cells. Moreover, definitions (58), (59)-(60), and (61) model, respectively, three distinct situations whereby, due to differential gene expression at different maturation stages:

the leukemic clone identified by $x=0.6$ has the highest self-renewal fraction among stem cells;

the leukemic clones corresponding to $x=0.35$, $x=0.55$, $x=0.7$ and $x=0.85$ have the highest self-renewal fraction among stem cells;

all leukemic clones in the stem-cell compartment have the same self-renewal fraction.

4.1.2 Model calibration 2

On the basis of considerations analogous to those presented in Section 4.1.1, we use the initial data (47) with

[TABLE]

and we choose

[TABLE]

Moreover, we define the functions $p_{1}(x)$ and $p_{2}(x)$ according to (54) with

[TABLE]

and

[TABLE]

Finally, we assume

[TABLE]

and

[TABLE]

i.e. we consider a biological scenario whereby the leukemic clones identified by $x=0.6$ and $x=0.4$ have the highest self-renewal fraction among stem cells and progenitor cells, respectively. The parameters and functions (64)-(67) are such that the values of the functions $p_{1}(x)$, $p_{2}(x)$, $a_{1}(x)$ and $a_{2}(x)$ at the point $x=0.6$ [i.e. the maximum point of the function $a_{1}(x)$] coincide with the parameter values for which the solutions of the ODE system (68) with $M=3$ are known to undergo periodic oscillations (Knauer, 2012, Knauer et al., 2019).

4.2 Main results

In agreement with the asymptotic results established by Corollary 3.2 and Corollary 3.3 of Theorem 3.1, the numerical solutions presented in Fig. 2 and Fig. 3 show that, under the parameter setting given in Section 4.1.1, the population density functions of stem cells $n_{1}(t,x)$, progenitor cells $n_{2}(t,x)$ and mature cells/leukemic blasts $n_{3}(t,x)$ become progressively concentrated at the maximum point(s) of the function $a_{1}(x)$. Moreover, the plots in the insets of Fig. 2 and Fig. 3 highlight the existence of leukemic clones which grow transiently before becoming ultimately extinct. Finally, the numerical solutions displayed in Fig. 4 illustrate how, in agreement with the asymptotic results of Corollary 3.4 of Theorem 3.1, if the function $a_{1}(x)$ is constant, the long-term limits of the population density functions $n_{1}(t,x)$, $n_{2}(t,x)$ and $n_{3}(t,x)$ are not concentrated at any particular point.

As mentioned earlier in the paper, these results communicate the biological notion that the self-renewal fraction of leukemic stem cells determines the fate of clonal selection. In fact, the leukemic clones characterised by the highest stem cell self-renewal fraction are selected, regardless of their properties in terms of stem cell proliferation rate, progenitor cell proliferation rate and progenitor cell self-renewal fraction. This supports the idea that interclonal variability in the self-renewal fraction of leukemic stem cells provides the necessary substrate for clonal selection to act upon.

Under the parameter setting given in Section 4.1.1, the total density of stem cells $\rho_{1}(t)$, progenitor cells $\rho_{2}(t)$ and mature cells/leukemic blasts $\rho_{3}(t)$ converge to some stable values (vid. Fig. 5). However, it is known that for given parameter choices the solutions of the ODE system (68) with $M=3$ undergo periodic oscillations, which result from the occurrence of Hopf bifurcation (Knauer, 2012, Knauer et al., 2019). The numerical results presented in Fig. 7 show that, under the parameter setting given in Section 4.1.2, oscillations emerge in the integrals of the solutions to the IDE system (5) with $M=3$ [i.e. in the total cell densities $\rho_{1}(t)$, $\rho_{2}(t)$ and $\rho_{3}(t)$]. In analogy with the previous cases, the numerical solutions presented in Fig. 6 indicate that the cell population density functions $n_{1}(t,x)$, $n_{2}(t,x)$ and $n_{3}(t,x)$ become progressively concentrated at the maximum point of the function $a_{1}(x)$ – i.e. the leukemic clone characterised by the highest stem cell self-renewal fraction is ultimately selected.

5 Discussion and conclusions

Recent progress in genetic techniques has shed light on the complex co-evolution of malignant cell clones in leukemias (Anderson et al., 2011, Belderbos et al., 2017, Ding et al., 2012, Ley et al., 2008). However, several aspects of clonal selection still remain unclear. In this work, we have adopted a mathematical modelling approach to study clonal selection in acute leukemias, with the aim of supporting a better understanding of the biological mechanisms which underpin observable clonal dynamics.

Our model consists of a system of coupled IDEs that describe the dynamics of cells at different maturation stages – from stem cells through increasingly mature progenitor cells to mature-cells/blasts – seen as distinct compartments. Each maturation compartment is structured by a continuous variable that identifies the clone of the cells. In order to incorporate into our model the high degree of interclonal heterogeneity which is observed in leukemia patients, we let the cellular proliferation rate and self-renewal fraction (i.e. the fraction of progeny cells adopting the same fate as their parent cell) in each maturation compartment be functions of the structuring variable.

In the framework of this model, we have established a number of analytical results which give answers to the open questions Q1-Q3 posed in the introduction of this paper. In summary:

• A1

Clonal selection is driven by the self-renewal fraction of leukemic stem cells. Theorem 3.1 rigorously shows that all leukemic clones with a non-maximal stem cell self-renewal fraction ultimately become extinct, independently of their proliferation rate.

• A2

Non-stem cell properties do not have a substantial impact on clonal selection. The result established by Corollary 3.4 formalises the idea that, in a scenario where the stem cells of all clones have the same self-renewal fraction, one should expect the stable coexistence of all clones to occur.

• A3

Only the clones whose stem cells are endowed with the highest self-renewal fraction can stably coexist in the presence of interclonal heterogeneity. Corollaries 3.2 and 3.3 put on a rigorous basis the notion that the leukemic clones whose stem cells have the highest self-renewal fraction are ultimately selected.

We have integrated our asymptotic results with numerical solutions of a calibrated version of the model based on real patient data. In agreement with the theoretical results by Stiehl et al. (2014a), our numerical results reveal the existence of leukemic clones which display transient growth before becoming extinct. Such an emergent behaviour was not captured by the IDE model considered by Busse et al. (2016), who studied the case where cells of different clones at the same maturation stage have the same proliferation rate. This suggests that interclonal heterogeneity in the cell proliferation rate may have an impact on transient clonal dynamics.

It has been shown using genetic techniques that, in most leukemia patients, the majority of leukemic cells is derived from a relatively small number of clones (Anderson et al., 2011, Ding et al., 2012). Our results indicate that this may be due to the fact that only a few leukemic clones are characterised by a high stem cell self-renewal fraction.

There is accumulating experimental evidence that many of the different genetic mutations involved in the development of leukemia increase the cell self-renewal. Possible examples include the TIM-3/Gal-9 autocrine stimulation (Kikushige et al., 2015) or alterations of Wnt/$\beta$-Catenin signalling (Wang et al., 2010). Moreover, it has been shown that genetic mutations can simultaneously affect the self-renewal and proliferation of cells due to crosstalk between different signalling pathways. For instance, the NUP98-Ddx10 oncogene increases both cell proliferation and self-renewal (Yassin et al., 2010); hyperactivation of the mTOR pathway leads to S6K1-mediated increase in self-renewal and reduction in proliferation (Ghosh et al., 2016); up-regulation of PLZF brings about increased self-renewal and reduced proliferation (Doulatov et al., 2009). The outcome of our analysis, which ascribes a pivotal role to increased self-renewal in orchestrating clonal selection, is in line with the observation that all the aforementioned genetic alterations lead to increased self-renewal, whereas they have divergent effects on cell proliferation.

The result that clonal selection is not influenced by progenitor cell properties is biologically meaningful, since it implies that mutations which affect only the properties of progenitor cells without altering the properties of stem cells cannot lead to the selection of leukemic clones. This result is new and not self-evident, as progenitor cells can expand independently of the influx from the stem cell compartment – although they possess smaller self-renewal fractions than stem cells (Roelofs et al., 2003, Stiehl et al., 2014c, Yamamoto et al., 2013). This insight is also clinically relevant. In fact, although it is known that progenitor cells are more sensitive to treatment interventions than stem cells, which are located in protective niches, our results suggest that manipulations of progenitor cells have no impact on clonal selection phenomena.

A rigorous mathematical understanding of clonal selection is needed since the potentially nonlinear interplay between different genetic and epigenetic hits (Bacher et al., 2008, Gale et al., 2008, Heuser et al., 2009) leads to complex fitness landscapes and non-trivial interdependencies between self-renewal fraction and proliferation rate of leukemic cells. In particular, there is evidence that combinations of leukemic mutations occur frequently in patients (Naoe and Kiyoi, 2013, Pui et al., 2015). An in silico approach can help to disentangle the impact of different mutations on the properties of leukemic stem cells. In this regard, our modelling approach can be further developed in several directions. For instance, in line with what was done for the ODE counterpart of the model presented here (Stiehl et al., 2014a, 2016, 2018), we plan to extend our model to take into account the effect of multiple feedback mechanisms and incorporate the occurrence of de novo mutations. Moreover, along the lines of the modelling method proposed by Doumic et al. (Doumic-Jauffret et al., 2011, Gwiazda et al., 2012), it may be interesting to replace the discrete maturation structure considered in this work by a continuous age structure, which would lead to the definition of a fully-continuously structured population model of clonal selection in acute leukemias. From a mathematical point of view, we also plan to carry out a systematic investigation of the conditions for the emergence of oscillations in the total cell densities, as shown by our numerical solutions.

Appendix A Multi-compartmental ODE model

In a number of previous papers (Stiehl et al., 2014a, 2015, 2018, 2016, Stiehl and Marciniak-Czochra, 2012), it was shown that mathematical models defined in the framework of the following ODE system can effectively recapitulate clinical data from leukemia patients

[TABLE]

with $i=2,\ldots,M-1$, $j=0,\ldots,J$ and $\displaystyle{Z_{M}(t)=\displaystyle{\sum_{j=0}^{J}}N^{j}_{M}(t)}$.

The index $i$ denotes the cell maturation stage while the index $j$ indicates to which leukemic clone the cells belong. In particular, the stem-cell compartment is labelled by the index $i=1$, the indices $i=2,\ldots,M-1$ correspond to increasingly mature progenitor-cell compartments and the mature-cell/blast compartment is labelled by the index $i=M$. Moreover, the index $j=0$ refers to healthy cells, whereas the different leukemic clones are labelled by the indices $j=1,\ldots,J$.

In the system of ODEs (68), the function $N^{j}_{i}(t)$ models the density of cells of clone $j$ at the maturation stage $i$ at the time instant $t\geq 0$. Cells in the compartment $i=M$ do not divide and are cleared from the system at rate $d>0$, which is assumed to be the same for healthy and leukemic cells (Busse et al., 2016, Malinowska et al., 2002, Savitskiy et al., 2003, Stiehl et al., 2014a, 2018). The parameters $p^{j}_{i}>0$ and $a^{j}_{i}>0$ model, respectively, the proliferation rate and the self-renewal fraction of cells of clone $j$ at the maturation stage $i$. Coherently with biological findings (Kondo et al., 1991, Layton et al., 1989, Shinjo et al., 1997), the terms $a^{j}_{i}$ are multiplied by the factor $\frac{1}{1+KZ_{M}(t)}$ to model the fact that the signal which promotes the self-renewal of dividing cells is absorbed by mature cells and leukemic blasts at a rate that depends on the total density of these cells $Z_{M}(t)$. The parameter $K>0$ is related to the degradation rate of the feedback signal by mature cells and leukemic blasts. This has proved to be a biologically consistent way of modelling the effects of feedback signals that control cell self-renewal (Marciniak-Czochra et al., 2009, Stiehl et al., 2014c, 2018, Stiehl and Marciniak-Czochra, 2011, 2012). In principle, the effects of feedback signals that control cell proliferation could also be included. However, it has been demonstrated that such signals have little impact on the dynamics of the blood system (Marciniak-Czochra et al., 2009, Stiehl et al., 2014b, c).

A version of the ODE model (68) with only one leukemic clone and three maturation stages (i.e. $M=3$ and $J=1$) has been fully analysed by Stiehl and Marciniak-Czochra (2012), while a two compartmental version of the model for healthy hematopoiesis (i.e. $M=2$ and $J=0$) has been analysed by Getto et al. (2013), Nakata et al. (2012), Stiehl and Marciniak-Czochra (2011). Possible applications of this model to clinical data can be found in the works by Stiehl et al. (2014a, 2015), whereas applications to healthy hematopoiesis are provided in the publications by Stiehl et al. (2014b, c). Finally, a version of this model with a continuous differentiation structure has been proposed and studied by Doumic-Jauffret et al. (2011).

Appendix B Proof of the uniform upper bounds (15)

In this appendix we prove the upper bounds (15) through a suitable development of the method of proof presented in (Busse et al., 2016). Using the system of IDEs (7), straightforward calculations give the following ODEs

[TABLE]

[TABLE]

and

[TABLE]

Under assumptions (2)-(4), for all $i=1,\ldots,M-1$ we have

[TABLE]

and

[TABLE]

for all $t\geq 0$ and all $x\in[0,1]$. Hence, estimating from above the right-hand sides of the ODEs (69)-(71) using the above estimates on $P_{i}$ and $Q_{i}$ along with the non-negativity of $n_{i\,\varepsilon}$ and $\rho_{i\,\varepsilon}$ for all $i=1,\ldots,M$ gives the following differential inequalities

[TABLE]

[TABLE]

for $i=2,\ldots,M-2$,

[TABLE]

The differential inequality (72) yields

[TABLE]

for all $t\geq 0$. Substituting the estimate (75) into the differential inequality (73) for $i=2$ we find

[TABLE]

for all $t\geq 0$. In a similar way, substituting the estimate (76) into the differential inequality (73) for $i=3$ gives

[TABLE]

for all $t\geq 0$. Using a bootstrap argument based on the method of proof that we have used for the case $i=3$, one can prove that

[TABLE]

for all $t\geq 0$ and for all $i=4,\ldots,M-2$. Finally, substituting the estimate (78) with $i=M-2$ into the differential inequality (74) we obtain

[TABLE]

for all $t\geq 0$. Combining the estimates (75)-(79) yields

[TABLE]

for all $t\geq 0$, which allows us to conclude that

[TABLE]

for all $t\geq 0$ and all $x\in[0,1]$. The latter inequality ensures that $P_{i}(\rho_{M\varepsilon}(t),x)$ satisfies the following relations for all $t\geq 0$ and each $i=1,\ldots,M-1$

[TABLE]

while

[TABLE]

Integrating over $[0,1]$ both sides of the IDEs (7) for $n_{1\,\varepsilon}$ and estimating from above gives the following differential inequality

[TABLE]

from which, using (B) and (B) with $i=1$, we find that for any $\varepsilon>0$

[TABLE]

that is, the upper bound (15) on $\rho_{1\,\varepsilon}$ is verified.

Furthermore, integrating over $[0,1]$ both sides of the IDEs (7) for $n_{2\,\varepsilon}$ and estimating from above using the fact that

[TABLE]

along with the upper bound (15) on $\rho_{1\,\varepsilon}$ gives

[TABLE]

The above differential inequality along with the estimates (B) and (B) with $i=2$ ensures that there exists $C>0$ such that for any $\varepsilon>0$

[TABLE]

for all $t\geq 0$. Hence, the upper bound (15) on $\rho_{2\,\varepsilon}$ is verified. The upper bounds (15) on $\rho_{i\,\varepsilon}$ for $i=3,\ldots,M-1$ can be proved in a similar way using a bootstrap argument.

Finally, integrating over $[0,1]$ both sides of the IDEs (7) for $n_{M\,\varepsilon}$ and estimating from above using the fact that

[TABLE]

along with the upper bound (15) on $\rho_{M-1\,\varepsilon}$ gives the following differential inequality

[TABLE]

which ensures that for any $\varepsilon>0$ we have

[TABLE]

that is, the upper bound (15) on $\rho_{M\,\varepsilon}$ is verified.

Acknowledgements.

TS and AM-C were supported by research funding from the German Research Foundation DFG (SFB 873; subproject B08). TL gratefully acknowledges support from the Heidelberg Graduate School (HGS).