Research article

An elementary mathematical modeling of drug resistance in cancer

  • Received: 07 September 2020 Accepted: 10 November 2020 Published: 02 December 2020
  • Targeted therapy is one of the promising strategies for the treatment of cancer. However, resistance to anticancer drug strongly limits the long-term effectiveness of treatment, which is a major obstacle for successfully treating cancer. In this paper, we analyze a linear system of ordinary differential equations for cancer multi-drug resistance induced mainly by random genetic point mutation. We investigate 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. Moreover, 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. Using numerical simulations, we also evaluate the response of the mutant cancer cell population as a function of time under different treatment strategies. At appropriate dosages, combination therapy produces significant effects for the treatment with low-turnover rate cancer. 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.

    Citation: Kangbo Bao. An elementary mathematical modeling of drug resistance in cancer[J]. Mathematical Biosciences and Engineering, 2021, 18(1): 339-353. doi: 10.3934/mbe.2021018

    Related Papers:

    [1] Cristian Tomasetti, Doron Levy . An elementary approach to modeling drug resistance in cancer. Mathematical Biosciences and Engineering, 2010, 7(4): 905-918. doi: 10.3934/mbe.2010.7.905
    [2] Natalia L. Komarova . Mathematical modeling of cyclic treatments of chronic myeloid leukemia. Mathematical Biosciences and Engineering, 2011, 8(2): 289-306. doi: 10.3934/mbe.2011.8.289
    [3] Rujing Zhao, Xiulan Lai . Evolutionary analysis of replicator dynamics about anti-cancer combination therapy. Mathematical Biosciences and Engineering, 2023, 20(1): 656-682. doi: 10.3934/mbe.2023030
    [4] Alexis B. Cook, Daniel R. Ziazadeh, Jianfeng Lu, Trachette L. Jackson . An integrated cellular and sub-cellular model of cancer chemotherapy and therapies that target cell survival. Mathematical Biosciences and Engineering, 2015, 12(6): 1219-1235. doi: 10.3934/mbe.2015.12.1219
    [5] Haifeng Zhang, Jinzhi Lei . Optimal treatment strategy of cancers with intratumor heterogeneity. Mathematical Biosciences and Engineering, 2022, 19(12): 13337-13373. doi: 10.3934/mbe.2022625
    [6] Hongli Yang, Jinzhi Lei . A mathematical model of chromosome recombination-induced drug resistance in cancer therapy. Mathematical Biosciences and Engineering, 2019, 16(6): 7098-7111. doi: 10.3934/mbe.2019356
    [7] Cameron Meaney, Gibin G Powathil, Ala Yaromina, Ludwig J Dubois, Philippe Lambin, Mohammad Kohandel . Role of hypoxia-activated prodrugs in combination with radiation therapy: An in silico approach. Mathematical Biosciences and Engineering, 2019, 16(6): 6257-6273. doi: 10.3934/mbe.2019312
    [8] Damilola Olabode, Libin Rong, Xueying Wang . Stochastic investigation of HIV infection and the emergence of drug resistance. Mathematical Biosciences and Engineering, 2022, 19(2): 1174-1194. doi: 10.3934/mbe.2022054
    [9] Pedro José Gutiérrez-Diez, Jose Russo . Design of personalized cancer treatments by use of optimal control problems: The case of chronic myeloid leukemia. Mathematical Biosciences and Engineering, 2020, 17(5): 4773-4800. doi: 10.3934/mbe.2020261
    [10] Ami B. Shah, Katarzyna A. Rejniak, Jana L. Gevertz . Limiting the development of anti-cancer drug resistance in a spatial model of micrometastases. Mathematical Biosciences and Engineering, 2016, 13(6): 1185-1206. doi: 10.3934/mbe.2016038
  • Targeted therapy is one of the promising strategies for the treatment of cancer. However, resistance to anticancer drug strongly limits the long-term effectiveness of treatment, which is a major obstacle for successfully treating cancer. In this paper, we analyze a linear system of ordinary differential equations for cancer multi-drug resistance induced mainly by random genetic point mutation. We investigate 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. Moreover, 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. Using numerical simulations, we also evaluate the response of the mutant cancer cell population as a function of time under different treatment strategies. At appropriate dosages, combination therapy produces significant effects for the treatment with low-turnover rate cancer. 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.


    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.

    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.

    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 x1μx0, where μ is the mutation rate, the subscript j in notation xj represents susceptibility (j=1 means susceptible, and j=0 not susceptible to this drug).

    Figure 1.  The four basic processes: death, division, mutation and transformation (applied to one of the "red" cells).

    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 μ1, μ2, ..., μn. Each mutation rate μi corresponding to a phenotype resistant to drug i. Thus, it is easy to find that different resistant phenotypes can be up to m=2n1. 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.

    Figure 2.  Mutation network for three drugs (see text for details, also see refs. [23,24]).

    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 μ, respectively. Assume that 0<μ1, the condition 0d<l corresponds to a clonal expansion. The ratio 0d/l<1 defines the turnover rate of cancer cells. We call the scenario where d/l1 a low-turnover, low death cancer. In contrast, the scenario where d/l1 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 ˜h depends on drug concentration (D) using Hill function ˜h=hD/(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:

    {N(t)=(ld)N(t),R(t)=(ld)R(t)+μN(t).tt. (2.1)
    {N(t)=(ldhD1K1+D1)N(t),R(t)=(ld)R(t)+μN(t)t>t. (2.2)

    where the drug therapy starts at time t. D1, K1 represent concentration and Michaelis constant of the first drug, respectively. The system (2.1) describes the pre-treatment dynamics with initial condition N(0)=N0>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 l1, d1 for sensitive cells, l2, d2 for resistant cells, respectively. Then the relative fitness is given by α=(l2d2)/(l1d1), resistant cells can be advantageous (α>1), neutral (α=1), or deleterious (α<1). Therefore, the relative fitness considered in this paper is only the neutral case. It can be modified according to different situations.

    Let us consider the treatment with two drugs. Denote by R1(t) and R2(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 ˜h=h(D1K1+D1+D2K2+D2), where D2, K2 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 hD2K2+D2; similarly, mutant cells that are resistant to only the second drug, they can still killed by the first drug with the rate hD1K1+D1.

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

    {N(t)=(ld)N(t),R1(t)=(ld)R1(t)+μN(t),R2(t)=(ld)R2(t)+μN(t),R(t)=(ld)R(t)+μR1(t)+μR2(t).tt. (2.3)
    {N(t)=(ldh(D1K1+D1+D2K2+D2))N(t),R1(t)=(ldhD2K2+D2)R1(t)+μN(t),R2(t)=(ldhD1K1+D1)R2(t)+μN(t),R(t)=(ld)R(t)+μR1(t)+μR2(t).t>t. (2.4)

    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 N0 (N01) 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

    N(t)=N0e(ld)tandR(t)=N0μte(ld)t

    By N(t)+R(t)=M and the mutation rate μ1, the time t can be estimated as

    t1ldlnMN0. (2.5)

    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.

    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

    R(t)=N0μte(ld)t.

    In view of (2.5) to evaluate t yields

    R(t)=N0μte(ld)tMμlnMN0l(1dl). (3.1)

    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 R1(t) and R2(t) are given by

    R1(t)=N0μte(ld)t,R2(t)=N0μte(ld)t.

    By using eigenvalue method, we obtain

    R(t)=N0(μt)2e(ld)tM[μlnMN0l(1dl)]2. (3.2)

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

    R(t)=N0(μt)ne(ld)tM[μlnMN0l(1dl)]n. (3.3)

    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/l1) such that μlnM/N0l(1d/l)>1, combination therapy is unlikely to yield an advantage than monotherapy.

    Figure 3.  The dependence of the amount of resistance present at time t on turnover rate d/l. The higher the turnover rate, the larger is the resistance. This dependence becomes stronger with the increase of the number of drugs. Parameter values are chosen as follows: l=0.2, M=109, N0=102, μ=104.

    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, μ. 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)).

    Figure 4.  The amount of resistance present at time t depends on the parameters of the model. (a) Dependence on the mutation rate, μ. The higher the value of μ, the larger is the resistance. The larger the number of drugs, the stronger this dependency. Parameter values are: l=0.2, d=0.1 M=109, N0=102. (b) Dependence on the detection size, M. The larger the value of M and smaller the number of drugs, the larger is the resistance. Parameter values are: l=0.2, d=0.1, N0=102, μ=104.

    ● 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.

    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, Ri(t) and Ri(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 Ri(t) (the symbol 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", Ri(t) (the symbol represents the contribution of the treatment phase to resistance generation). Now, we calculate and compare the quantities Ri(t) and Ri(t).

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

    R(t)=MμlnMN0ld,

    thus the expression of R1(t) is given by

    R1(t)=MμlnMN0lde(ld)(tt). (3.4)

    The resistance that is generated only during treatment R1(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

    R1(t)=Mμ(K1+D1)hD1[e(ld)(tt)e(ldhD1K1+D1)(tt)]. (3.5)

    As we can see if R1(t)>R1(t) is satisfied for any t>t, there must be

    MμlnMN0ldMμ(K1+D1)hD1.

    Therefore, we have

    D1K1(ld)h+dl, (3.6)

    This result holds under the assumption that MN0e. On the other hand, in order to ensure treatment intensity is in the biologically relevant range, we have h>(ld)(K1+1). For convenience, in the analysis below, we choose h2(ld) 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 D1K1(ld)h+dl (assuming that M/N0e), we have R1(t)>R1(t), that is, resistance is mainly generated before the beginning of the treatment.

    In the case of two drugs, R2(t) is the solution of R(t) in system (2.4) with the initial condition

    R(t)=M[μlnMN0ld]2,

    and therefore

    R2(t)=M[μlnMN0ld]2e(ld)(tt). (3.7)

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

    {N(t)=M,R1(t)=0,R2(t)=0,R(t)=0,

    is given by

    R2(t)=Mμ2(K1+D1)(K2+D2)h2D1D2[e(ld)(tt)e(ldhD1K1+D1)(tt)e(ldhD2K2+D2)(tt)]+Mμ2(K1+D1)(K2+D2)h2D1D2e(ldhD1K1+D1hD2K2+D2)(tt). (3.8)

    Thus for any t>t, if R2(t)>R2(t) is true, then

    M[μlnMN0ld]22Mμ2(K1+D1)(K2+D2)h2D1D2.

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

    D1K1(ld)h+dlandD2K2(ld)h+dl, (3.9)

    we have R2(t)>R2(t). This result holds under the assumption that MN0e2. 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 R2(t)>R2(t) holds under the realistic treatment intensity D1K1(ld)h+dl and D2K2(ld)h+dl (assuming that M/N0e2). Hence, the generation of resistance in pre-treatment phase always plays the dominant role than during-treatment.

    Indeed, the result of Ri(t)>Ri(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 DiKi(ld)h+dl such that hDiKi+Di+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.

    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.

    Figure 5.  Evaluation of different treatment strategies effects. We compare the amount of resistance present at time t with (a) no treatment and drug-1 monotherapy, (b) a combination of 150 mg drug-1 and 1 mg drug-2 (Drug-1 & Drug-2, 1x), and a combination of 150 mg drug-1 and 2 mg drug-2 (Drug-1 & Drug-2, 2x). The time t=0 refers to t=t. l=0.2, d=0.1, other parameter values are listed in Table S.1.

    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.9951. 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.

    Figure 6.  Evaluation of different treatment strategies effects. We compare the amount of resistance present at time t with (a) no treatment and drug-1 monotherapy, (b) a combination of 150 mg drug-1 and 1 mg drug-2 (Drug-1 & Drug-2, 1x), and a combination of 150 mg drug-1 and 2 mg drug-2 (Drug-1 & Drug-2, 2x). The time t=0 refers to t=t. l=0.2, d=0.199, other parameter values are the same as above.

    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.

    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(ld) (see section 3.2). However, the condition for this result in [28] to hold was h>ld (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.

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

    The author declares there is no conflict of interest.



    [1] F. Bray, J. Ferlay, I. Soerjomataram, R. L. Siegel, L. A. Torre, A. Jemal, Erratum: Global cancer statistics 2018: globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries, CA Cancer J. Clin., 68 (2020), 394-424.
    [2] C. Sawyers, Targeted cancer therapy, Nature, 432 (2004), 294-297.
    [3] C. Brown, Targeted therapy: an elusive cancer target, Nature, 537 (2016), S106-S108.
    [4] X. Sun, B. Hu, Mathematical modeling and computational prediction of cancer drug resistance, Briefings Bioinf., 19 (2018), 1382-1399. doi: 10.1093/bib/bbx065
    [5] X. Hu, Z. Zhang, Understanding the genetic mechanisms of cancer drug resistance using genomic approaches, Trends Genet., 32 (2016), 127-137. doi: 10.1016/j.tig.2015.11.003
    [6] R. Brown, E. Curry, L. Magnani, C. S. Wilhelm-Benartzi, J. Borley, Poised epigenetic states and acquired drug resistance in cancer, Nat. Rev. Cancer, 14 (2014), 747-753. doi: 10.1038/nrc3819
    [7] O. S. Rukhlenko, F. Khorsand, A. Krstic, J. Rozanc, L. G. Alexopoulos, N. Rauch, et al., Dissecting RAF inhibitor resistance by structure-based modeling reveals ways to overcome oncogenic RAS signaling, Cell Syst., 7 (2018), 161-179.
    [8] X. Sun, J. Bao, Y. Shao, Mathematical modeling of therapy-induced cancer drug resistance: connecting cancer mechanisms to population survival rates, Sci. Rep., 6 (2016), 22498. doi: 10.1038/srep22498
    [9] J. Zhang, F. Zhou, X. Wu, X. Zhang, Y. Chen, B. S. Zha, et al., Cellular pharmacokinetic mechanisms of adriamycin resistance and its modulation by 20(S)-ginsenoside Rh2 in MCF-7/Adr cells, Br. J. Pharmacol., 165 (2012), 120-134.
    [10] C. B. Gambacorti-Passerini, R. H. Gunby, R. Piazza, A. Galietta, R. Rostagno, L. Scapozza, Molecular mechanisms of resistance to imatinib in Philadelphia-chromosome-positive leukaemias, Lancet Oncol., 4 (2003), 75-85. doi: 10.1016/S1470-2045(03)00979-3
    [11] M. E. Gorre, M. Mohammed, K. Ellwood, N. Hsu, R. Paquette, P. N. Rao, et al., Clinical resistance to STI-571 cancer therapy caused by BCR-ABL gene mutation or amplification, Science, 293 (2001), 876-880.
    [12] F. McCormick, New-age drug meets resistance, Nature, 412 (2001), 281-282. doi: 10.1038/35085665
    [13] P. Bajger, M. Bodzioch, U. Foryś, Singularity of controls in a simple model of acquired ´ chemotherapy resistance, Discrete Contin. Dyn. -B, 24 (2019), 2039-2052.
    [14] J. Foo, F. Michor, Evolution of acquired resistance to anticancer therapy, J. Theor. Biol., 355 (2014), 10-20. doi: 10.1016/j.jtbi.2014.02.025
    [15] N. Kumar, G. M. Cramer, S. Dahaj, B. Sundaram, J. P. Celli, R. V. Kulkarni, Stochastic modeling of phenotypic switching and chemoresistance in cancer cell populations, Sci. Rep., 9 (2019), 10845. doi: 10.1038/s41598-019-54346-0
    [16] A. Hodgkinson, L. Le Cam, D. Trucu, O. Radulescu, Spatio-Genetic and phenotypic modelling elucidates resistance and re-sensitisation to treatment in heterogeneous melanoma, J. Theor. Biol., 466 (2019), 84-105. doi: 10.1016/j.jtbi.2018.11.037
    [17] J. Cosgrove, J. Butler, K. Alden, M. Read, V. Kumar, L. Cucurull-Sanchez, et al., Agent-based modeling in systems pharmacology, CPT: Pharmacometrics Syst. Pharmacol., 4 (2015), 615-629.
    [18] P. Sudalagunta, M. C. Silva, R. R. Canevarolo, R. R. Alugubelli, G. DeAvila, A. Tungesvik, et al., A pharmacodynamic model of clinical synergy in multiple myeloma, EBioMedicine, 54 (2020), 102716.
    [19] A. Kaznatcheev, J. Peacock, D. Basanta, A. Marusyk, J. G. Scott, Fibroblasts and alectinib switch the evolutionary games played by non-small cell lung cancer, Nat. Ecol. Evol., 3 (2019), 450-456. doi: 10.1038/s41559-018-0768-z
    [20] J. H. Goldie, A. J. Coldman, A mathematic model for relating the drug sensitivity of tumors to their spontaneous mutation rate, Cancer Treat. Rep., 63 (1979), 1727-1733.
    [21] A. J. Coldman, J. H. Goldie, A stochastic model for the origin and treatment of tumors containing drug-resistant cells, Bull. Math. Biol., 48 (1986), 279-292. doi: 10.1016/S0092-8240(86)90028-5
    [22] J. H. Goldie, A. J. Coldman, Drug Resistance in Cancer: Mechanisms and Models, Cambridge University Press, New York, 2009.
    [23] N. L. Komarova, D. Wodarz, Drug resistance in cancer: Principles of emergence and prevention, Proc. Natl. Acad. Sci., 102 (2005), 9714-9719. doi: 10.1073/pnas.0501870102
    [24] N. Komarova, Stochastic modeling of drug resistance in cancer, J. Theor. Biol., 239 (2006), 351- 366. doi: 10.1016/j.jtbi.2005.08.003
    [25] Y. Iwasa, M. Nowak, F. Michor, Evolution of resistance during clonal expansion, Genetics, 172 (2006), 2557-2566. doi: 10.1534/genetics.105.049791
    [26] J. Foo, F. Michor, Evolution of resistance to targeted anticancer therapies during continuous and pulsed administration strategies, PLoS Comput. Biol., 5 (2009), e1000557.
    [27] I. Bozic, J. G. Reiter, B. Allen, T. Antal, K. Chatterjee, P. Shah et al., Evolutionary dynamics of cancer in response to targeted combination therapy, ELife, 2 (2013), e00747.
    [28] C. Tomasetti, D. Levy, An elementary approach to modeling drug resistance in cancer, Math. Biosci. Eng., 7 (2010), 905-918. doi: 10.3934/mbe.2010.7.905
    [29] E. Afenya, C. Calderón, Diverse ideas on the growth kinetics of disseminated cancer cells, Bull. Math. Biol., 62 (2000), 527-542.
    [30] J. Gallaher, A. Babu, S. Plevritis, A. R. A. Anderson, Bridging population and tissue scale tumor dynamics: a new paradigm for understanding differences in tumor growth and metastatic disease, Cancer Res., 74 (2014), 426-435. doi: 10.1158/0008-5472.CAN-13-0759
    [31] A. Obenauf, Y. Zou, A. L. Ji, S. Vanharanta, W. Shu, H. Shi, et al., Therapy-induced tumour secretomes promote resistance and tumour progression, Nature, 520 (2015), 368-372.
  • MBE-18-01-018-supplementary.pdf
  • This article has been cited by:

    1. Zakir Hussain, Malaya Dutta Borah, 2022, A Computational Aspect to Analyse Impact of Nutritional Status on Drug Resistance, 978-1-6654-7100-8, 1, 10.1109/SILCON55242.2022.10028912
    2. A. Sa’adah, D. A. Kamil, G. E. Setyowisnu, 2022, 2501, 0094-243X, 020004, 10.1063/5.0091002
    3. Zhihao Zhong, Chao Liu, Yatao Xu, Weili Si, Wenjun Wang, Liping Zhong, Yongxiang Zhao, Xiaochen Dong, γ ‐Fe 2 O 3 Loading Mitoxantrone and Glucose Oxidase for pH‐Responsive Chemo/Chemodynamic/Photothermal Synergistic Cancer Therapy , 2022, 11, 2192-2640, 2102632, 10.1002/adhm.202102632
    4. Hao Jiang, Jingxin Liu, You Song, Jinzhi Lei, Quantitative Modeling of Stemness in Single-Cell RNA Sequencing Data: A Nonlinear One-Class Support Vector Machine Method, 2024, 31, 1557-8666, 41, 10.1089/cmb.2022.0484
    5. Sachin G. Nair, Sonu Benny, Wesley M. Jose, T.P. Aneesh, Epigenetics as a strategic intervention for early diagnosis and combatting glycolyis-induced chemoresistance in gynecologic cancers, 2024, 358, 00243205, 123167, 10.1016/j.lfs.2024.123167
    6. Evgenii Khailov, Ellina Grigorieva, Optimal Melanoma Treatment Protocols for a Bilinear Control Model, 2023, 11, 2227-7390, 3289, 10.3390/math11153289
    7. Prathibha Ambegoda-Liyanage, Sophia R.-J. Jang, Resistance in oncolytic viral therapy for solid tumors, 2024, 469, 00963003, 128546, 10.1016/j.amc.2024.128546
    8. Muse Ji, Hongbing Liu, Xinxin Liang, Mingli Wei, Dongmei Shi, Jingxin Gou, Tian Yin, Haibing He, Xing Tang, Yu Zhang, Tumor cells are under siege from all sides: Tumor Cell-Mimic Metal−Organic framework nanoparticles triggering Cuproptosis/Ferroptosis/Apoptosis for Chemo-Chemodynamic-Photothermal-Immunological synergistic antitumor therapy, 2024, 485, 13858947, 149640, 10.1016/j.cej.2024.149640
    9. Prathibha Ambegoda, Hsiu-Chuan Wei, Sophia R-J Jang, The role of immune cells in resistance to oncolytic viral therapy, 2024, 21, 1551-0018, 5900, 10.3934/mbe.2024261
    10. Cordelia McGehee, Yoichiro Mori, A mathematical framework for comparison of intermittent versus continuous adaptive chemotherapy dosing in cancer, 2024, 10, 2056-7189, 10.1038/s41540-024-00461-2
  • Reader Comments
  • © 2021 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5575) PDF downloads(518) Cited by(10)

Figures and Tables

Figures(6)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog