Existence, continuous dependence and finite time stability for Riemann-Liouville fractional differential equations with a constant delay

1.   Introduction

The notion that cancer is an evolutionary process, driven by random mutations and selection of the fittest clones, has become so widespread to be, nowadays, commonplace. Indeed random mutations are the source of genetic diversity of clonal population which has a great relevance to explain the emergence of cancer types resistant to therapies. On the other hand the growth of tumor species depends on how effectively they are able to access resources and, on the other hand, on how successfully they develop mechanisms to prevent detection and elimination by the immune system ^[17]. Taking a step further one can say that cancer is an ecosystem ^[19,27] formed by coexisting populations, embedded in an environment comprising normal and immune cells ^[23].

The multiplicity of species in cancer populations has a clear relevance for the design of therapies as heterogeneity is a major factor in cancer drug resistance, see e.g. ^[29]; even though a therapy can decimate a cancer type, one or more variants of the tumor population exist which are resistant, driving to the resurgence of treatment-refractory disease ^[16]. This observation has led to the idea of combination therapy, in which agents with different actions are combined, thus increasing the likelihood of synergistic anti-tumor effects, ^[6,10,29].

A more radical approach maintains that it is impossible to eliminate all cancer species, of which one will be anyway resistant, but it is more appropriate to devise therapies to keep cancer cells numbers under control: this is the so called adaptive therapy. Not surprisingly, given the relevance of the problem, there exist a large literature on mathematical models of cancer dynamics. For a review see for instance ^[2,11,32] and references therein. In particular several contributions can be found in the framework of population dynamics: among others on tumor immune interaction, several populations are considered in ^[9,33] and spatial-temporal dynamics in ^[1]. The effects of therapies are studied in the context of evolutionary dynamics ^[14,15], while the immuno-therapy are considered in ^[8,13].

In most models it is assumed that cancers grows by cell-autonomous proliferation, manifested as an initial exponential increase in cell number, eventually limited by the carrying capacity, see e.g. ^[3]; in the following, models of this type will be called exponential. However, there is a increasing body of works showing that, in some cases, cancer growth at low densities deviates from the exponential growth and that cell population kinetics are best described by considering processes models involving the Allee effect, see for instance ^[20,22,26], a biological phenomenon characterized by a positive correlation between population size or density and the mean individual fitness of a population or species.

In particular, two types of Allee effect can occur: weak or strong. In the weak effect the rate of growth increases with the population density and remains always positive; the strong effect, in addition, is characterized by a population threshold below which the population decreases until it reaches extinction, ^[7]. Models incorporating the Allee effect have been developed to fit data from clinical outcomes following tumor resection ^[26] and the in vitro BT-474 breast cancer cell line data, ^[20]. In addition, observations of in vitro cell growth have shown that low seeding density may have an adverse effect on cancer growth, ^[22].

Allee effect can originate from a variety of mechanisms including cooperative interactions which can potentially be relevant to cancer ^[22]. Cooperation is common to a wide variety of organisms, see for instance ^[7,25]; in cancer it might be needed to produce a sufficient density of growth factors required for tumour proliferation or pro-angiogenic growth factors such as vascular endothelial growth factor A (VEGFA), ^[22].

Some specific mathematical models of mechanisms leading to an Allee effect have been presented in ^[4,21] and the role of the Allee effect in the spreading of cancer cells has been investigated in ^[12,30].

All these findings are very relevant both to understand how cancer evolves and to determine the more suitable type of therapy. In particular, as mentioned before, it has been suggested recently that the maximum tolerated dose therapy (MTD), commonly used, should be replaced by an adaptive therapy, whose aim is to control, rather than to eradicate, cancer. The rationale of is that MTD therapy is not able, in general, to completely destroy cancer populations, and in some cases attempts to eradication can lead to worse outcomes, ^[15]. Simulations and calculation to support this idea are based on different variations of the exponential model. However, in situations where the Allee effect occurs, MTD could be a better option, in that it can eliminate permanently cancer populations, if the effect is strong or, at least, delay significantly cancer regrowth, as in case of weak effect.

Motivated by these observations, in this work it will be investigated how the Allee effect can influence the outcome of different therapies. In details, the paper is organized as follows. In the next section, population dynamics under Allee effect will be analyzed, with special reference to the times required to reach a given concentrations. These results will be extended, in Section 3, in case when different cancer populations are present. In Section 4, some simulations will be shown of the effect of different therapies if the Allee effect is present or otherwise.

