Mathematical Analysis of a Clonal Evolution Model of Tumour Cell Proliferation

Mathematical models of tumour growth provide insight into the dynamic characteristics of tumour cell populations. Important issues are cell proliferation, cell heterogeneity, and cell differentiation. These issues are currently examined in two hypothesized models of tumour evolution: the cancer stem cell (CSC) model and the clonal evolution (CE) model[27, 40]. Both models have scientific support, as well as therapeutic implications, and in fact, both models may be involved in the development of a tumour. The essential distinction of the two models is the role of self-renewal in specific cells, and the fraction of the total cell populations that these cells comprise. Here self-renewal means the ability of a cell to inherit a specific function through an unlimited number of successive cell generations.

The CSC model hypothesizes that a very small sub-population of tumour stem cells generate the entire tumour cell population [43, 49]. In mathematical treatments, these stem cells have the ability to self-renew indefinitely, and through sequential mutations generate all the heterogeneous and differentiated cell types comprising the tumour [1, 6, 13, 18, 23, 24, 34, 37, 38, 39, 41, 44, 47, 50, 51, 54]. This stem cell population lies at the apex of a hierarchical structure of cell types, and tumour evolution is dependent on their unconstrained self-renewing ability [30].

The CE model hypothesizes that a tumour population is composed of multiple genetically identical clones, which have the possibility of mutation, selection, and expansion [8, 43, 48]. In the CE model all undifferentiated cells have a self-renewing capacity for contributing to the tumour evolution. These cells, however, are not at the apex of a hierarchal tree, but rather dispersed widely throughout the tumour cell population as a large fraction of the total tumour cell count.

A central element of cell self-renewal is the Hayflick limit, which constrains differentiated cell lines to a finite number of divisions [28]. As differentiated cells divide, telomeres (nucleotide sequences at the ends of chromosomes) shorten until a critical limit is reached, and further divisions are prohibited [33, 46]. The existence of a mechanism which reverses telomere shortening was predicted several decades ago. The CSC model hypothesizes that cancer stem cells circumvent telomere shortening by using the enzyme telomerase to replace their telomeres, and thus obtain the ability to divide indefinitely [43]. Recently, it has been shown that around 90% of all types of human cancer exhibit a form of telomerase activation [27]. In the CSC model this property resides in an extremely small sub-population, from which all the differentiated tumour cells derive. In contrast, the CE model hypothesizes that a large number of undifferentiated cells, in clonal sub-populations, possess the ability to restore telomeres, and thus sustain the tumour evolution. These two models differ greatly in their fraction of self-renewing cells within the total population. Mathematical models provide a way to compare these self-renewal properties in proliferating cell populations.

In order to model the dynamics of self-renewing cells lineages, which correspond physiologically to chromosomal telomere lengths, is it necessary to track all cells through successive generations. Many mathematical treatments of telomere structure in cell population dynamics have been developed [4, 5, 7, 16, 31]. The CSC model has been treated for example in [31], where telomere shortening is investigated in continuum differential equation models. The key focus of these treatments is that mother stem cells produce two daughter cells, one of which is a stem cell, and the other, a differentiated cell, with a limited number of future divisions. These models describe a finite number of sub-populations, each with an idealised precise telomere length. The self-renewal property of a stem cell is captured by one daughter cell remaining in the highest telomere class and the other daughter cell transiting to a lower telomere class. A descent of telomere length shortening continues in each class, with one daughter cell retaining the length of the mother and the other a lower length. This model yields a minute fraction of stem cells at the apex of the total population, consistent with the CSC literature [43].

The objective of this paper is to develop a mathematical analysis for the alternative CE model and to quantify its dynamic properties. Our model here incorporates continuum, rather than discrete, telomere lengths in cells. This idealised continuum telomere length is assumed for convenience, to avoid unwieldy compartmentalisation with each telomere length subclass. In our CE model of telomere shortening there is an unbalanced division of a mother cell to two daughter cells in terms of telomere length. The telomere length of a daughter cell may be less (corresponding to a differentiated cell), or may be equal or greater (corresponding to a self-renewing cell) than the mother cell. In this way the restoration of telomeres and the self-renewal property of cells is distributed through a large proportion of the tumour cell population, consistent with the CE literature[43]. The distribution of daughter cell telomere lengths is governed by a mathematically formulated rule that assigns the telomere restoration property to cells which may be viewed as those capable of indefinite divisions in each of the diverse clonal sub-populations.

