Averaging principle for space-fractional stochastic partial differential equations driven by Lévy white noise and fractional Brownian motion

1.   Introduction

Cancer is a major public health problem and a leading cause of death in every country of the world ^[1]. Promising clinical approaches that could counter cancer are urgently needed, such as targeted therapy. After decades of development, it has become an effective and standard way for the anticancer treatment ^[2]. However, drug resistance to targeted therapy often strongly limits the long-term effectiveness of treatment and reduces the survival rates for cancer patients ^[3]. Therefore, drug resistance is one of the major obstacles to improving the chances of clinical anticancer treatment success.

Diverse forms of approaches have been proposed that aim to better understand and eventually overcome drug resistance ^[4], one of which, of course, is the experimentally reveal mechanisms of drug resistance. Although the cancer resistance mechanisms are very complex, several of them have been revealed ^[5,6,7,8,9]. One mechanism of resistance generation is the result of genetic events such as mutations. Studies have shown that two main types of mutations, gene amplification and point mutation, lead to cells drug resistance ^[10,11,12]. Gene amplification is the result of an overproduction of one or more specific genes, which means that a parts of the genome is copied to a much greater extent than the replication of DNA composing the rest of the genome. This defect expands the phenotype that the gene gives to the cells, which in turn, induces resistance by giving the cells with more copies of a specific gene than the drug can handle. However, the frequency of gene amplification occurs is much lower than point mutation^[10]. Point mutation prevents the drug from successfully binding to the target protein, which is a permanent phenomenon and makes its progeny cells resistant to the drug. It is this type of drug resistance that we consider in this paper.

Experimental studies in vitro and in vivo have made rapid progress in revealing the mechanisms of drug resistance. One may pay attention to such problems as: how does the generation of drug resistance depend on the tumor size, mutation rate and the number of drugs? Does resistance arise before the beginning of the treatment or mainly during treatment? How does the possibility of resistance generation change depend on the dosage of the drug? Whether combination therapy have an advantage over monotherapy? Independent experimental studies are expensive and time-consuming because of multiple cell types and drug dosages need be considered except for other experimental conditions and challenges. As an alternative, mathematical modeling is a powerful tool that may helps us address these problems.

Many mathematical models have been proposed to study drug resistance based on biological mechanisms, such as ordinary differential equation (ODE) model ^[13], stochastic model ^[14,15], stochastic differential equation (SDE) model ^[8], partial differential equation (PDE) model ^[16] or other models ^[17,18,19]. The first model of drug resistance due to point mutation can be traced back to the classical study by Goldie and Coldman ^[20] in 1979, who analyzed the probability of the appearance and evolution of resistant by using stochastic processes (see ^[21,22]). After that, more results on point mutation are obtained by the authors for the case of one or multi-drug resistance ^{[23,24,25,26]}. For example, in ^[24], a stochastic model is developed by Komarova to analyze the dependence of the probability of resistance generation on the initial tumor load, mutation rate and turnover rate. This model showed how the pre-existence of resistance is more important than the generation of during-treatment. Using birth and death processes, Foo and Michor ^[26] described the evolution of resistance during treatment due to a single (epi)genetic alteration. They further provided a methodology for predicting the risk of resistance under various dosing strategies. Another recent work is by Bozic et al.^[27], in which branching processes are used to describe the evolutionary dynamics of lesions in response to treatment. They improved on the results previously obtained in these works and provided a closer match to simulations. Although these models have greatly improved our understanding of cancer progression, the stochastic methods used are often have their complexity. Tomasetti and Levy ^[28] built a simple ODE-based model to obtain comparable results to those using complex techniques of stochastic model. The advantage of this model is that it is possible to analyse the results of resistance to any number of drugs. However, this framework did not address the relationship among drug resistance and dosage of the drug and evaluate the efficiency of combination drug therapy. It is desirable to develop mathematical models that can fully evaluate the evolution of drug resistance and the efficiency of combination therapy.