2.   Single cancer clone

This section will study the case of a single cancer species or clone, pointing out the differences between the usual growth equations and those incorporating the Allee effect (weak and strong). Even if the treatment is quite straightforward, it can be useful as introduction to more complex cases.

There are different ways to model the weak and strong Allee effects, see e.g. ^[7]: for instance denoting by $x$ the number of elements in the population (cancer cells in our case), the evolution of the population con be written as:

where $r$ is the reproduction rate and $K$ corresponds to the carrying capacity, i.e. the maximum value that $x$ can take. Parameter $A$ is a threshold: $A > 0$ corresponds to the strong effect whereas $A \leq 0$ models the weak case. In the following, it will be assumed $A = 0$ to model the weak Allee effect, so that Eq (2.1) becomes

The parameter $c$ is related to the intra-specific cooperative interactions in a population. More precisely $c$ measures the amount of cooperation that a population needs to grow: if $c$ is very large an appreciable growth occurs only for high population density (compare ^[5]): on the contrary if $c = 0$ species can increase even without cooperation and Eq (2.2) reduces to

describing an initial exponential growth eventually limited by a logistic term; in the sequel this process will be called exponential, as usually done in the literature. It is obvious that Eqs. (2.3) and (2.2) have the same stationary states, $x_s = 0$ (unstable), and $x_s = K$ (stable); thus logistic and weak Allee processes differ mainly on the rate of growth and this difference is large when $x$ is small. In order to elucidate this point, let's suppose $x < < K$ so that the logistic term $x/K$ can be disregarded: then Eq (2.2) can be easily integrated and the results is

where $x_0$ is the initial value and $x > x_0$. On the other hand, the integration of Eq (2.3), again neglecting the term $x/K$, leads to

Thus, to attain a given value $x^* < < K$ (which could be, for instance a detectability threshold) the weak Allee process requires, with respect to the logistic case, an additional time

which increases with $c$, as expected. In Figure 1 exponential growth is compared with the dynamics characterizing the weak Allee effect, for different value of $c$.

In Figure 1 exponential dynamics growth is compared to the case of weak Allee effect, for increasing values of $c$ to which correspond a decreasing rate of growth.

Consider now the application of a therapy: its effect can be modelled by adding a destruction term $-g x$, to Eqs (2.2) and (2.3), see e.g. ^[28], where $g$ may be either a constant or a function of time. It is straightforward to show that if $g > r$, $x = 0$ is a globally stable point of Eqs (2.3) and (2.2) (indeed $x = 0$ is the only stationary point in $R^+$); in addition in case of Eq (2.2), $x = 0$ is locally stable for any $g > 0$ even though, obviously, its basin of attraction depends on the actual value of $g$. Note that, in both models, cancer can decrease only during medical treatment, to grow again when the therapy is stopped, albeit with different rates.