Our CE model belongs to the class of continuum structured population models, with age and time as dynamic variables, and telomere length as a population structure variable. In the past three decades physiologically structured population models have been increasingly utilised to shed light on some important phenomena of cell populations [4, 5, 15, 16, 29, 31, 36]. The power of structuring a population with respect to physiological variables is of great value in understanding the evolution of biological populations. There is an increasing literature of physiologically structured population models with more than one (physiological) structuring variable. The development of a unified mathematical framework, where structuring variables play substantially different roles, promises to be extremely challenging from the theoretical point of view.

We first consider the following linear model for an age and telomere length structured proliferating cancer cell population.

Above $p(a,l,t)$ stands for the density of cells of age $a$, having telomere length $l$ at time $t$. We assume a maximum cell age denoted by $a_{m}$ and a maximum telomere length denoted by $l_{m}$. The population count at time $t$ of cells with age between $a_{1}$ and $a_{2}$ and telomere length between $l_{1}$ and $l_{2}$ is

and the total population of all cells at time $t$ is

$\mu(a,l)$ quantifies the natural mortality of cells of age $a$ and telomere length $l$. A mother cell of age $a$ and telomere length $l$ divides into two daughter cells of age $0$ having (possibly) different telomere lengths, at a rate determined by the function $\beta(a,l)$. The function $r(l,\hat{l})$ describes the distribution of daughter telomere lengths $l$ from a mother cell of telomere length $\hat{l}$. The boundary condition (1.2) accounts for cell division of a mother cell into two daughter cells and requires that

From a probabilistic interpretation of $r(l,\hat{l})$ it would be natural to normalise the maximal telomere length $l_{m}$ to $1$. Throughout we retain $l_{m}$ as a general parameter. We will impose further regularity assumptions on the model ingredients later on. Note that our model (1.1)-(1.3) can be considered as a continuous telomere length structured counterpart of the model recently introduced in [31].

Our starting model (1.1)-(1.3) is a linear one, moreover the telomere length $l$ only plays an important role in the somewhat unusual boundary condition (1.2). Hence to establish the existence of the governing linear semigroup (and therefore the existence of mild solutions of the PDE (1.1)-(1.3)) we use a boundary perturbation result due to Greiner, see Theorem 2.3 in [26]; see also [25] and [14] for similar general results. It is very natural to apply the boundary perturbation result of Greiner, since the unperturbed generator (arising from equation (1.1) with zero flux boundary condition) is readily shown to generate a translation semigroup. Moreover, it has the added advantage that the spectrally determined growth behaviour of the semigroup follows almost instantaneously, see Proposition 3.1 in [26]. For basic definitions and results not introduced in the section we refer the reader to [3, 19].

To apply the perturbation result of Greiner we set the framework as follows. Assume that $\beta\in C^{1}_{+}([0,a_{m}]\times[0,l_{m}])$ and $\mu,r\in L^{\infty}_{+}((0,a_{m})\times(0,l_{m}))$, and that all of the model parameters are non-negative. In particular, it is natural to assume that $\beta$ does not vanish (in $a$) identically for any $l>0$. We also impose the natural assumption that cells reaching the maximal age do not reproduce anymore, i.e. $\beta(a_{m},\cdot)\equiv 0$. This assumption is completely natural from the biological point of view. If cells would still divide upon reaching the maximum age then it would be natural to extend the age-interval beyond $a_{m}$. In other words, we have set the maximal age $a_{m}$ such that cells of this age do not divide for any more. On the other hand, we will see later when we introduce nonlinearities in model (1.1)-(1.3), that non-reproducing cells still play a role in the population dynamics, in particular they may affect competition induced mortality. Also note that alternatively, we could have assumed that the mortality is locally integrable with respect to age, but $\int_{0}^{a_{m}}\mu(a,\circ)\,\mathrm{d}a=\infty$ holds, i.e. no individuals survive the maximal age $a_{m}$.

