D (2010) An Elementary Approach to Modeling Drug Resistance

Resistance to drugs has been an ongoing obstacle to a successful treatment of many diseases. In this work we consider the problem of drug resistance in cancer, focusing on random genetic point mutations. Most previous works on mathematical models of such drug resistance have been based on stochastic methods. In contrast, our approach is based on an elementary, compartmental system of ordinary differential equations. We use our very simple approach to derive results on drug resistance that are comparable to those that were previously obtained using much more complex mathematical techniques. The simplicity of our model allows us to obtain analytic results for resistance to any number of drugs. In particular, we show that the amount of resistance generated before the start of the treatment, and present at some given time afterward, always depends on the turnover rate, no matter how many drugs are simultaneously used in the treatment.


Introduction
One of the main reasons for the failure of chemotherapy is the development of drug resistance.Drug resistance strongly limits the success of the therapy and reduces the survival rates.It is a fundamental problem since cancer that becomes resistant to all available drugs, may leave the patient with no therapeutic alternatives.
There are multiple mechanisms by which drug resistance may develop.Drug resistance may be permanent or temporary (in which case it can reverse over time).Resistance is either relative or absolute, where relative resistance refers to a resistance that depends on the dosage.Typically, in relative resistance, the higher the dosage is, the smaller is the resistance that develops.On the other hand, in absolute resistance no matter what the dose is, a resistant cell will not be affected by the drug.In general, contemporary literature seems to consider drug resistance to be mainly relative rather than absolute (see, e.g., [6,15,28,33]).
Drug resistance appears to be both a spontaneous phenomenon caused by random genetic mutations, as well as an induced one (which means that patients that take a particular drug increase the chances of developing resistance to it).Indeed, there is clear experimental evidence that at least in the case of some drugs, known as "mutagenic drugs", the drug itself can induce resistance to itself (e.g., see [34,35,37]).Such resistance is referred to as "drug-induced resistance".
There are various classifications of drug resistance, which we briefly summarize below.For a more comprehensive picture we refer to the book by Teicher [38] and to the references therein.
One type of resistance is "kinetic resistance".This term refers to the reduction in the effectiveness of the drug which is caused by the cell division cycle.Such resistance is generally only temporary.Many drugs (such as methotrexate, vincristine, and citosine arabinoside, to name a few) are mainly effective during only one specific phase of the cell cycle, e.g., during the S phase, when the DNA is synthesized.Thus, in the case of a short exposure to the drug, the cell will not be affected if during that time it is in a different phase.Even more importantly, the cell will be substantially invulnerable if it is out of the cell division cycle, i.e., in a "resting state" or in the G 0 state.This means that the number of cells that are affected by the drug, is lower for cell populations that have low proliferation rates.
Resistance to drugs may develop as a consequence of genetic events such as mutations, rather than developing due to kinetic reasons.This category includes both "point mutations" and "chromosomal mutations", also known as "gene amplifications".Point mutations are random genetic changes that occur during cell division.These mutations cause the replacement of a single base nucleotide or pair, with another nucleotide or pair in the DNA or RNA.This is a random event with a very small (yet nonzero) probability that modifies the cellular phenotype, making any of its daughter cells resistant to the drug.Gene amplification is the consequence of an overproduction of a particular gene or genes.This means that a limited portion of the genome is reproduced to a much greater extent than the replication of DNA composing the remainder of the genome.Such a defect amplifies the phenotype that the gene confers on the cell, which, in turn, induces resistance by essentially providing the cells with more copies of a particular gene than the drug is able to cope with.While gene amplification (and consequently the resistance induced by it) may be a temporary phenomenon, point mutations are permanent.
Drug resistance has been extensively studied in the mathematical literature.We would like to comment on some of the works and refer interested readers to the references therein.Kinetic resistance has been mathematically studied in [1] and [32].The models in these papers are based on ODEs.An alternative approach, using age-structured models can be found in [2,5,7,8,41].The first model of resistance due to point mutations is the celebrated model by Goldie and Coldman (see [3,4,9,10,11,12,13]). Using stochastic processes, Goldie and Coldman show how the probability of having no drug resistance is larger in smaller tumors.By considering the effects of various treatment regimes another conclusion of their study is that the best strategy to administer drugs is to use all available drugs simultaneously.A more recent work on point mutations is by Komarova et al. [20,21,22,23,24].For example, in [20,21], probabilistic methods and a hyperbolic PDE are used to show how the pre-treatment phase is more significant in the development of resistance than the treatment phase.Another recent work on point mutations is by Iwasa et al., [18], in which continuous-time branching processes are used to calculate the probability of resistance at the time of detection of the cancer.Modeling of resistance due to gene amplification can be found, e.g., in [16,17,19].Drug resistance in these works is studied using stochastic branching models.Finally, for mathematical models and experimental findings on drug induced resistance, we refer, e.g., to [14,37].
In this paper we introduce a new and very simple model for absolute drug resistance, which we assume is caused only by random genetic point mutations.For this particular setup, most existing mathematical models are based on stochastic methods.In contrast, our approach is based on a compartmental system of ODEs, whose variables are the normal cancer cell population (that is susceptible to the drug), and the population of cancer cells that are resistant to the drug due to point mutations.The purpose of this paper is to show that elementary ODE-based techniques can be successfully used to obtain comparable results to those that were previously derived using the substantially more complex mathematical machinery of stochastic methods.
The structure of the paper is as follows: In Section 2, we outline the basic modeling assumptions and develop our mathematical model of drug resistance.We start by considering the single drug case, and proceed with the more general case of two or more drugs.The main results that we obtain from the model are presented in Section 3. We show that the amount of resistance generated before the beginning of the treatment, and present at some given time afterward, always depends on the turnover rate, no matter how many drugs are used simultaneously.We also use our model to compare the amount of resistance that originates before and after the treatment starts.Section 4 focuses on studying the differences between our results and other works in the field.Concluding remarks are provided in Section 5.
Once again we would like to stress that the main contribution of this work is in providing an elementary way to derive comparable results to those that were previously obtained using substantially more complex mathematical machinery.