Obviously the situation is very different if a strong Allee effect is relevant: in this case if $x < A$ cancer population will tend to zero, in other word $x = 0$ is a stable point and its basin of attraction is the segment $[0, A]$. During the therapy, the basin of attraction increases to $[0, A']$, with $A' > A$. Thus, if at the beginning of the therapy $x < A'$, the medical treatment is able to decrease $x$ below the threshold $A$, then the cancer population will continue to decrease even after the therapy has been discontinued.

3.   Two cancer clones

This section will focus on the effects of interaction between populations both subjected to Allee effect and compared with the case of traditional exponential growth, that will be dealt first. In the following $x_1$ and $x_2$ will denote ambiguously the two population and their size.

3.1.   Exponential growth

Basic elements determining the evolution of the cell populations are proliferation, predation and competition for resources: the growth of the cancer species is limited by the amount of available resources as in the single species case and further constrained by inter-specific competition.

The model can be formalized by a system of ordinary differential equations, see for instance ^[25]:

Consider tumor clones $x_1$: the first two terms in the RHS of the equation represent the growth of $x_1$ in isolation, i.e. in absence of other cancer species and medical treatment, compare Eq (2.3). Development of $x_1$ is constrained by the competition with clone $x_2$ ({measured by the competition rate $b_{12}$}); furthermore the medical treatment (for instance chemotherapy) can act on $x_1$ and its effects are represented in the model by the term $g_1(t)x_1$ where $g_1$ takes into account the drugs kinetics in the organism, see ^[18]. This is equivalent to rewrite the growth term as $f_1(t)x_1 = (r_1-g(t))x_1$: when $g(t) > r_1$ then $dx_1/dt < 0$, meaning that the cancer can be eradicated by a given treatment. The same considerations apply mutatis mutandis to tumor species $x_2$. In the following, for simplicity's sake, $f_i(t)$, will be denoted by $f_i$. If $f_i > 0$ the asymptotic states of the system are well known from population theory ^[25]: there exist four possible such states depending on the values of parameters

as follows.

1. $a_{12} < 1$, $a_{21} > 1$: species $x_1$ reaches its carrying capacity $K_1$, whereas $x_2$ goes to $0$ and this trend does not depend on initial conditions. Stationary points are $(K_1, 0)$ (stable), $(0, 0)$ (unstable), and $(0, K_2)$ (unstable).

2. $a_{12} > 1$, $a_{21} < 1$: $x_2$ dominates, irrespective of initial condition. The only stable state is obviously $(0, K_2)$.

3. $a_{12} > 1$, $a_{21} > 1$: again just one species survives depending now on the initial conditions. There are $4$ stationary states : $(K_1, 0)$, $(0, K_2)$ (locally stable) and $(0, 0)$, $(x_1^*, x_2^*)$ unstable with $x_i^* \in (0, K_i)$.

4. $a_{12} < 1$, $a_{21} < 1$: the two species now coexist and the only stable state is of the form $(x_1^*, x_2^*)$, with $x_1^*, x_2^* \neq 0$.

Cases $1-3$ embody the principle of competitive exclusion by which if the competitive interaction are strong enough just one species survives.

Turning now to the effect of therapies, it is clear the only way both $x_1$, $x_2$ can decreases is that $f_i < 0$, $i = 1, 2$ that is $r_i < g_i$; obviously when the medical treatment ends at least one species starts to grow again.

3.2.   Weak Allee effect

Equations for competing species affected by Allee effect can be derived in a way analogous to (3.1), with the obvious modifications; in particular for the weak effect:

The corresponding nullclines are $x_1 = 0, x_2 = 0$ and

In case of weak effect, the main difference with the exponential case of Eq (3.1), is that, in absence of therapy, there exist just two types of asymptotic states, namely dominance depending on the initial conditions or coexistence, as can be shown with standard methods: in other words no species completely dominate the other. Values of $a_{12}$, $a_{21}$ determine the occurrence of the two cases and, if the principle of competitive exclusion applies, the areas of the respective basins of attraction.

A coexistence of two species is depicted in Figure 2, where the two species have the same $c$, but different reproduction rates $r$, with $r_1 > r_2$. The figure also shows the vector field of system (3.2).

When the therapy is applied, the state $(0, 0)$ becomes locally stable for any $g_i > 0$, so, at least for certain initial conditions, both cancer types decrease; discontinuation of therapy results in the growth of at least one species, albeit with a sensibly slower rate than in the exponential case.

3.3.   Strong Allee effect

Equations for populations subjected to a strong Allee effect are similar to the system (3.2), but for the presence of a threshold parameter $A > 0$. A complete analysis of competing populations subjected to the strong Allee effect can be found in ^[31]:

Nullclines of system (3.4) are $x_1 = 0, x_2 = 0$ and

From the analysis of the nullclines, it is straightforward to show that, for any possible set of parameter values, the state $(0, 0)$ is stable; obviously, its basin of attraction, denoted by $B_0$, depends on the parameter values and increases with $A_1$, $A_2$, $b_{12}$, $b_{12}$, whereas decreases for increasing $r_1$, $K_1$, $r_2$, $K_2$.

If the interactions between the two species are strong, there are two possible asymptotic states $(K_1, 0)$, $(0, K_2)$, with basins of attractions $B_1$, $B_2$ respectively, and which species survives depends on the initial conditions, whereas weak competition results in a stable coexistence state with basin $B_c$, compare also ^[31].

The results of a medical treatment, $g_i \neq 0$, is to increase the threshold values $A_1$, $A_2$ and, correspondingly, to widen the basin of attraction of $(0, 0)$, to a new set $B_0'$: clearly $B_0' \supset B_0$. If the therapy can move a point $(x_1(0), x_2(0)) \in B_0'$ to $(x_1(t), x_2(t)) \in B_0$, for some time $t$, then, after discontinuation of the treatment, the cancer population will keep decreasing.

Some speculative examples of the effects of therapies in case of cancer population subjected to a strong Allee effect will be shown the next section in Figure 5, where the same parameters of Figure 3 have been used, with the sole exception of thresholds $A_1, A_2$.

4.   Applications

The experimental evidence for the Allee effect concerns so far the weak case ^[26,20] and in this section we will look for other possible examples of weak effect occurrence. To this aim, we have applied the system (3.2) to data found in ^[24]. As in ^[24], we assume cancer population to be composed of two types, $x_1$ susceptible to a therapy and $x_2$ resistant, while $x = x_1+x_2$ is the total cancer load; note that in ^[24] the data have been fitted with a pair of exponential derived from a probabilistic model: our aim here is not to compare the two fits, but, rather, to show that a model of weak Allee effect can also fit the data, and to obtain some information on the parameter values, and in particular on $c$.

In the following $x_1$ denotes the clone susceptible to medical treatment and $x_2$ is the resistant one; we assume $r_1 > r_2$, representing an evolutionary trade off between rate of growth and resistance to therapy. Initial conditions for all simulations are set as following: $x_1 = 1.2 \cdot 10^8$ and $x_2 = 6 \cdot 10^7$.

The first simulation presents results of the fit of system (3.2) to data of colon-rectal cancer progression reported in ^[24] to which the reader is referred for further details: here, as in ^[24] cancer population is supposed to be composed of two clonal types with different sensitivities to therapies. Results are presented in Figure 3. In the left hand panel, the fit concerns the control case (no therapy). Cancer evolution in which just a type of cancer clones is affected by the therapy is represented in the center panel: here after an initial decline of the total load $x$, cancer grows again, due to the unrestricted increase of the resistant species even after the second application of the same treatment. Finally, the right hand panel depicts a situation in which a combination therapy is adopted and it is able to completely eradicate the cancer. In general, there is a good agreement between data and the numerical prediction of the model. Note that in the left hand panel, in particular, is suggestive of a weak Allee effect in the interval between therapies as the cancer load takes some time before starting to regrow.

A further example consistent with a weak Allee effect is represented by data again reported in ^[24] concerning cancer cells treated with the drug centuximab. As shown in Figure 4, at the beginning of the therapy the total cancer load decreases due to the declining of the sensitive population, but that allows the resistant cells to growth unchecked by the competition and this growth is very slow at low concentration. This fact and the excellent agreement of the results of the model with the data, is, in our opinion, a clear indication of the presence of a weak Allee effect.

Finally we have used the system (3.4) to simulate speculative results of a therapy acting on populations in which a strong Allee effect is present. As it can seen in Figure 5, after the initial increase, if the therapy is able to bring the population under their thresholds, after the end of medical treatment, the total cancer load will keep decreasing to cancer eradication.

5.   Conclusion

Results of the present work clearly support the idea that the occurrence, or otherwise, of the Allee effect at the onset of cancer is relevant not just for a better understanding of cancer dynamics but also for the choice of the most appropriate therapy, ^[22]. Indeed, in presence of a strong effect the administration of a maximum tolerated dose can be preferable to a therapy aiming just to control cancer in that the former can push the number of tumor cells below the threshold thus ensuring eradication. In case of weak effect tumor will necessarily grow again after therapy discontinuation, however MTD should ensure a long period of relatively low count of cancer cells.

The problem is, therefore, to find suitable means to determine the use whether a given type of cancer is affected by an Allee effect. Obviously solving this problem is simpler if it is possible to perform "in vitro" experiments; in particular, it is easier to demonstrate that a threshold exists (strong effect) rather than there is no threshold, and not that the study fails to reveal one ^[7]. In principle, changing the initial density of the cells plated in vitro one can try to pinpoint the size of the population for which the per capita population growth rate becomes negative. However, estimates of the Allee threshold are usually affected by errors due to fluctuations which are particularly relevant at low cells abundance (compare ^[7] and citations therein). Discovering a weak effect is more complicated as it depends on the analysis of time series data, as shown in ^[20,26], and also by the simulations presented here, but of course results are less clear cut than in case of a strong effect, as they depend crucially on the choice of models, and alternative interpretations are possible. Obviously, the presence of an Allee effect is more difficult to find in clinical situations: to find the occurrence of the effect requires the observation of cancer at low population density and this is not always possible during clinical trials and furthermore the effect can be masked by other factors, such as the action of the immune system ^[28].

Acknowledgments

This work is supported by the "Departments of Excellence 2018 - 2022" Grant awarded by the Italian Ministry of Education, University and Research (MIUR) (L. 232/2016).

Conflict of interest

All authors declare no conflicts of interest in this paper.