For the linear problem (1.1)-(1.3), since we are dealing with density functions, the natural choice of state space is the following Lebesgue space.

Elements of $\mathcal{X}$ above can be understood as equivalence classes of measurable functions $f(\cdot)(\circ)$ on the square $(0,a_{m})\times(0,l_{m})$ such that $\int_{0}^{a_{m}}\left|f(a)(\circ)\right|\,\mathrm{d}a\in L^{1}(0,l_{m})$. We further set $\mathcal{Y}=L^{1}(0,l_{m})$. We define the operators $\mathcal{A}$ and $\mathcal{B}$ as follows

with

We further introduce the norm

With the $||\cdot||_{A}$ norm $D(\mathcal{A})$ is complete, and the maximal operator $\mathcal{A}\,:\,(D(\mathcal{A}),||\cdot||_{A})\to\mathcal{X}$ is continuous and linear. Furthermore, we define

Then $\mathcal{L}$ is also continuous and linear, and we have Im$(\mathcal{L})=\mathcal{Y}$. We denote by $\mathcal{A}_{0}$ the restriction of $\mathcal{A}$ to Ker$(\mathcal{L})$. It is then clear that $\mathcal{A}_{0}$ generates the strongly continuous and nilpotent shift semigroup $\mathcal{S}$, explicitly given as

In particular note that for $t>a_{m}$ we have $(\mathcal{S}(t)\,u_{0})\equiv 0$ for any initial condition $u_{0}\in\mathcal{X}_{+}$. We now define the bounded linear perturbing operator $\Phi\,:\,\mathcal{X}\to\mathcal{Y}$ as follows

We also define the corresponding perturbed generator $\mathcal{A}_{\Phi}$ as

We recall the main result from [26] for the reader’s convenience.

We now apply Theorem 2.1 to establish the existence of the governing linear semigroup. In our setting we have $\mathcal{X}^{*}=L^{\infty}((0,a_{m})\times(0,l_{m}))$, and $\mathcal{Y}^{*}=L^{\infty}(0,l_{m})$. Furthermore, to compute the adjoint $\Phi^{*}\,:\,D(\Phi^{*})\subset\mathcal{Y}^{*}\to\mathcal{X}^{*}$, we note that

if we let

Next we note that (see [32, Sect.III.5]) $g\in D(\mathcal{A}_{0}^{*})$ if there exists an $f\in\mathcal{X}^{*}$ such that

For any $u\in D(\mathcal{A}_{0})$ integration by parts yields

for $g\in\left\{W^{1,\infty}\left((0,a_{m});L^{1}(0,l_{m})\right)\,|\,g(a_{m},\cdot)\equiv 0\right\}$, and if we let $f=\frac{\partial g}{\partial a}$. Hence it follows from the regularity assumptions on $\beta$ and $r$ that we have

Hence Theorem 2.1 implies that $\mathcal{A}_{\Phi}$ generates a strongly continuous semigroup.

Next we note that for $\lambda\in\rho(\mathcal{A}_{0})$, the operator $\mathcal{L}_{|\,Ker(\lambda-\mathcal{A})}$ is a continuous bijection from $\left(ker(\lambda-\mathcal{A}),||\cdot||_{A}\right)$ onto $\mathcal{Y}$, hence its inverse

is continuous, for $\lambda\in\rho(\mathcal{A}_{0})$. In particular, a straightforward calculation shows that in our setting we have $\mathcal{L}_{\lambda}\,:\,y(\circ)\to e^{-\lambda\,\cdot}y(\circ)$, and therefore for $\lambda$ large enough $(\mathcal{I}-\mathcal{L}_{\lambda}\,\Phi)$ is invertible and positive. Hence by Lemma 1.4 in [26] we have that $R(\lambda,\mathcal{A}_{\Phi})$ is positive for $\lambda$ large enough. Since $\mathcal{B}$ is a bounded multiplication operator, we draw the following conclusion.

In this section we study the asymptotic behaviour of solutions of the linear model (1.1)-(1.3). In particular, as we will see, we are going to characterise the spectral bound of the linear semigroup generator $\mathcal{A}_{\Phi}+\mathcal{B}$ implicitly via the spectral radius of an associated (bounded) integral operator. We will then obtain estimates for the spectral radius of this integral operator.