An Elementary Model for Drug Resistance
In this section we develop a simple mathematical model for absolute drug resistance in the presence of random genetic point mutations.Our approach is to model the process using a linear system of ODEs.In its essence, our model enjoys similarities with the well known model of Goldie and Coldman [11] (see also [7]).The main difference between our approach and [11], is that we do not make any use of probability theory.Instead, we rely on a purely deterministic system.The advantages and disadvantages of using such an approach will be examined below.

The single-drug case
We start with the case of resistance to a single drug.Accordingly, we follow two populations: The first group is composed of wild-type cancer cells (cells that are sensitive to the drug).We denote the number of wild-type cancer cells at time t, by N (t).The second group are cells that have undergone a mutation.These cells are resistant to the drug.The number of mutated cells at time t is denoted by R(t).
We assume that cancer grows exponentially according to the Skipper-Schabel-Wilcox model, also known as the log kill model (see [36]).We also assume that the drug therapy starts at t * .Our model can then be written as: (2.1) The system (2.1) describes the pre-treatment phase, while the system (2.2) follows the dynamics after the treatment starts.The difference between both systems is the introduction of H, the drug-induced death rate, a term that appears only after the treatment starts, i.e., in (2.2).
In both systems, L, D, and u denote the birth, death, and mutation rates, respectively.We assume that 0 ≤ D < L and 0 < u 1.The system (2.1) is written assuming that mutations occur as a result of a wild-type cell differentiating into one wild-type and one mutant cell.This is a standard assumption, see, e.g., [21].
The initial conditions for the pre-treatment system (2.1) are given as constants N (0) = N 0 = 0 and R(0) = 0.The initial conditions for the system (2.2) are N (t * ) and R(t * ), which are the solutions of (2.1) at t = t * .

Remarks:
1.In this model we assume that both the wild-type and the resistant (mutated) cells have the same birth and death rates.This assumption is made in order to simplify the initial presentation, and can be easily modified to model situations where the resistant cancer cells R have a relative fitness advantage or disadvantage with respect to the wild-type cancer cells N .
2. Another modification could be to replace the exponential growth of cancer by a different model such as the gompertzian growth, which was empirically shown to provide a better fit (see the Norton-Simon hypothesis in [29], [30], [31]).
3. We assume that mutations happen only in one direction, i.e., wild-type cells mutate and become resistant but not vice versa.This seems to be a reasonable assumption in the case in which the focus is on point mutation resistance and not on resistance caused by gene amplification.Indeed, the probability of reversal of a point mutation is much smaller than the probability of the point mutation itself, and can therefore be neglected.
4. The time of the beginning of the treatment, t * , can be related to the size of the tumor at that time.If we assume that the total number of cancer cells at time t * is M , we can use the exponential growth of cancer and the relatively small mutation rate u to estimate t * as

The 2-drug case
We now consider the case of a treatment in which two drugs are being simultaneously used.We denote by R 1 (t) and R 2 (t), the mutant cancer cell populations that mutated by time t so that they are resistant only to the first or to the second drug, respectively.We reserve the notation R(t) for the population of cells that are resistant to both drugs at time t.With these notations, the model for drug resistance with two drugs can be written as: Similarly to the single-drug case, also in the 2-drug case we distinguish between the pretreatment dynamics (described by (2.4)) and the dynamics after the treatment started, which is given by (2.5).

Remarks:
1. Note that we assume that if a cell is already resistant to one drug it does not make it any less vulnerable to the combination of drugs (i.e., H remains unchanged).Such an assumption can be justified given that for cells that are resistant only to one drug, the second drug is still effective.This assumption can be easily modified to account for matters of physical or chemical nature that can influence the level of resistance that a cell (that is already resistant to one drug) may develop to the second drug.
2. We assume that the probability of a point mutation, and consequently the mutation rate u, is the same for any non fully resistant state in which a cell may be, an assumption that can also be easily modified.
Extending the 2-drug model (2.4)-(2.5) to the setup of n drugs is straightforward.The resulting general n-drug case is given in Appendix A.

Analysis and Results
In this section we discuss the main results that can be obtained by analyzing the models for drug resistance that were introduced in Section 2.

Dependence on the turnover rate D/L
Our first result is that the amount of resistance present at the time when the treatment starts, t * , always depends on the turnover rate D/L no matter how many drugs are simultaneously used.This result is easily obtained by considering the solutions of R(t * ) in the single-drug systems (2.1), the two-drug system (2.4), and the n-drug system (A.1).For details on the derivation of the solutions we refer to Appendix B.
For the single drug case the solution is (using (2.3) to evaluate t * ) where M is the total number of cancer cells when the therapy begins.Similarly, for the two-drug therapy we have the following The extension to the general n-drug case is obvious: The various expressions for R(t * ), (3.1)-(3.3),contain the turnover ratio D/L.We can thus conclude that the slower the growth of the cancer is (i.e., the closer the turnover rate D/L is to 1) the larger is the pre-treatment drug resistance.Conversely, the faster the tumor grows (i.e., the closer the turnover rate is to zero) the smaller is the resistance that develops prior to the beginning of the treatment.This result is independent of the number of drugs.

How much resistance originates before the treatment?
Assume that mutations could be terminated after time t * so that the only drug resistance that is present after t * would be the "progeny" of the resistance generated before therapy started.
We refer to such resistance as the "pre-treatment resistance at time t" and denote it by R p (t).We would like to compare R p (t) with the resistance that is generated exclusively by mutations that occur during treatment, which we denote by R d (t), and refer to as the "during-treatment resistance at time t".Mathematically we can stop the mutations by setting u = 0. Our second result is that for any t > t * , the pre-treatment resistance is greater than the during-treatment resistance, i.e., R p (t) ≥ R d (t).This result holds under the assumptions that (L − D) < H and M/N 0 ≥ C, where the constant C depends on the number of drugs.For example we will see that in the one-drug case, C = e, and in the two-drug case C = exp(1 + (3)).
Indeed, in the single drug case, the solution of the system (2.2), subject to the initial conditions N(0) = M and R(0) = 0, is given by Here, to simplify the notations, time is measured from the beginning of the treatment, i.e., t = 0 refers to what we previously considered to be t = t * .On the other hand, R p (t) is the solution of (2.1) at time t * that is then multiplied by an exponential term e (L−D)t that accounts for the growth of this resistance during treatment, Thus, clearly R p (t) ≥ R d (t) for any t ≥ 0.Moreover, for sufficiently large t, we have which nicely illustrate the key players in determining the proportion between the two populations.We would like to stress that in the case of mutagenic drugs, where the mutation rates are much higher during treatment, the result may be easily reversed, i.e., the resistance generated after the beginning of the treatment may exceed the pre-treatment resistance.
For two drugs, the solution of the system (2.5), subject to the initial conditions Similarly to the single drug case, R p (t) is the solution of (2.4) multiplied by the exponential term e (L−D)t to account for the growth of such resistance during treatment, It is now easy to verify that R p (t) ≥ R d (t) for any t ≥ 0, unless t is large and which is impossible given the assumptions.

Remarks:
1.A similar result can be obtained in the n-drugs case, and the details can be found in Appendix A.
2. Equations (3.5) and (3.8) clearly show how the amount of resistance generated before the beginning of the treatment and present, including its progeny, at any given time afterward depends on the turnover rate.
3. We note that in the two-drug case, the interval in which the pre-treatment resistance is smaller than the during-treatment resistance, i.e., R p (t) < R d (t), is slightly larger than in the single drug case.
4. If we allow the drug-induced death rate to take different values, for example, H/2 for the populations R 1 and R 2 and H for the population N , (which models a situation in which only one of the two drugs can kill these cells), then the interval for which R p (t) < R d (t) can be made larger.For the particular choice of these values, the interval becomes

A discontinuous administration of drugs
Thanks to the simplicity of the model, it is possible to consider the dynamics of resistance also in the more realistic case of treatment regimes where the drugs are administered intermittently.In such a scenario, one interesting issue is to study the connection between the frequency of the administration of the drug (e.g. when the total dosage is fixed) and the development of resistance.Such a study would be particularly interesting if we assume a different model for the cancer growth.This setup is beyond the scope of this paper and it is left for future research.