In this paper, we improve the mathematical model proposed by Tomasetti and Levy ^[28] to connect drug resistance to dosage of the drug. The basic modeling assumptions and ODE-based model are introduced in section 2. In section 3, we show the main results obtained from the model analysis: the dependence of resistance present at time of the beginning of the treatment on turnover rate, the relative roles between pre-treatment phase and treatment phase, and the efficiency of combination therapy are discussed. Finally, we finish the paper with conclusions and future directions in section 4.

2.   Mathematical model

In this section, we propose a linear system which consists of a set of ordinary differential equations for cancer drug resistance in the context of multi-drug therapy. This model is an extension of previous work by Tomasetti and Levy ^[28], the main difference is that we improve the model by a more accurate Hill function dosing model, so as to analyze the influence of drug concentration changes on the amount of resistance generated in pre-treatment phase and during treatment.

2.1.   The conceptual framework

In the course of cancer progression, each cell has a certain probability to undergo death, division, mutation, or transformation. These different processes can lead to an increase or decrease in the number of cells of the given type, see Figure 1. However, mutations may give rise to different types of cells that are resistant to drugs. To better understand how the process of generation of resistance to the drug, we present the following conceptual framework. If we only consider treatment with one drug, then there are two types of cells in the system that are either susceptible or resistant to this drug. We can describe this process by the network $x_1\xrightarrow{\mu} x_0$, where $\mu$ is the mutation rate, the subscript $j$ in notation $x_j$ represents susceptibility ($j = 1$ means susceptible, and $j = 0$ not susceptible to this drug).

For the treatment with $n$ different drugs, a cell has to accumulate $n$ mutations to become resistant to all drugs. We assume that there are $n$ types of mutations, and the mutation rates are denoted by $\mu_1$, $\mu_2$, ..., $\mu_n$. Each mutation rate $\mu_i$ corresponding to a phenotype resistant to drug $i$. Thus, it is easy to find that different resistant phenotypes can be up to $m = 2^n-1$. As an example, the mutational processes for $n = 3$ can be presented by a combinatorial mutation network, which is shown in Figure 2. Each node corresponds to a phenotype, which is described by the binary number of length $n$: "0" represents susceptible to the drug corresponding to its position, "1" represents resistant. The symbol below the nodes represents the corresponding variable. For example, 010 denotes this phenotype is resistant to drug 2 but not to drug 1 and 3. 000, 111 denote fully susceptible or resistant to all drugs, respectively.

2.2.   The case of one drug

Let us suppose that we have one drug. Denote by $N(t)$ and $R(t)$ the number of wild-type cancer cells (drug-sensitive cells) and resistant cells (resistance is acquired by means of a mutation) at time $t$, respectively. Our model is based on the following assumptions:

$(1)$ According to the Skipper-Schabel-Wilcox model, we assume that cancer cells follow exponential growth as in ref. ^[28].

$(2)$ Denote the birth, death, and mutation rate as $l$, $d$ and $\mu$, respectively. Assume that $0 < \mu\ll1$, the condition $0\leq d < l$ corresponds to a clonal expansion. The ratio $0\leq d/l < 1$ defines the turnover rate of cancer cells. We call the scenario where $d/l\ll1$ a low-turnover, low death cancer. In contrast, the scenario where $d/l\approx1$ describes extremely high-turnover, high death cancer. To simplify the initial presentation, we consider a symmetrical case, assume that the birth rate, death rate and mutation rate are the same for all types cells.

$(3)$ In this paper, we only consider the scenario of non-mutagenic drugs. The drug-induced death rate $\tilde{h}$ depends on drug concentration (D) using Hill function $\tilde{h} = h\cdot D/(K+D)$ as in ref. ^[8], where $h$ represents maximal death rate, and $K$ is a Michaelis constant representing the drug concentration associated with reaching the half-maximal inhibition effect.

$(4)$ Assume that resistance is acquired by means of a point mutation, which happens only in one direction. As a result, a drug-sensitive cell differentiating into one drug-sensitive and one mutant cell ^[24].

With these assumptions, the model given by Tomasetti and Levy ^[28] can be modified as follows:

where the drug therapy starts at time $t^*$. $D_1$, $K_1$ represent concentration and Michaelis constant of the first drug, respectively. The system (2.1) describes the pre-treatment dynamics with initial condition $N(0) = N_0 > 0$ and $R(0) = 0$, and system (2.2) describes the dynamics after the treatment starts with initial conditions $N(t^*)$ and $R(t^*)$, which are the solutions of system (2.1) at $t = t^*$.

Remark 1. In the present model, we assume that resistant cells behave in the same way as the sensitive cells. This may not be the case, it is possible that resistant cells have a relative fitness advantage or disadvantage compared to the sensitive cells. See for example ^[25], if the birth and death rate are $l_1$, $d_1$ for sensitive cells, $l_2$, $d_2$ for resistant cells, respectively. Then the relative fitness is given by $\alpha = (l_2-d_2)/(l_1-d_1)$, resistant cells can be advantageous ($\alpha > 1$), neutral ($\alpha = 1$), or deleterious ($\alpha < 1$). Therefore, the relative fitness considered in this paper is only the neutral case. It can be modified according to different situations.

2.3.   The case of two drugs

Let us consider the treatment with two drugs. Denote by $R_1(t)$ and $R_2(t)$ the number of mutant cancer cells that they are resistant only to the first or to the second drug at time $t$, respectively. The number of mutated cells that are resistant to both drugs at time $t$ is denoted by $R(t)$. Before presenting our model, we have the following additional assumptions in addition to the case of one drug:

$(5)$ Assume that two drugs have similar effects when inducing the death of drug-sensitive cancer cells, namely, both two drugs can increase the death rate of drug-sensitive cells. Using Hill function, the death rate of drug-sensitive cancer cells following treatment can be given by $\tilde{h} = h\cdot\big(\frac{D_1}{K_1+D_1}+\frac{D_2}{K_2+D_2}\big)$, where $D_2$, $K_2$ represent concentration and Michaelis constant of the second drug, respectively.

$(6)$ Assume that for any cells that are non fully resistant state, if they are resistant to only one drug, the second drug is still effective and can kill the cells with the maximum rate. That is, the mutant cells that are resistant to only the first drug, they can still killed by the second drug with the rate $\frac{h\cdot D_2}{K_2+D_2}$; similarly, mutant cells that are resistant to only the second drug, they can still killed by the first drug with the rate $\frac{h\cdot D_1}{K_1+D_1}$.

Under these assumptions, the model for drug resistance with two drugs is given by

Similarly, in the case of treatment with two drugs, the system (2.3) and (2.4) describe the dynamics of pre-treatment phase and after the treatment starts, respectively.

Remark 2. To model the clinically relevant situation where we start treatment when the tumor reaches a detectable size, we use the following trick to estimate the time of the beginning of the treatment $t^*$. Assume that the initial number of cancer cells is $N_0$ ($N_0\gg1$) instead of one cell and the total number of cancer cells at time $t^*$ is $M$. Using eigenvalue method, we can easily obtain the solutions of system (2.1) are

By $N(t^*)+R(t^*) = M$ and the mutation rate $\mu\ll1$, the time $t^*$ can be estimated as

In this way, we can extend our model to the case of treatment with any number of drugs are being simultaneously used. As an example, we show the model of treatment with three drugs in detail in S.1 of Supplementary information. While the basic mathematical tools that we used allow us to study the case of $n$-drug, it is undeniable that any deterministic model has its complexity and limited validity in mathematical processing as the number of drugs increases. Moreover, the combination of more than three drugs is less practical in clinical application and may cause unknown side effects. Therefore, we mainly focus on the combination of two drugs in this paper.

3.   Results

3.1.   Dependence on the turnover rate $d/l$

In this section, we investigate how does the generation of resistance at time of the beginning of the treatment $t^*$ depends on the turnover rate $d/l$. Furthermore, we analyse the effect of the mutation rate and detection size on resistance.

First, consider the case of one drug only. For the system (2.1), the solution of $R(t^*)$ with initial condition $R(0) = 0$ is given by

In view of (2.5) to evaluate $t^*$ yields

Therefore, we can see for one-drug treatment, the resistance present at time of the beginning of the treatment is dependent of the turnover rate $d/l$.

For the case of two drugs, the solutions of $R_1(t^*)$ and $R_2(t^*)$ are given by

By using eigenvalue method, we obtain

Similarly, if we extend to the case of treatment with $n$ drugs, then

Another expression for the number of cells resistant to $n$ drugs with no cross resistance can be found in ^[27] (Eq (11) in the Supplement). Now we can obtain the same result, that is, the amount of resistance present at time $t^*$ always depends on the turnover rate $d/l$. This dependence becomes increasingly stronger with the increase of the number of drugs, see Figure 3. For low-turnover rate, lower-death cancer, the resistance in the pre-treatment is smaller. In contrast, for high-turnover rate, higher-death cancer, the larger is the resistance present before treatment. It should be noted that for extremely high turnover rate ($d/l\approx1$) such that $\frac{\mu\ln M/N_0}{l(1-d/l)} > 1$, combination therapy is unlikely to yield an advantage than monotherapy.

In addition, with the increase of the number of drugs, we find two important differences:

● The amount of resistance now depends on the mutation rate, $\mu$. If the mutation rate is higher, the larger is the resistance present when the tumor reaches a detectable size. The larger the number of drugs, the stronger this dependency (Figure 4(a)).

● An increase in the detection size, $M$, results in a larger resistance. Increasing the number of drugs, the smaller the resistance is generated (Figure 4(b)).

In the simulations above, the parameters we used within biologically reasonable ranges, which are collected or estimated from previous studies. The details are described in S.2 of Supplementary information. The biological meaning underlying the parameters along with their units, estimated values and reference sources are listed in Supplementary Table S.1.

3.2.   Generation of resistance before and after treatment

In this section, we compare the generation of resistance before treatment and during treatment. In other words, at what stage is resistance mainly generated: before treatment starts, or during treatment? How does the possibility of resistance generation change depend on the dosage of the drug? In order to do this, we introduce two quantities, $R_i^\uparrow(t)$ and $R_i^\downarrow(t)$. If we artificially set the mutation rate to zero after time of the beginning of the treatment, then the only drug resistance that is present at time $t$ $(t > t^*)$ comes from before treatment. We call such resistance as the "pre-treatment resistance at time $t$", this is denoted by $R_i^\uparrow(t)$ (the symbol $\uparrow$ represents the contribution of the pre-treatment phase to resistance generation, and the subscript $i$ indicates the number of drugs we use). Next, we set the mutation rate to zero in the pre-treatment phase, and turn it back after the treatment starts. Thus, resistance develops only during treatment. We define such resistance as the "during-treatment resistance at time $t$", $R_i^\downarrow(t)$ (the symbol $\downarrow$ represents the contribution of the treatment phase to resistance generation). Now, we calculate and compare the quantities $R_i^\uparrow(t)$ and $R_i^\downarrow(t)$.

For one drug, it is clear that the pre-treatment resistance $R_1^\uparrow(t)$ is the solution of $R(t)$ in system (2.2) with the initial condition

thus the expression of $R_1^\uparrow(t)$ is given by

The resistance that is generated only during treatment $R_1^\downarrow(t)$, is the solution of $R(t)$ in system (2.2), subject to the initial conditions $N(t^*) = M$ and $R(t^*) = 0$, we obtain

As we can see if $R_1^\uparrow(t) > R_1^\downarrow(t)$ is satisfied for any $t > t^*$, there must be

Therefore, we have

This result holds under the assumption that $\frac{M}{N_0}\geq e$. On the other hand, in order to ensure treatment intensity is in the biologically relevant range, we have $h > (l-d)(K_1+1)$. For convenience, in the analysis below, we choose $h\geq2(l-d)$ is realistic. With this, our result agree with ^[24,23] also on the selection of threshold value of treatment intensity (see page 361 of ^[24] and page 9715 of ^[23]). There are two conclusions from these calculations:

● In the case of one drug, the generation of resistance before and at a specific time after the beginning of treatment depends on the turnover rate $d/l$ (see Eq (3.4)).

● Under the realistic treatment intensity $D_1 \geq \frac{K_1(l-d)}{h+d-l}$ (assuming that $M/N_0\geq e$), we have $R_1^\uparrow(t) > R_1^\downarrow(t)$, that is, resistance is mainly generated before the beginning of the treatment.

In the case of two drugs, $R_2^\uparrow(t)$ is the solution of $R(t)$ in system (2.4) with the initial condition

and therefore

On the other hand, the solution of the system (2.4) with the initial conditions

is given by

Thus for any $t > t^*$, if $R_2^\uparrow(t) > R_2^\downarrow(t)$ is true, then

Now, it is easy to verify that under the treatment intensity

we have $R_2^\uparrow(t) > R_2^\downarrow(t)$. This result holds under the assumption that $\frac{M}{N_0}\geq e^{\sqrt{2}}$. We conclude the following:

● Eq (3.7) shows that the dependence of resistance generated before and at a specific time after the treatment starts on the turnover rate $d/l$ becomes increasingly stronger with the increase of the number of drugs.

● For two drugs, the $R_2^\uparrow(t) > R_2^\downarrow(t)$ holds under the realistic treatment intensity $D_1 \geq \frac{K_1(l-d)}{h+d-l}$ and $D_2 \geq \frac{K_2(l-d)}{h+d-l}$ (assuming that $M/N_0\geq e^{\sqrt{2}}$). Hence, the generation of resistance in pre-treatment phase always plays the dominant role than during-treatment.

Indeed, the result of $R_i^\uparrow(t) > R_i^\downarrow(t)$ can be given a intuitive explanation. Before the beginning of the treatment, cells have a clonal expansion because of $l > d$, which certainly facilitates the generation of resistance. However, during treatment, cells have a negative growth rate since the treatment intensity $D_i \geq \frac{K_i(l-d)}{h+d-l}$ such that $\frac{h\cdot D_i}{K_i+D_i}+d > l$. This makes it less likely that more mutations can be acquired. Therefore, the generation of resistance before the beginning of the treatment is greater compared to the resistance created during treatment.

It is worth noting that in the discussion above, we considered that all mutations occurred before or after treatment. There is another situation that cannot be ignored, that is, some mutations happen before treatment and others happen after. The question might becomes: what happens when the secondary mutation occurs before or after treatment? In fact, we still have the same conclusion, that is, resistance is primarily created before treatment. The details on the derivation of this question are given in S.3 of Supplementary information.

3.3.   Evaluation of drug combination effect

The focus of this section is to evaluate the efficiency of drug combination therapy. In the following study, we use numerical simulations to examine three treatment strategies: drug-1 alone (Figure 5(a)) or with a combination of drug-1 and drug-2 at doses of 150 mg and 1 mg (1x; Figure 5(b)) or 150 mg and 2 mg (2x; Figure 5(b)). Note that the drug doses have been normalized and we use dimensionless concentrations of drugs in the simulations. On the other hand, we would like to stress that the parameter values used in this section are purely for numerical illustration, and do not represent specific model fits or therapies.

As shown in Figure 5, either treated with single drug or a combination of two drugs, an appropriately elevated dosage of the drug can reduce the number of resistant cells. Furthermore, by comparing Figure 5(a) with 5(b), we find that the number of resistant cells $R(t)$ in Figure 5(b) is several orders of magnitude lower than in Figure 5(a). Therefore, our model predict that combination therapy has a significant advantage over single drug therapy.

However, this result may not hold for very high-turnover cancers. For example, if we assume the death rate $d = 0.199$, then the turnover rate $d/l = 0.995\approx1$. The other parameter values are consistent with above. Next, we use the same way as in Figure 5 to compare different treatment strategies. The simulation results are shown in Figure 6.

From Figure 6, it can be seen that the number of resistant cells $R(t)$ in Figure 6(b) is more than in Figure 6(a), although they are of the same order of magnitude. Therefore, in this case, combination therapy is unlikely to yield an advantage than monotherapy, which also illustrates the result in section 3.1.

4.   Discussions