As we have pointed out earlier one of the advantages of the perturbation theorem of Greiner (Theorem 2.1) is that it allows us to establish a desirable regularity property of the semigroup in a straightforward fashion. To this end, for the rest of the paper we further assume that $r$ satisfies the following regularity condition.

Note that, for example if $r$ is continuous on the square $[0,l_{m}]\times[0,l_{m}]$, then condition (3.14) clearly holds. We recall now from [26] the result we are going to apply, for the readers convenience. In particular, Proposition 3.1 in [26] reads as follows.

In the proposition above $\rho_{+}(\mathcal{A}_{0})$ stands for the component of the resolvent set of $\mathcal{A}_{0}$, which is unbounded to the right. We apply now Proposition 3.3 in our setting. In what follows, with a slight abuse of notation, we will denote the semigroup generated by $\mathcal{A}_{\Phi}+\mathcal{B}$ by $\mathcal{T}_{\Phi}$.

We shall point out that we really needed to utilise Proposition 3.3 by Greiner, to obtain that the asymptotic behaviour of the semigroup is determined by the leading eigenvalue of its generator (if it exists). This is due to the distributed structuring with respect to the telomere length of cells, which implies that the governing semigroup is not necessarily eventually differentiable; hence we could not conclude, as for example in [22] for a classic size-structured models, that the semigroup is eventually compact, and henceforth the spectral mapping theorem holds true.

As we have seen earlier the semigroup $\mathcal{T}_{\Phi}$ generated by $\mathcal{A}_{\Phi}+\mathcal{B}$ is clearly positive, but it may not necessarily be irreducible. To see this, recall from [19], that the semigroup generated by $\mathcal{A}_{\Phi}+\mathcal{B}$ is irreducible if and only if for every $f$, $0\not\equiv f\in\mathcal{X}_{+}$, we have that $R(\lambda,\mathcal{A}_{\Phi}+\mathcal{B})f\gg 0$, i.e. the resolvent is strictly positive, for some $\lambda>s(\mathcal{A}_{\Phi}+\mathcal{B})$.

Let $f\in\mathcal{X}_{+}$, and note that the solution of the resolvent equation
$R(\lambda,\mathcal{A}_{\Phi}+\mathcal{B})f=u$ is

Applying the boundary operator $\Phi$ on both sides of equation (3.16), it is easily shown that $u(0,\cdot)$ satisfies the inhomogeneous integral equation

Above in (3.17) we introduced the notation

Note that since telomere length is preserved during the lifetime of an individual, it is intuitively clear that the semigroup is irreducible if offspring of all telomere length is produced by some individuals. In other words, if there were individuals of particular telomere lengths who would not produce individuals of any other lengths, then the semigroup would be reducible. We formulate now a rigorous condition for the irreducibility of the semigroup, as this will play an important role later in the qualitative analysis of model (1.1)-(1.3).

Proposition 3.4 implies that the asymptotic behaviour of solutions of model (1.1)-(1.3) is determined by the eigenvalues of the generator (if there are any), hence we study this eigenvalue problem now. In particular, the solution of the eigenvalue problem

is given by

Applying the boundary operator $\Phi$ on both sides of the equation above we have

where we defined

Hence $\lambda\in\mathbb{C}$ is an eigenvalue of $\mathcal{A}_{\Phi}+\mathcal{B}$ if and only if, for the given $\lambda$, the integral equation (3.21) has a non-trivial solution $\psi(0,\cdot)$. Then, the eigenvector corresponding to $\lambda$ is given by (3.20). We are in particular interested in the leading eigenvalue (if it exists), which is the spectral bound of the generator, since the semigroup $\mathcal{T}_{\Phi}$ is positive. This dominant real eigenvalue, together with the corresponding eigenspace, determines the asymptotic behaviour of solutions. Also note that equation (3.21) gives naturally rise to define a parametrised family of (positive and bounded) integral operators $\mathcal{O}_{\lambda}$ as