Discussion
The discussion in this section is devoted to comparing the results obtained with our model and the results obtained in [18,20,21].The problem studied in all papers is similar, i.e., drug resistance caused by random mutations, though the mathematical techniques being used are different.Here we use a deterministic system of ODEs as opposed to the various stochastic methods used in [18,20,21].The specific goals of each paper are somewhat different.In [18] the probability to have developed resistance by the time of clinical detection of the cancer is calculated.Instead, in [20,21] the focus is on the probability to have resistance as t → ∞.Finally, in our paper we calculate the "average" amount of resistance developed at any given time before and after the beginning of the treatment.We would like to note that while in this work we deal with numbers of cells, in [18,20,21] the quantities of interest are probabilities.
Comparing our results with some of the main results of [18,20,21], the outcome is as follows: 1. Our results agree with [18], in the sense that both works clearly show the dependence of resistance on the turnover rate at the time of detection (see equation (3.5) and table 1 of [18, page 2563], also found below as equation (4.1)).Furthermore, our results agree with [18] also on the effect that the mutation rate and the detection size have on resistance.
2. Most of our results agree with the results of [20,21] which are summarized on page 352 of [21], except for one issue.In our work, we show that no matter how many drugs are used in the treatment, the amount of fully resistant mutants (generated before the beginning of the treatment and present, including their progeny, at some given time afterward) depends on the turnover rate.See (3.1)-(3.3)and Remark 2 in Section 3. In contrast, one of the main results of [20,21] was that this dependence holds only in the multi-drug setup, and not for the single drug case (see [21, page 360]).The reason for this difference is due to the fact that [21] studies the probability to have such resistance in the limit, as t → ∞.It is only at t = ∞ that the results of [21] show a lack of dependence of the resistance on the turnover rate.This result does not hold at any finite time.Furthermore, the conclusion in [21] that in the single drug case, the probability of treatment success does not depend on the turnover rate (see [21, page 352]), is related to the definition of a successful treatment as a complete extinction of the tumor as time becomes infinite.Different definitions of a successful treatment (such as allowing tumors not to exceed a certain size or simply considering finite times) will lead to a dependence on the turnover rate also in the single drug case.
3. The elementary mathematical tools that we used in formulating our model, enable us to analytically study the n-drug case.Of course, with an increasing number of drugs, any deterministic model that considers averaged quantities has a limited validity.For two drugs or more, the results obtained using the approach in [20,21] rely on numerical simulations.
4. Due to the simplicity of our model is it easy to consider the dynamics of resistance also in the more realistic case of treatment regimes where the drugs are administered intermittently.It is also straightforward to compare different therapy regimes.
5. Since our approach is deterministic, it has a definite disadvantage when compared with stochastic models as it does not provide the time-dynamics of the probability distribution of the generated resistance.At the same time, it is important to note that the stochastic approach may yield results of a limited value.Indeed, according to [18], the probability that resistance develops by the time of detection in the one-drug case is given by Here, M is the detection size.Hence, for values of M = 10 9 (which, at present, is approximately the lower limit of clinical detection, see [33, page 31]) and u ≥ 10 −8 [25,27,40], the probability that resistance develops by the time a tumor is detected is always greater than .9999.This is the case since (L/D) ln(L/(L − D)) is always greater than 1 if D < L, where the last inequality is supported by experimental data [39].Hence, in this case the actual calculation of the probability distribution of the generated resistance, while showing the explicit dependence on the parameters, may not provide useful information if the product of M and u is greater than 1.
6.While from a mathematical point of view, it is a common practice to compute asymptotics as t → ∞, it is more desirable in the problem of drug resistance (and its related concept of treatment success) to study the dynamics for finite time, a time that is at most of the order of several years.

Conclusions
In this paper we developed a simple compartmental ODE model for absolute drug resistance, which is assumed to be caused only by random genetic point mutations.Our model was derived and analyzed for an arbitrary number of drugs.The simplicity of our model is the main strength of our approach.Alternative approaches to the problem that are based on stochastic methods, become very complex with an increasing number of drugs.Of course, with an increasing number of drugs, any deterministic model that considers averaged quantities has a limited validity.One of our results is that the amount of resistance that is generated before the beginning of the treatment, and which is present at some given time afterward, always depends on the turnover rate, no matter how many drugs are used.This result seems to contradict some results in the literature, as discussed in Section 4.
Additional results that were obtained with our model, validate the corresponding results in the literature.In particular, we demonstrate that at any given point in time, the resistance whose origin is before treatment is greater than the resistance created after the therapy starts.
where I is the identity matrix.The last matrix is nilpotent of order 2, thus the fundamental matrix is given by e (L−D)t 1 0 ut 1 .
From this (3.1)follows.Similarly, by writing the ODE system (2.4) in matrix form, we obtain the following matrix: The last matrix is nilpotent of order 3, thus the fundamental matrix is given by e (L−D)t which (3.2) follows.