This section is designed to comment on the results obtained with our model and those obtained in ^{[23,24,27,28]}. In order to model the evolution of drug resistance caused by stochastic genetic point mutations, Tomasetti and Levy ^[28] proposed a basic ODE model. The advantage of their method lay in the simplicity of the model, so that it could obtain drug resistance results comparable to those that were obtained by using complicated stochastic methods before (e.g, in ^[23,24,27]). However, the drug-induced death rate considered in ^[28] is a constant, also in ^[23,24,27]. Drug-induced death rate depends on the concentration of the drug actually, different death rate is caused by different drug concentration. In this paper, we improve the model proposed by Tomasetti and Levy ^[28], establish a more accurate Hill function dosing model, and analyse the effect of drug concentration changes on the level of drug resistance before and after the treatment. Compared with the previous work, the main differences of this paper are as follows:

(i) Our results agree with ^[28], which clearly show that no matter how many drugs are used in the treatment, the amount of resistance (generated at the beginning of the treatment and within a certain period of time after the treatment) always depends on the turnover rate, this result is independent of the drug concentration (see Eq (3.1)–(3.3)). In contrast, one of the main results of ^[23,24] was that this dependence holds only in the case of multi-drugs, but not in the case of single drug (see page 9715 of ^[23] and page 360 of ^[24]). Furthermore, we also analyse the effect of the mutation rate and detection size on the amount of resistance. Both of these effects become stronger with the increase of the number of drugs (see Figure 4).

(ii) In this work, we investigate the relative role of the pre-treatment phase and during-treatment in the development of drug resistance. Our results indicate that the resistance generated before the beginning of the treatment is greater than that developed during-treatment. This result depends on the concentration of the drug, which holds only when the concentration of the drug reaches a lower limit (see Eq (3.6)). From this perspective, the drug-induced death rate $h$ must satisfy $h > 2(l-d)$ (see section 3.2). However, the condition for this result in ^[28] to hold was $h > l-d$ (see page 6). The reason for this difference is that ^[28] did not consider the effect of drug concentration on drug resistance, neither in ^[24,27].

(iii) We improve the work of ^[28] and evaluate the effect of drug combination. The results show that regardless of single-drug therapy or two-drug combination therapy, the amount of resistant mutants can be reduced by appropriately increasing the drug concentration (see Figure 5). More importantly, our model predict that for cancer with low-turnover rate, combination therapy has significant advantage over monotherapy, which can further reduce drug resistance. However, it should be noted that for cancer with very high-turnover rate (close to 1), combination therapy can not significantly reduce the amount of resistant mutants compared to monotherapy, so in this case, combination therapy would not have advantage over monotherapy (see Figure 6). This result is consistent with ^[24,27], although the mathematical methods used are different. We use an elementary, compartmental ODE system, while ^[23,24,27] adopted various stochastic methods.

In the final, we would like to reveal several interesting extensions associated with the model and drug resistance. In our model, the cancer cells are modeled by an exponential growth as in ref. ^[28]. It is a reasonable assumption for small number of cancer cells with early growth. In fact, the growth of cancer cells slows when the population is large, this growth can be modeled by a more realistic model such as the Gompertz growth, as Afenya and Calderón stated that this is best for describing chronic myelogenous leukemia (CML) growth ^[29]. Another growth can be replaced by logistic growth, which has been found to provide a better fit in numerous studies ^[8,30].

As assumed in section 2.2 this work is carried out for non-mutagenic drugs. However, some resistant is actually generated by the drug-induced. For example, Obenauf et al. ^[31] showed that drug-sensitive cancer cells secrete a variety of resistance factors (such as IGF and HGF) under the effect of targeted therapy, which can facilitate the growth, dissemination and metastasis of cancer cells and further increase drug resistance. In this case, the balance of the resistance generated between pre-treatment and during-treatment may be reversed. Therefore, it is worth considering the drug-induced resistance into our framework.

Drug resistance is a common obstacle during targeted therapy, which may be associated with eventual treatment failure. Although one of the most effective treatment methods currently is the use of a variety of targeted therapy strategies, including monotherapy or multi-drug combination therapy, intermittently or continuously therapy, the design of the optimal treatment therapies (e.g., timing, sequence) and the selection of the optimal dosages remains a challenge. It is our hope that a more reasonable mathematical prediction model will be proposed to facilitate the design of clinical therapies and we left it for future research.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Nos. 11871238).

Conflict of interest

The author declares there is no conflict of interest.