With this, the characteristic equation (3.21), which is notably a functional equation, in contrast to a scalar equation in case of a model with age-structure only; can be viewed as an eigenvalue problem for a bounded linear integral operator. More precisely, $\lambda$ is an eigenvalue of the generator $\mathcal{A}_{\Phi}+\mathcal{B}$, if and only the integral operator $\mathcal{O}_{\lambda}$ has eigenvalue $1$. Note that, since $K$ defined in (3.22) is strictly positive, the integral operator $\mathcal{O}_{\lambda}$ (for every $\lambda\in\mathbb{R}$) is irreducible if and only if condition (3.18) holds [45, Ch.V]. Hence, rightly so, the irreducibility conditions of the semigroup $\mathcal{T}_{\Phi}$ and the integral operator $\mathcal{O}_{\lambda}$ coincide.

Also note that the function $[0,\infty)\ni\lambda\to r(\mathcal{O}_{\lambda})$ is continuous and strictly monotone decreasing. These properties can be established by using perturbation results from [2] and [32], respectively; see also [9, 11] for similar developments. Also note that if $r$ satisfies condition (3.14) then $\mathcal{O}_{0}$ is shown to be compact exactly in the same way as the operator $\Phi$ was shown to be compact earlier. Hence the spectrum of $\mathcal{O}_{0}$ may contain only eigenvalues and $0$.

Therefore, in the case when the spectrum of $\mathcal{A}_{\Phi}+\mathcal{B}$ is not empty, we have the following complete characterisation of the asymptotic behaviour of solutions of model (1.1)-(1.3).

1. (1)

   If $r(\mathcal{O}_{0})<1$, then solutions of model (1.1)-(1.3) decay exponentially.

2. (2)

   If $r(\mathcal{O}_{0})>1$, then solutions of model (1.1)-(1.3) grow exponentially. Moreover, if $r$ satisfies condition (3.18), then solutions exhibit asynchronous exponential growth.

3. (3)

   If $r(\mathcal{O}_{0})=1$, then for any eigenvector $\psi(0,\cdot)$ corresponding to the spectral radius $1$ of $\mathcal{O}_{0}$, the function $\psi$ in (3.20) determines a one-parameter family of positive steady states of model (1.1)-(1.3). If $r$ satisfies (3.18), then there is only one such family of steady states, moreover they are strictly positive, too.

After the previous general analysis of the asymptotic behaviour of our model next we intend to study the effect of the well-known capacity of telomere restoring of cancer cells on the dynamics. In particular we are going to show, that at least for some general classes of the model ingredients, the telomere length restoring capacity of cancer cells may have a drastic effect on the asymptotic behaviour of solutions already in the linear model (1.1)-(1.3). First we note that in the absence of telomere restoring capacity of cells, the function $r$ necessarily vanishes on the half square $\hat{l}\leq l$. That is, when a mother cell divides, it only gives birth to daughter cells with shorter telomeres. This is a well-known mechanism observed in healthy cell populations. In particular this telomere shortening of healthy cells results in apoptosis (when reaching the celebrated Hayflick limit) and prevent the possibility of drastic mutations caused by a very large number of iterations of faulty DNA replication. In this case the boundary condition (1.2) can be rewritten as

This in particular implies that the eigenvalue problem (3.21) now reads

where we also introduced the notation $\Psi(l):=\psi(0,l)$, for simplicity. Let us assume now that the the function $r$ is separable, i.e. $r(l,\hat{l})=r_{1}(l)r_{2}(\hat{l})$ holds for some functions $r_{1},r_{2}$. For example we may assume that $r_{1}$ is continuously differentiable and $r_{2}$ is bounded, and that $r_{1}$ is positive while $r_{2}$ is non-negative. Then assumption (3.14) clearly holds. In this case differentiating equation (3.25) (assuming that an eigenvector $\psi$ with a smooth $\Psi$ component exists) yields the differential equation

together with the initial condition

The solution of (3.26) is

which, utilising (3.27), leads to the following characteristic equation

It is clear that equation (3.29) does not admit any solution $\lambda\in\mathbb{R}$, which, together with the positivity of the semigroup, implies that the spectrum of $\mathcal{A}_{\Phi}+\mathcal{B}$ does not contain any eigenvalue with a corresponding eigenvector $\psi$ with a continuously differentiable $\psi(0,l)$. On the other hand it is clear from equation (3.25) that any eigenvector is continuous with respect to its second variable. From the biological point of view this phenomenon is associated with the constant loss of telomere length of newborn cells. Indeed, in the absence of telomere restoring capacity we may expect that the cell population accumulates at the minimal length. From the mathematical point of view, the non-existence of differentiable eigenvectors is associated with the fact that the governing semigroup cannot shown to be eventually differentiable due to the telomere length structuring.

Next we consider what happens if we account for the telomere restoring capacity of cancer cells. In particular we are going to show that in this case even exponential growth of the cancer cell population is possible. As before, we start with the case of a separable $r$, i.e. we assume that $r(l,\hat{l})=r_{1}(l)r_{2}(\hat{l})$ for some functions $r_{1},r_{2}$. In this case the integral operator $\mathcal{O}_{0}$ is of rank one, and it is rather straightforward to exactly determine the spectral radius of $\mathcal{O}_{0}$, which, as we have seen earlier, determines the asymptotic behaviour of the semigroup $\mathcal{T}_{\Phi}$. In particular, we have

From (3.30) we can see that depending on the functions $r_{1},r_{2},\mu$ and $\beta$, the spectral radius $r(\mathcal{O}_{0})$ may be greater than $1$. For example in the case of constant functions $r_{1},r_{2},\mu,\beta$, and setting $l_{m}=a_{m}=1$ (note that we can always normalise the maximal age and telomere length), we have

To obtain estimates for the spectral radius of $\mathcal{O}_{0}$ in the more general and difficult non-separable case, note that the Krein-Rutman theorem asserts that if $\mathcal{O}_{0}$ is compact and positive, and its spectral radius is positive, then it has a positive (not necessarily strictly positive) eigenvector corresponding to its spectral radius. On the other hand de Pagter proved in [12] that if the operator is also irreducible then its spectral radius is strictly positive. As we noted earlier if $r$ satisfies condition (3.14) then $\mathcal{O}_{0}$ is shown to be compact in the same way as the operator $\Phi$. Assuming now that the spectral radius is positive, let $u\in L^{1}_{+}(0,l_{m})$ denote the positive eigenvector corresponding to the spectral radius. Then we have

which yields

This observation allows us to obtain immediately the following estimates for the spectral radius

To obtain different estimates, in particular when $\mathcal{O}_{0}$ may not be irreducible we are going to utilise some minimax principles established in [35]. Recall that if $\mathcal{X}$ is a Banach lattice with positive cone $K$, and with dual space $\mathcal{X}^{*}$ and dual cone $K^{*}$, respectively; then a set $H^{\prime}\subseteq K^{*}$ is called $K$-total if and only if from $\langle x,x^{\prime}\rangle\geq 0,\,\forall\,x^{\prime}\in H^{\prime}$ it follows that $x\in K$. Then for any positive linear endomorphism $\mathcal{O}$ on a Banach lattice $\mathcal{X}$ and for any $x\in K$ one defines

Recall that by Lemma 3.1 in [35] for any $K$-total set $H^{\prime}\subseteq K^{*}$ we have

Moreover, Lemma 3.3 in [35] asserts that for any $0\not\equiv x\in K$ we have $r_{x}(\mathcal{O})\leq r(\mathcal{O})$. If we let $x\equiv 1$, then we have for any $x^{\prime}\in K^{*}$

Similarly, recall from [35] that if we define

then for every $x\in K$ we have $r(\mathcal{O})\leq r^{x}(\mathcal{O})$. Again, choosing $x\equiv 1$, we have for any $x^{\prime}\in K^{*}$

Hence we obtain the following estimates for the spectral radius of $\mathcal{O}_{0}$

Note that the estimates (3.33) and (3.38) are quite different, in general.

We next provide hypotheses on $r(l,\hat{l})$ that yield specific growth behavior of the solutions in the presence or absence of highest telomere class renewal. It is known that in the discrete telomere length case, the cell population can have polynomial growth or decay with cells with shortest telomere length having the highest power growth over time (see e.g. [4],[5],[16]). Similar results hold in the continuum telomere length case if we divide the population into telomere length classes. We first consider the case of no self-renewal within any class, that is, all cell divisions result in daughter cells in a shorter telomere length class.