|Home | About | Journals | Submit | Contact Us | Français|
We review a large volume of literature concerning mathematical models of cancer therapy, oriented towards optimization of treatment protocols. The review, although partly idiosyncratic, covers such major areas of therapy optimization as phase-specific chemotherapy, antiangiogenic therapy and therapy under drug resistance. We start from early cell-cycle progression models, very simple but admitting explicit mathematical solutions, based on methods of control theory. We continue with more complex models involving evolution of drug resistance and pharmacokinetic and pharmacodynamic effects. Then, we consider two more recent areas: angiogenesis of tumors and molecular signaling within and among cells. We discuss biological background and mathematical techniques of this field, which has a large although only partly realized potential for contributing to cancer treatment.
Mathematical modeling of anticancer chemotherapy has more than four decades of history (see e.g., Eisen (1979) and Swan (1990) for the survey of early research in this area). It contributed to the development of ideas of chemotherapy scheduling, multidrug protocols, recruitment and synchronization of cancer and normal cells. It also helped in refinement of mathematical tools of control theory applied to the dynamics of cell populations. However, few of these works found way into clinical application. The reasons for this failure are that, on one hand, important biological processes were unrecognized and crucial parameters were not known, and on the other, that the mathematical intricacy of the models was not understood. In this paper, we review the existing literature and discuss several topics, which may help improve the situation and realize the potential existing in the mathematical approach:
The philosophy of this paper is related to our professional experience. The authors have been involved for two decades in attempts to develop a satisfactory theory of optimal control of models of chemotherapy and antiangiogenic therapy using ordinary differential equations. This included collaboration with cancer researchers in the areas of data processing and mathematical modeling of results of DNA microarray experiments. Another focus was development of models of the cell cycle for the purpose of estimation of cell-cycle-phase specific action of anticancer drugs, as well as investigation of gene amplification as the mechanism of resistance of cancer cells. Still another focus was modeling and identification of signaling and regulatory pathways in cell-gene networks in collaboration with different biologists and clinicians. This latter approach currently seems to be on the rise.
Table 1 includes an idiosyncratic selection of models discussed in the present review in reverse chronological order. The purpose of this table is to provide an idea of predominant trends in modeling, without a pretense of completeness. Many more models are found in the body of the paper.
Carcinogenesis may be viewed as a result of loss of control over the cell cycle leading to an altered behavior of a part of cell population. Any treatment directed against this process may be viewed as a control action for which the cell cycle is an object. Control models are formulated to enable analysis and optimization of the diverse dynamic protocols of drug administration. In that, they differ from models which merely include therapy through changes in kinetic parameters, even though some of these latter are successful in matching experimental data (e.g. Michor et al., 2005; Komarova and Wodarz, 2005; Kim et al., 2008c). The control system approach, particularly when it is combined with pharmacokinetics and pharmacodynamics (PK and PD) models, allows finding an optimal schedule and dose of an anticancer drug.
The cell cycle is a sequence of phases traversed by each cell from its birth to division. These phases are: G1, or growth; S, or the DNA synthesis; G2, or preparation for division (gap phase); and M, or division. After division, the two progeny cells usually re-enter G1. However one or both progeny may become dormant or resting, or in other words, enter the quiescent G0 phase. From G0, after a variable and usually long time, the cells may re-enter G1 (Baserga, 1985).
This idealized description is a simplification. For example, in solid tumors there exists a geometric gradient of oxygen and nutrients. This causes stratification in viability of cells. Usually, cycling cells are located near the surface or near blood vessels; more distant layers are occupied by dormant cells, while the remotest regions form a necrotic core. This may lead to self-limiting growth phenomena, which may be described by biologically based nonlinear models including Pearl-Verhulst, Cox-Woodbury-Meyers, or Michaelis-Menten equations (see e.g. Murray, 1990; Swan,1985, 1986; Swan and Vincent, 1977) or Gompertz equations (Gyllenberg and Webb, 1989; Speer et al., 1984; Sullivan and Salmon, 1972; Wheldon, 1988; Zietz and Nicolini, 1979). It may be noted that Gompertz (1825) proposed his model for demographic purposes, so that the biological meaning of its parameters is not immediate. However, growth data from experimental tumors can be approximated by the Gompertz curve, as first suggested by Laird (1964). Although such models have been used to model the growth of cancer cell population and after some modifications allow including therapy scheduling, they do not include cell cycle structure.
One of the features, which are particularly interesting for modeling of chemotherapy scheduling, is cell-cycle-specificity of majority of anticancer agents. This feature is essential for the initial period of chemotherapy, when at issue is the most efficient reduction of the cancer burden. This is of practical importance in nonsurgical cancers such as leukemias but may be also of interest in solid tumors. Cell-cycle specificity of cytotoxic drugs is important since it makes sense to apply anticancer drugs when cells gather in the sensitive phases of the cell cycle. It can be approached by dissection of the cell cycle into disjoint compartments, with drug focused in a subset.
One class of mathematical models based on cell-cycle kinetics was introduced by us (Swierniak, 1989; Swierniak and Kimmel, 1984) and has been analyzed in numerous papers since (e.g. Swierniak and Polanski, 1994, Swierniak, 1995, Swierniak and Duda, 1994, Kimmel and Swierniak, 2006, Swierniak et al., 2003b), both from computational as well as analytical perspective.
In this approach, the cell cycle is modeled using compartments, which describe different cell cycle phases or clusters of phases. The transit times through the phases of the cell cycle are variable, particularly in malignant cells. Usually it is assumed that most of this variability is concentrated in the G1 phase (and in G0 whenever it exists). The simplest models arise if the transit times through each compartment are assumed exponentially distributed (Jansson, 1975; Takahashi, 1966, 1968). In this case we are led to the system of linear ordinary differential equations. The control of the cell cycle via the cytotoxic and cytostatic agents is introduced into these equations as a multiplicative term perturbing the parameters, which represent the average transfer times through different cell cycle phase. The resulting equations include linear combinations of products of state and control variables (see e.g. Mohler 1973). This leads to a number of analytical and computational problems. Similar models in discrete-time were used for optimization of radiotherapy protocols (Kim et al., 1974, Bahrami and Kim, 1975).
Although the control problem of cancer chemotherapy may be formulated for the entire cell cycle treated as a single compartment (e.g. Swierniak and Duda, 1992), solutions to these single- compartment models are oversimplified and not very informative. Of the more detailed multi-compartment models, the simplest and most natural ones are models with two and three compartments (Swierniak and Kimmel, 1984, Swierniak, 1995, Swierniak and Polanski, 1994, Swierniak et al., 2003a). In these models the G2 and M phases are combined into one compartment. In the two-compartment model, G0, G1 and S form the other compartment, while a range of three-compartment models arise by separately treating, either the DNA synthesis phase S or the dormant phase G0. These alternative models allow considering various drugs used in chemotherapy such as killing, blocking or recruiting agents. The killing agents are represented by G2/M specific agents, which include the so-called spindle poisons such as Vincristine, Vinblastine or Bleomycin ,which destroy the mitotic spindle (Calabresi and Schein, 1993) and Taxol (Fister and Panetta, 2000) or 5-Fluorouracil (Chabner and Logo, 1996) affecting cells mainly in division. Killing agents also include S-specific drugs such as Cyclophosphamide (Fister and Panetta, 2000) or Metotraxate (Clare et al., 2000) acting mainly in the DNA replication phase, Cytosine Arabinoside - Ara-C, rapidly killing cells in phase S through inhibition of DNA polymerase by competition with deoxycytosine triphosphate (Coly et al., 1984). Among the blocking drugs we can mention antibiotics Adriamycin, Daunomycin, Dexorubin and Idarudicin, which cause the progression blockage on the border between the G1 and S by interfering with formation of the polymerase complex or by hindering the separation of the two polynucleotide strands in the double helix (Alison and Sarraf, 1997). Another blocking agent is Hydroxyurea (Dibrov et al., 1983) which is found to synchronize cells by causing brief and invisible inhibition of DNA synthesis in the phase S and holding cells in G1. The recruitment action was demonstrated (Andreeff et al., 1992) for Granulocyte Colony Stimulating Factors - G-CSF, Granulocyte Macrophage Colony Stimulating Factors - GM-CSF, Interleukin-3 - IL-3, especially when combined with the Human Cloned Stem Cell Factor - SCF.
There is some controversy in the literature concerning the site and the mode of action of some drugs. For example, although mostly active in specific phases, Cyclophosphamide and 5-Fluorouracil kill cells also in other phases of the proliferation cycle. This makes it possible to combine them with other cycle specific agents (Bonadonna et al., 1995). On the other hand, some antimitotic agents such as curacin A (Kozusko et al., 2001) act by increasing the S phase transition and decreasing the M phase transition, which results in a G2 block.
Killing agents in our models, are assumed to affect the G2/M phase, which makes sense from a biological standpoint for several reasons. First, in mitosis, the cell wall becomes thin and porous and therefore the cell is more accessible. Second, chemotherapy during mitosis prevents the creation of progeny cells. While the killing agent is the only control considered in the two-compartment model, in the three-compartment models, also a blocking agent is considered which slows down the transit through S and then releases cells when another G2/M specific anticancer drug has a maximum killing potential (the so-called synchronization; Brown and Thompson, 1975). This strategy may have the additional advantage of protecting the normal cells which would be less exposed to the second agent, e.g., due to less dispersion and faster transit through G2/M; (Dibrov et al., 1983; Agur, 1988). This cell cycle model includes separate compartments for the G0/G1, S and G2/M phases. One of the problems in chemotherapy of leukemias are leukemic stem cells, which remain quiescent and are not sensitive to most cytotoxic agents (Michor, 2008; Komarova and Wodarz, 2007; Chabner and Longo, 1996). Similar findings were reported for breast and ovarian cancers, e.g. in (Fister and Panetta, 2000; Clare et al., 2000). As indicated by these authors, the insensitivity of dormant cells to the majority of anticancer drugs and high percentage of tumor cells at rest is a fact which, if ignored, leads not only to clinical problems but also to erroneous theoretical considerations. Experiments with Ara-C (Coly et al., 1984), indicated that while injected twice during cell cycle or combined with Adriamycin or anthracyclines, Ara-C led to serious reduction of leukemic burden without an evident increase of negative effect on normal tissues. This gain was attributed to the recruitment of leukemic cells from G0 by Ara-C. It also became possible to efficiently recruit quiescent cells into the cycle using cytokines (Tafuri and Andreef, 1990). Cytokines are substances playing a role in the regulation of normal hemopoiesis, such as like G-CSF, GM-CSF, and especially IL-3 combined with SCF. Then, a cytotoxic agent such as Ara-C or anthracyclines may be used. The three-compartment model that describes such therapy uses separate compartments for the G0, G1 and S+G2/M phases and includes such a recruiting agent. Moreover, it also enables analysis of the change of the transit time through the G0 phase due to the feedback mechanism that recruits the cells into the cycle when chemotherapy is applied. In a similar way we may model other types of manipulation of the cell cycle, such as for example the use of triterpenoids to inhibit proliferation and induce differentiation and apoptosis in leukemic cells (Konopleva et al., 2002).
In the models considered by us, the problem of finding an optimal cancer chemotherapy protocol is formulated as an optimal control problem over a finite time-interval, the fixed therapy horizon. The state variable is the average number of cancer cells and the control is the effect of the drug dosages on cell subpopulation. The goal is to maximize the number of cancer cells, which the agent kills, or minimize the number of cancer cells at the end of the therapy session, while keeping the toxicity to the normal tissues acceptable. The latter aspect is implicitly modeled by including an integral of the control over the therapy interval in the objective so that minimizing controls balances the amount of drugs given to the patient with the conflicting objective to kill cancer cells.
The idea was used in a mathematical model (Swierniak et al., 2003b) which has an arbitrary number of compartments. The models mentioned above as well as other compartmental models, whose dynamics arises from balance equations with constant transition rates, fall into this class. For example, drug resistance can be modeled using the same framework by adding additional compartments representing drug-resistant sub-populations of cancer cells. Analytical approaches to find optimal control in these models are based on application of the Pontryagin Maximum Principle (Pontryagin et al., 1962), which usually results in two types of strategies as candidates for optimality. The first one, called the bang-bang control, corresponds to a treatment protocol, which alternates maximum doses of chemotherapy with rest periods when no drug is administered. The second one, called the singular control, consists of application of varying doses at less than maximum level. Bang-bang controls, which are widely used as protocols in medical treatments, are the more natural choice as candidates for optimality, and it even has been observed numerically that singular protocols lead to a worse performance (Swierniak et al., 1996; Duda, 1994). Rigorous mathematical treatment of such models, using high-order necessary conditions for optimality, presented by Ledzewicz and Schaettler (2002a, 2002b, 2006) demonstrated that singular arcs were indeed excluded from optimality for these models. However, in general, the Maximum Principle only provides first order necessary conditions for optimality and therefore the trajectories it identifies may not be optimal. It is therefore important to further investigate optimality of these conditions. In Swierniak et al. (2003a), we formulated an algorithm, which allows to determine whether bang-bang controls, which satisfy the conditions of the Maximum Principle, are locally optimal or not.
One of the conceptual problems that are encountered concerns the "resonances". The resonance effect, postulated by several authors (e.g. Dibrov et al., 1985; Agur, 1988; Webb, 1990) as the way to either maximize the efficacy of treatment or to spare the organism’s normal cells, is similar to synchronization using blocking and killing agents, except that the resonance is usually mentioned in the context of a periodically administered single agent. At the resonance frequency and dosage, maximum efficiency is achieved. This tactics may be combined with an attempt to spare the normal cells by taking advantage of the difference in cell cycle duration of cancer and normal cells. Having this in mind, let us review several models using variants of synchronization and resonance. A similar effect is also observed in anticancer radiotherapy and is potentially of use in improving treatment protocols (Tarnawski et al., 2000).
Models in the form of ordinary differential equations ignore the fact that not all the distributions of transit times in the cell cycle are exponential. A more realistic description is using functional differential or integral equations. A good example is the model used to solve the simplest control problems arising in chemotherapy, concerning cells synchronization in normal tissues (Dibrov et al., 1983). The model of the cell cycle assumes shifted exponential or gamma distributions of the interdivision time. This leads to a functional differential equation with infinite delay. Since the purpose of drug administration is to arrest as many cells as possible in the specific phase, and not allow them to enter the phase in which they are exposed to the anticancer drug, the performance index should be maximized. This problem, which is mathematically more complicated, may be still solved using the abstract version of the Maximum Principle leading to bang-bang candidates for optimal chemotherapy strategies. Contrary to the optimization problems discussed previously, where side effects were considered only in the form of the integral term in the performance index, in this model the side effects of the drug are central. A similar line of reasoning can be applied to model toxicity to the bone marrow, which is one of the main limiting factors in chemotherapy. Such model was introduced by Panetta and analyzed as an optimal control problem by Fister and Panetta (2000). They proposed an objective quadratic in the control, which led to the protocols with gradually increasing doses after a no-dose interval, which seems to be nonrealistic. The use of an objective linear in control (Ledzewicz and Schaettler, 2007b) leads once more to bang-bang strategies, which remain structurally the same even if pharmacokinetics of the drug is taken into account. Another model used by Dibrov et al. (1985) includes both "probabilistic" and "deterministic" compartments. Probabilistic compartments are modeled by ordinary differential equations, while the cell density function in the deterministic phase is a solution of a partial differential equation of von Foerster (1959). Solving and introducing this solution to the equations for the probabilistic compartment leads to a system of equations with time-delay. Dibrov et al. (1985) examine the effect of multiple periodic treatments and find that the optimal period of drug administration almost coincides with the mean intermitotic time for normal cells and with this choice the fraction of normal cells in the sensitive phase at moment of drug injection is minimized.
A deterministic model which represents kinetic heterogeneity of cancer cells belongs to the class of age –structure models (Arino,1995; Bertuzzi et al., 1986; Arino et al., 2007; Webb, 1985). Such models allow a simple representation of cell populations with variable but uncorrelated cell cycle times. They exhibit the asynchronous exponential growth property, which means that population grows asymptotically exponentially growth with a steady distribution of the structure variable. More complex models, involving subpopulations with age structure, are required to take into account the different cell cycle phases which enables modeling of the effect of chemotherapy in the age structured context (see e.g. Lebowitz and Rubinow, 1974). Although this class of models has been developed by some authors to find recommendation for chemotherapy planning, it does not lead to general statements since no qualitative methods are known for analysis of such models. Simulation remains the only tool of their investigation. More promise could be related to the approach proposed by Shackney (1999) in which randomness occurs in the transitions between different cell maturation stages. Further developments of this idea enable to find an approximate formula for the distribution of cell cycle time as a function of parameters which may be dependent on the chemotherapy schedules (Arino et al., 2007). The results are however far from being applicable to some practical conclusions.
The problem of cytotoxity reduction by proper scheduling of phase specific drugs is also considered in a number of papers by Agur and coworkers (e.g. Agur, 1988; Agur et al., 1988; Cojocaru and Agur, 1992).Their line of reasoning is based on the so called Z-method and utilizes a much simpler model in the form of functional equation which may be considered as a linear difference equation with variable sampling time.
The elimination coefficient Z measuring the treatment efficacy is defined by the ratio of the elimination time of the malignant population and the one of critical host population. Agur et al. (1988) find that treatment efficacy is a nonmonotonous function of the relation between the cell generation time and the period of drug administration with the maxima occurring when the critical host cell cycle length is a multiple of the chemotherapeutic period. The results imply that short drug-pulses at appropriate intervals may be more efficient than a drug administered at arbitrary intervals or a continuously released drug. Under the condition that the cell cycle parameters of malignant cells have a relatively large variation, the drug protocol could be determined by the host temporal parameters alone and should reduce cytotoxity even in the case of similar mean cell cycle times for cancer and normal tissues.
A mathematically elegant model of resonance has been provided by Webb (1990) who used a distributed cell cycle model and the notion of Floquet exponents to demonstrate rigorously the existence of multiple resonances.
All these attempts, including our own, have to be viewed with caution, because of the existence of aftereffects in the action of many cytotoxic agents (Kimmel and Traganos, 1986). Actions of these drugs may extend beyond the span of a single cell cycle. For example, cells blocked in the S-phase of the cell cycle and then released from the block, may proceed apparently normally towards mitosis but then fail to divide, or divide but not be able to complete the subsequent round of DNA replication. If such effects are substantial, they are likely to disrupt or complicate the resonances.
Drug absorption, distribution, metabolism and excretion are regulated by molecular processes taking place inside all cells. It is now widely recognized that these regulatory processes can be significantly altered by mutations that occur in cancer cells, leading to drug resistance. Emergence of clones of cancer cells resistant to chemotherapy is important in treatment and prevention of systemic spread of disease. Models incorporate drug resistance by dividing cancer cell populations into separate compartments, representing subpopulations with different drug sensitivity. The number of compartments and the rate of flow among them depend on the assumptions about the mechanisms causing drug resistance.
Drug resistance can either be intrinsic, exhibited by a subset of cells from the very beginning of a treatment, or acquired, triggered by the chemotherapeutic agent (Baer, 2005). Both types of resistance can be taken into account in models, where they determine the number and form of equations as well as, in the case of intrinsic resistance, initial conditions.
A mutation can modify gene expression within a cell, which results in a changed cell behavior. Apart from regulation of DNA transcription, the best-studied form of regulation, the expression of a gene, may be controlled during RNA processing and transport, RNA translation, and the posttranslational modification of proteins. The degradation of proteins and intermediate RNA products can also be regulated. Proteins performing the regulatory functions mentioned above are produced by other genes. This gives rise to genetic regulatory systems structured by networks of regulatory interactions between DNA, RNA, proteins, and small molecules (de Jong, 2000). If a mutation disrupts such regulatory interaction at any stage, it can give rise to single-step drug resistance, as in earlier models (Coldman and Goldie, 1983, Panetta, 1998). On the other hand, if the so-called gene amplification process is involved, in which additional gene copies are acquired leading to an increase in transcription products, multicompartmental models describing populations with various levels of drug resistance are more appropriate (Kimmel et al.1998, Kimmel and Swierniak, 2006, Swierniak et al., 1999). Similar approach can be used if resistance to different drugs in one population is concerned (Michor et al., 2006; Iwasa et al., 2006).
Moreover, a mutation of a single gene can result in a resistance not only to a distinct cytotoxic agent, but also in cross-resistance to a whole range of drugs with different structures and cellular targets. This phenomenon is called multiple drug resistance (MDR) (Krishna and Mayer, 2000; Liscovitch and Lavie, 2002; Doherty and Michael, 2003; Ozben, 2006).
Taking into account the biological background presented in this section, it is clear that a mathematical model combining cancer chemotherapy with drug resistance should be based on the following assumptions:
These effects are typical for cancer cells because of their high mutability, which is much higher than in normal cells. Regardless of the particular molecules that play a role in acquired drug resistance, it can be assumed to begin with a mutation of a specific gene, triggered by a drug. Lists of those genes and drugs can be found in many papers, e.g. in (Takemura et al., 1997; Doherty and Michael, 2003). Depending on cell population and genes involved, this can lead either to classical irreversible mutation models (e.g., Goldie and Coldman, 1979) or models of reversible, multistage mechanisms such as gene amplification or other transformations of cancer cell gene expression (Kimmel and Axelrod, 1990; Kimmel et al., 1992; Axelrod et al., 1993). Both approaches seem applicable, depending on the biochemical mechanism of resistance and lead to new mathematical problems.
The emergence of resistance to chemotherapy has been first considered in a point mutation model of Coldman and Goldie (1983, 1986), Coldman et al. (1985), and then in the framework of gene amplification by Harnevo and Agur (1991, 1992, 1993). The main idea was that there existed spontaneous or induced mutations of cancer cells towards drug resistance and that the scheduling of treatment should anticipate these mutations. The consequences of the point mutation model can be translated into simple recommendations, which have even been recently tested in clinical trials.
One way the genome of cancer cells may rapidly evolve is by an increase in copy number of specific genes, referred to as gene amplification. Gene amplification can be enhanced by conditions that interfere with DNA synthesis and is increased in some mutant and tumor cells. Increased number of gene copies may produce an increased amount of gene products and, in tumor cells, confer resistance to chemotherapeutic drugs. Amplification of oncogenes has been observed in many human tumor cells and also may confer growth advantage on cells which overproduce the oncogene products (for an overview see surveys by Stark, 1993, and Windle and Wahl, 1992). Therefore, more realistic models should take into account these processes and, among other, predict the observed pattern of gradual loss of resistance in cancer cells in vitro placed in a non-toxic medium.
In the classical experiments of Schimke and his coworkers (Brown et al., 1981; Kaufmann et al., 1981), the anticancer drugs served to select for cells with amplified genes. In some of cell lines, when the selective agent was removed, the cells with amplified genes gradually disappeared from the population. It was observed that in such cases the amplified genes were located on extrachromosomal fragments of DNA called Double Minute Chromosomes (DM's). In these experiments, the dihydrofolate reductase (DHFR) gene was amplified after exposing murine 3T6 cells or mouse sarcoma S-180 cells to methotrexate.
In other cases, the amplification was stable, i.e., it persisted after the selective agent had been removed. In such cases, the amplified genes usually were located on elongated chromosome arms. The most regular of these elongated arms exhibited a regular band structure (the so called Homogeneously Staining Regions or HSR's), but other less regular structures also were observed. They were either caused by reintegration of extrachromosomal genes as proposed by Windle et al. (1991) or they arose by a separate mechanism as proposed by Smith (1976) and Stark (1993).
The population distribution of numbers of gene copies per cell can be estimated by flow cytometry after staining gene products. In the experiments mentioned, two features of these distributions were notable. (1) As expected, under nonselective conditions, the proportions of resistant cells (with amplified genes) decreased with time. (2) Less obvious, the shape of the distribution of the number of gene copies limited to the resistant cell subpopulation seemed to remain stable during the loss of resistance. A mathematical model of unstable drug resistance should take into account (i) stochastic changes in number of gene copies from one generation to another and (ii) the stochastic variability in cell lifetimes. A stochastic process which accommodates both (i) and (ii) is a random walk superimposed on the time-continuous branching process of cell proliferation, i.e., a branching random walk (Kimmel and Axelrod, 2002). A population is considered as a set of abstract particles of different types.
The simplest models of gene amplification in (Kimmel and Stivers, 1994) assume the above process. Cells with 2 j−1 gene copies are said to belong to type j (with 0 gene copies, to type 0). The parameters b and d are the probabilities of gene amplification and deamplification, respectively.
There is no doubt that the parameters of spontaneous cancer cells populations existing in vivo in humans differ considerably from these of the test-tube ”transformed“ cells and from those of cells in the induced animal tumors. However, much information regarding cell-cycle-phase specificity of anticancer agents has been obtained using in vitro experimental models. The numerical values of the probabilities of gene amplification and deamplification can be estimated based on data from the Schimke’s experiments. The probabilities of deamplification are of the order of 0.10, while the probabilities of amplification are about 5 times lower. In consequence, the process is strongly subcritical. This implies that in the absence of selection, the amplified phenotype disappears from the population. It can be revived by rare primary events, such as amplification of extrachromosomal genes following a deletion of the target gene from the chromosome arm (see further on). Further models of unstable amplification are discussed in Kimmel and Stivers (1994).
Stable amplification of DHFR gene was observed in the experimental system of Windle et al. (1991) in a Chinese Hamster Ovary (CHO) cell line, which contained only a single DHFR gene. Cells were challenged by methotrexate. Amplified genes residing on extrachromosomal elements were observed in cell cultures 8–9 generations later, while predominantly chromosomally amplified genes were seen after about 30 generations (only these two time points were investigated). This can be interpreted as an indication that some extrachromosomal elements containing amplified gene copy numbers are eventually reintegrated into chromosomes. In the model devised to reproduce these observations (Kimmel et al., 1992), the basic indivisible unit which serves as the template for the production of additional gene copies is the amplicon, which contains at least one copy of the target gene. The size of such structures could range from submicroscopic to an entire arm of a chromosome and they may be circular or linear. The acentric (replicating) element (ARE) is an extrachromosomal molecular structure containing one or more amplicons but no centromere. A centromere is required for regular segregation to progeny cells. The reintegrated element (RE) is the ARE after it has reintegrated into a chromosome. In our studies, different cell types correspond to different levels of cells’ sensitivity to a cytotoxic agent. Cells of type 0 are sensitive to the agent, whereas the other types consist of resistant cells of increasing level of resistance (for example, cells with increasing count of DHFR or CAD gene copies per cell).
The gene amplification model was extensively simulated and analyzed, employing methods briefly desribed below. In our works we have presented a model of chemotherapy based on a stochastic approach to evolution of cancer cells. Asymptotic analysis of this model lead to understanding of its dynamics and explanation of stable and unstable behaviors, depending on model parameters (Kimmel et al., 1998; Swierniak et al., 1999). Application of optimization techniques, in turn, resulted in recommendations for optimized therapy (Smieja et al., 2001; Swierniak and Smieja, 2005).
In the simplest case of cell cycle control with drug resistance evolving, by gene amplification, the resistant cells are totally insensitive to drug’s action. As a consequence, the first equation of the system is bilinear, since the treatment (control) term multiplies the growth rate coefficient of sensitive cells, whereas the remaining equations describing resistant cells are linear with constant parameters. Such simplification enables qualitative analysis of the asymptotic behavior of the system and explanation of some apparently surprising phenomena. Using Laplace transforms and its inversions we explained why, although the process is subcritical (amplification less than deamplification) and the probability of the primary mutation event is very low, the resistant population may grow exponentially. More precisely, we proved that the population decays only if the average proliferation time (the inverse of the parameter of the exponential distribution) of the resistant subpopulation is long enough, compared to the difference between deamplification and amplification probabilities. This means that if proliferation rate is the only parameter affected by control, the resistant subpopulation may maintain itself in the subcritical case.
Moreover the special structure of the model enables its treatment as a closed loop system with positive feedback. This, in turn, allows finding of the minimal dose of the anticancer agent which while administered sufficiently long ensures decay of the tumor population. It is also very helpful in formulation of the modified model which may be useful for optimization of chemotherapy protocols in this time. This model has a form of an integro-differential equation. Using an abstract version of the Maximum Principle, we found candidates for optimal treatment strategies (Smieja et al., 2001). These candidates are bang-bang strategies (as explained earlier on) with periodic protocols as suboptimal solutions easy for numerical consideration. Refinements including more realistic assumptions, such as partial drug resistance of a part of drug resistant subpopulation, multidrug therapy and phase specificity of anticancer agents were published by Smieja and Swierniak (2003) and Smieja (2008). In fact, most patients undergo the currently approved treatment regimens that may include multidrug chemotherapy (Kellen, 2003; Merino et al., 2004). There exist a large number of models of multidrug therapies (see e.g. Komarova and Wodarz, 2005). However, since they do not include explicit variables representing therapy, their use for developing strategies for optimal therapies requires further work.
Although a finite dimensional version of the models of drug resistance does not demonstrate all the phenomena discussed above it allows for rigorous mathematical treatment of optimal control problem resulting from optimization of drug protocols. Such analysis was performed in (Ledzewicz et al., 2004), leading to conclusion that while bang-bang strategies are usually optimal the optimality of singular controls (intermediate doses) depends on the proportion of resistant cells in the cancer cell population.
In the models described above a commonly made simplification is an identification of the drug dose with its concentration in plasma or even with its effects on cancer or normal cells. In reality these clearly are different phenomena and their relations are defined by the so called pharmacokinetics (PK) and pharmacodynamics (PD) (Collins and Dedrick, 1982). The first term describes the transition between the dose and the concentration in the plasma, while the second represents the effect of the given concentration on the cell viability. Roughly speaking PK describes what the cells do to the drug (e.g transportation, metabolization, withdrawal) and PD illustrates what the drug does to the cells (e.g. DNA damage, mitotic spindle destruction, myelosuppression). These phenomena are largely influenced by the diurnal rhythms of proteins and the mRNA levels of which vary according to a 24-hour period. They are also disturbed significantly by both acquired and genetically determined resistance of the cells.
The simplest model of PK was proposed by Bellman (1983) in the form of the first order linear differential equation (the first order inertia). It is equivalent to the assumption that for continuous infusion the growth and decay of the plasma concentration are exponential with the same time constant (Martin, 1992). It has been proved (Swierniak and Duda, 1992 for the simplest models, and Ledzewicz and Schaettler, 2005 for the general compartment models) that appending the compartmental models of the cell cycle by the linear PK compartment does not change qualitatively the results of optimization. The situation complicates when we want to take into account different rates in concentration build up and the drug clearance by the system. One way to overcome this difficulty is to relate the speed of the buildup with the dose of the drug. It leads to bilinear model of PK which formally is similar to the model of cell cycle kinetics (Ledzewicz and Schaettler, 2005) and by the same does not change the structure of the general multicompartmental model and in turn the properties of the solution to optimization problem are still in due. It is a first attempt of introducing nonlinearities into the PK model and it could be complicated by introducing more sophisticated relations. The model of pharmacodynamics (PD) describes the effectiveness of the drug. In some sense the models previously analyzed define the effectiveness by the fraction of ineffective cell divisions. Thus it changes between 0 and 1 and PD is modeled by a function linear in some range of concentration and equal to 0, and to 1, respectively outside this range. For low and high concentrations this model may be too simple. More realistic functions which describe the saturation effects include Michaelis–Menten formula or the sigmoid function. The former is reasonable for fast acting drugs which then saturate at the high concentrations while the latter more accurately approximate the effectiveness at both lower and higher concentrations. In multidrug treatment such formulations are possible only if the drug act in different way. If the drugs act similarly and in the same time their combined effectiveness depends on the concentrations of all of them and all their synergistic properties. Even for the most commonly used drugs these types of interaction are far from being recognized and understood. One way to avoid this problem is to bundle all such drugs into one control action in mathematical models. One important finding (Ledzewicz and Schaettler, 2005) is that geometric properties of PK/PD models have a direct effect on the qualitative solutions of optimization problems. For linear PK and linear or Michaelis-Menten PD models qualitatively results of optimizations do not change and quantitatively the results do not differ significantly. In the case of sigmoidal PD models nothing changes for convex regions but for the concave ones singular arcs (corresponding to intermediate doses of the drug) could be optimal.
Recently it can be observed an increasing interest in PK/PD modeling related to the mutual dependence between these effects from one side and resistance of cancer cells, especially multi drug resistance, drug-drug interactions including synergy/antagonism relations, and combined therapy from the other (Basdevant et al., 2005). Moreover above mentioned relationships between PK/PD mechanisms and the circadian clock as well as their influence on the micro-environmental phenomena in cancer development and response to therapies have been studied intensively (e.g. Levi and Schibler, 2007,Oklejewicz et al., 2008, Roskelley and Bissel, 2002, Panetta et al., 2008).
Angiogenesis is a complex process which leads to the formation of new vessels from existing ones and it is stimulated and controlled by molecular factors called activators (stimulators) and inhibitors (blockers) of angiogenesis. During progression of tumor these factors are released by the tumor to develop its vascular network which enables its growth and subsequently determines the possibility of cancer metastasis. Since vascularization is necessary for tumor development in late sixties of the last century a new anticancer therapy was proposed, which would target not directly the cancer cells, but the emerging vascular network. This so-called antiangiogenic therapy was for the first time hypothesized by Folkman (1971, 1972) more than thirty years ago. Folkman’s suggestions follow:
As it is known, the most important obstacle against successful chemotherapy is drug resistance acquired by cancer cells while the normal tissues retain sensitivity to the drug. This negative feature of chemotherapy may be used as an advantage in the antiangiogenic therapy which is directed towards a particular subset of normal tissues and only indirectly destroys tumor cells and it is why it has been called by Kerbel (1997) a therapy resistant to drug resistance. Antiangiogenic therapy does not exploit tumor cell sensitivity, relying instead on tumor suppression consequent to inhibition of associated vasculature. For more than ten years Folkman’s ideas were not followed by experimental or clinical investigations. However now, tumor angiogenesis belongs to the most inspiring areas of research in oncology. Kerbel (2000) presents 10 significant reasons for the explosive growth in tumor angiogenesis research and development of antiangiogenic drugs.
Complexity of the process of vascularization, as well as the way in which inhibitors, stimulators and antiangiogenic drugs act, results in complex models (see e.g. Mantzaris et al., 2004) applicable for simulation of the process but less useful in analysis or synthesis of therapy protocols. An exception are models proposed by Hahnfeldt et al. (1999) who suggested that the tumor growth with incorporated vascularization mechanism can be described by a Gompertz-type or logistic-type equation with variable carrying capacity, which defines the dynamics of the vascular network. Roughly speaking, these models incorporate the spatial aspects of the diffusion of factors that stimulate and inhibit angiogenesis into a non-spatial two-compartmental model for cancer cells and vascular endothelial cells. The models considered here belong to this class.
The important feature of the mathematical models of angiogenesis is the self-limiting growth phenomena. As mentioned in the Section 3.1, such properties are exhibited by Gompertz and logistic models in which the growth is limited by the carrying capacity. Hahnfeldt et al. (1999) proposed to define the Gompertz and logistic carrying capacity, which constrains the tumor growth, as a varying tumor volume sustainable by the vessels, approximately proportional to the vessel volume.
Although the equation proposed by Hahnfeldt to model the tumor growth appear similar to the Gompertz equation, the carrying capacity is not constant but varies with changes of the volume of the vessels. The dynamics of the growth of the vessel volume represented by its doubling time depend on the stimulators of angiogenesis, inhibitory factors secreted by tumor cells and natural death rate of the endothelial cells. In Hahnfeldt, et al. (1999) it has been assumed that the inverse of the doubling time is the sum of these three factors.
Considering factors influencing growth and decay of tumor vasculature, Hahnfeldt et al. (1999) asserted that tumor-driven inhibitors from all sites act systemically, whereas tumor-derived stimulators act locally. On the other hand analyzing a diffusion-consumption equation for the concentration of stimulator or inhibitor inside and outside the tumor, Hahnfeld et al. (1999) concluded that the inhibitor influences target endothelial cells in the tumor at a rate that grows ultimately as the area of the active surface between the tumor and the vascular network which in turn is propotional to the square of the tumor diameter. This leads to the conclusion that the inhibitory factor is proportional to the power 2/3 of tumor volume since the volume is proportional to the cube of the diameter. The modification of this model proposed by D'Onofrio and Gandolfi (2004) assumes that the effect of stimulators and natural mortality on the inverse of the doubling time is constant while the effect of inhibitors is proportional to the active surface of the area of tumor being in contact with the vascular network and therefore to the square of the tumor radius. Combinations of tumor growth models given by Gompertz-type and logistic-type equations with vascular network models proposed by Hahnfeldt et al. (1999) and D'Onofrio and Gandolfi (2004) result in four nonlinear models of tumor angiogenesis. Ergun et al. (2003) proposed yet another simplified model, : in which the growth of the vascular network is independent of the tumor size.
Nevertheless to obtain a complete model of the tumor growth in the vascular stage we should add one of the two proposed previously models of growth (Gompertz or logistic type). The interesting finding is that all these models have the same nontrivial equilibrium point. The models are strongly nonlinear but by a logarithmic change of variables and scaling transformations, it is possible to simplify them and find their asymptoptic properties using the standard Lyapunov type analysis of stability (local and global) following the line of reasoning presented in D'Onofrio and Gandolfi (2004).
Application of antiangiogenic therapy can be incorporated in the model using a factor multiplicatively increasing the death rate of the vessels. For the constant dose of the drug it is possible to find a dose such that the equilibrium point is equal to zero. According to the conditions of stability given in D'Onofrio and Gandolfi (2004) this leads to the conclusion that the vascular network and in turn the tumor can be eradicated. This conclusion is crucial for the analysis. It is enough to ensure that the population of endothelial cells responsible for angiogenesis behaves in a required way, because the size of tumor population tracks similar transients. In D'Onofrio and Gandolfi (2004) it has been proven that the same effect might be reached for periodic therapy with mean value satisfying an equality or nonequality condition, which is however generally only necessary. It is not sufficient, since the original Hahnfeldt model of eradication of the tumor depends on the shape of pulses in the periodic protocol. For some other models, this condition is necessary and sufficient.
Although during simulation all the models discussed lead to similar behavior if uncontrolled, their behavior in the presence of control, corresponding to different therapeutic protocols, may differ significantly. Moreover, clinical interpretation of the results also depends on the choice of the model.
Another class of models based on ordinary differential equations was proposed and analysed by Agur and co-workers (Agur et al., 2004; Arakelyan et al., 2002, 2005; Forys et al., 2005). The main purpose of these models is to reflect instability of the newly formed vessels structure. The models consist of a module of tumor dynamics and two others, the angiogenesis and vessel maturation modules, coupled through the action of regulatory proteins. The model simulations demonstrate the role of microenvironmental conditions in the metastatic potential of the tumor. The simplified version of the model enables analytical considerations, which reflect instability and cycles observed in the angiogenesis process. The drawback of this class of models is that they do not reproduce the stable behavior observed in less aggressive tumors.
Constant or periodic therapies discussed previously, which ensure tumor eradication, have an important drawback. They should be applied for a long therapy period. Duration of the treatment protocols and cumulated dose of the drugs should be limited because of shortage of the antiangiogenic drugs, their costs, and side effects. In Ergun, et al. (2003) and more rigorously in Ledzewicz and Schattler (2007a) the optimal control problem for the Ergun model with Gompertz type tumor growth and free terminal time is solved. The authors found that the optimal strategy consists of bang-bang and singular intervals. Swierniak (2009) proposed, for the model by d'Onofrio and Gandolfi (2004), to optimize the protocol in the fixed finite time of therapy with the primary goal of finding the treatment protocol that maximizes tumor cure probability. Taking into account the integral constraint on the dose of the drug, the optimization problem has isoperimetric form, which facilitates its solution. Moreover we may introduce a component responsible for the minimization of vascular network related to the tumor development.
The additional term related to the volume of vascular network may be regarded as a constraint imposed on the possible dynamics of the system. On the other hand by the choice of the weighting coefficients we obtain a new possibility of analysis of the mutual dependence between the tumor growth and the volume of the vascular network. Thus it is reasonable to provide an extensive analysis of their effect on the solution of the optimal control problem. Ledzewicz and Schattler (2007a, 2008) solved the optimal control problem for the Hahnfeldt original model in the similar way as for the Ergun model and once more proved that in the optimal strategy some elements are singular. However, a reformulation of the optimization problem for the Ergun model enables avoidance of singular arcs and leads to pure bang-bang solutions. The only model in which singular arcs as parts of optimal trajectory could not be eliminated is the standard Hahnfeldt model with Gompertz type tumor growth as it has been proved rigorously in Ledzewicz and Schattler (2008). Tumor angiogenesis belongs to the most inspiring areas of research in oncology. The most important constrain in efficient antiangiogenic therapy is the accessibility of antiangiogenic agents. This is why the rational anticancer therapy should contain combination of antiangiogenic therapy with more standard modalities of anticancer treatment such as radiotherapy. The effect of radiotherapy should be included in both compartments of the model because radiation destroys both cancer and normal tissues. The classical LQ model (e.g. Thames and Hendry, 1987) assumes that the damage to DNA, which is the principal target for the radiation, has two components: a linear one that is a consequence of a simultaneous break in both DNA strands caused by a single radiation particle, and a quadratic one that is the result of two separate but adjacent breaks in different strands caused by two different particles. In Ergun et al. (2003) the authors used this scheme combined with their model of vascular cancer growth which lead to a complicated optimal control problem. The resulting strategy contains both bang-bang and singular intervals. We discuss a model of such combined therapy (Swierniak, 2008) for the D'Onofrio-Gandolfi (2004) model. In our model we omit this second term and introduce only linear (in dose) effect into both equations. Accompanying radiotherapy supports the effect of antiangiogenic therapy. Moreover the effect of tumor eradication may be achieved easier and faster, although the theoretical results based on the theory of stability still have an asymptotic form only.
In Ergun, et al. (2003) the optimal control problem for the Ergun’s model was presented for a free terminal time. In our study, we proposed the same optimization scheme as discussed above for the antiangiogenic therapy. The integral constraints on control variables although similar in form, have different meaning. For radiotherapy they measure the feasible cumulated negative effect of the radiation, while for the antiangiogenic agent they mostly represent the limited availability of the agent and in part the potential side effects of the drugs. Once more, singular arcs for antiangiogenic therapy are not feasible since there are no finite intervals of constant solutions to the adjoin equation and singular arcs for radiotherapy do not satisfy higher-order necessary conditions of optimality. This leads to the conclusion that intermediate doses of the drug and radiation are not optimal and that the optimal protocol contains only switches between maximal dose and resting intervals.
Yet another possibility, arguably more feasible from clinical point of view, is to combine chemotherapy and antiangiogenic therapy. The first attempt in analysis and optimization of the model of such combined therapy was presented by Swierniak (2008). In this case however a two-compartment model seems nonadequate. The realistic model should take into account drug resistance of the cancer population treated by chemo-toxic agents and phase dependence of the drugs (see e.g. Kimmel and Swierniak 2006). The simplest model, which takes into account emergence of drug resistance in tumor chemotherapy and parallel treatment by antiangiogenic agents may be defined in the form of the three-compartment model which combines a very simple model of drug resistance (Ledzewicz et al., 2004) with one of the previously discussed models of antiangiogenic therapy. Although the class of models and optimization problems looks very similar to the ones previously analysed nevertheless we are led to the third order system of differential equations and therefore the synthesis of an optimal control law is much more difficult. Analysis of this problem will be the subject of our further studies. Yet another observation found by Ledzewicz and Schaettler (2007b) is that including PK effect by appending the Hahnfeldt model with first order linear dynamics equation does not change qualitatively the optimal treatment strategies.
Among the difficulties that arise in applying mathematically designed treatment protocols, two play a major role. First, the coefficients in the performance index for protocol optimization should be individually chosen taking into account patient’s condition. Second, the effects of a drug on cancer cell population vary, depending on the genotype and phenotype of patient’s cells. Both these factors call for including intracellular processes in the models. This is consistent with the current shift from general cytotoxic drugs towards targeted and personalized therapies. Although only few of these have been implemented, such as for example Gleevec many are under development (see references in Brehme et al., 2009). Theoretical studies are also under way, such as, for example, the model for chronic myeloid leukemia by Michor (2008).
However, so far, only several models of this type have been built. In most cases, the part of the model describing dynamics of intracellular processes is simplified in such way that it accommodates an a priori chosen treatment (Evans et al., 2004; Chappell et al., 2008; Kronik et al., 2008). Some approaches helped to suggest therapy protocols different from those actually in use (Kronik et al., 2008; Shochat and Rom-Kedar, 2008). In order to personalize treatment, more sophisticated models of intracellular processes should be applied and combined with optimization of treatment based on control theory (Cappuccio et al., 2007). An example of a sophisticated model of cell dynamics under treatment, which implicitly refers to within-cell processes, is the stochastic agent-based model of Roeder et al. (2006). The model describes the interaction between imatinib treatment and chronic myelogenous leukemia cells. The model accounts for the progression of normal and leukemic cells through three stages of the myeloid lineage: stem cells, progenitor cells, and mature cells. The agent-based formulation captures the inherent diversity of individual members of a large population and accounts for the probabilistic behavior of individual cells. The original algorithm is computationally demanding. However, Kim et al. (2008a) derived an equivalent system of difference equations that captures the dynamics of the original model which makes computations much more efficient and subsequently, a partial differential equation model Kim et al. (2008b).
One of the approaches to combat cancer involves understanding and amplifying organism’s own defenses. This can be achieved, for example, through intra-tumor administration of ex vivo activated alloreactive cytotoxic-T-lymphocytes, use of oncolytic viruses or treatment with cytokine-based drugs (Cohen and Kaufman, 2007; Moschos et al., 2007). Early work in this area includes construction of mathematical models, which reproduce dynamics of cancer cells targeted by oncolytic viruses (Paiva et al., 2008) and analysis of properties of these models (Wodarz and Komarowa, 2009). Other models involve injected CTLs (Kronik et al., 2008) or Interleukin-21-based therapy (Capuccio et al., 2006). Several models of signaling pathways playing important roles in cells’ responses to cytokines viral infections have been developed. For example, interferon (IFN) -β activated pathway has been studied by Smieja et al. (2008), and the pathway of IFN-γ by Zi et al. (2005). There exist several techniques of merging models describing processes on intracellular and population level (Haseltine et al., 2005 or Roose et al., 2007, Paiva et ai., 2009, Venkatasubramanian et al., 2008). A comprehensive but mathematically tractable model of interaction of dynamics of interferon-based defenses in a liquid cell culture was published by Getto et al. (2008).
Combining models at single-cell level with control-based approaches will help explain why some cancer cells escape chemotherapy while other do not, which may help making rational therapeutic decisions. Models involving within-cell processes, if they are mathematically or computationally tractable, have an increased chance of accomplishing this. We will outline steps needed to develop a mathematical model of the differentiation pathway, accounting for transcriptional and receptor-related noise. The framework is that in response to treatment, cancer cells and their environment release molecules (the „factors”), which activate signaling pathways leading either to programmed death (apoptosis) or to activation of defense mechanisms, such as differentiation or emergence of drug resistance. The mathematical form of the control and the specific model of resulting defense mechanisms depend on the mode of action of the drug. In order to understand the dynamics of the signaling pathway in response to the anticancer agent, it is necessary to construct a mathematical model linking the known proximal (cell membrane receptor-related) and distal (transcription factor-related) and usually unknown (intermediate) elements of the pathway, as well as the membrane receptor-related and transcriptional noise. We will describe the methodology developed in Lipniacki et al. (2007).
The first step is development of deterministic, ordinary differential equations (ODE) -based model, valid for averaged quantities, for the purpose of estimation of the coefficients of interaction in the signal transduction pathway. As mentioned above, the model has the following components: (i) cell membrane receptor-to-proximal kinase module, (ii) proximal kinase-to-transcription factor, and (iii) transcription factor-to-effector protein. The next step involves taking into account transcription noise and exploiting a stochastic differential equation model. In eukaryotic cells, most of the molecules such as mRNA and proteins occur in number of copies that is sufficiently large to justify a deterministic approximation. One process, in which stochasticity may occur, is the initiation of transcription, since the number of molecules of a transcription factor attached to its binding motif in a gene promoter, is either 0 or 1 (0,1, or 2 in diploid genome). Stochastic differential equations driven by such random process are known as Davies equations and have been known but rarely used, until they were introduced in the context of stochastic transcription (Lipniacki et al., 2006) A prediction of this model is that in individual cells, mRNA will be produced in bursts of random duration, causing fluctuations in pathway responses in individual cells (Raj et al., 2006). Both mRNA bursts and variability among single cells responses were experimentally observed.
Other stochastic factors include_receptor-related noise and Poisson process forcing. When the concentration of the factor becomes very low, not all the receptors on the cell membrane bind with the factor. Under such circumstances, there is a constant probability μ that a factor particle becomes bound to the receptor. Therefore, the number of saturated receptors at any time is a random variable with binomial distribution, which is approximated by a Poisson distribution when μ is small and N (the total number of receptors) is large. The parameter of the Poisson will be equal to λ = Nµ. In particular, a proportion of p0 =exp(−λ) cells will have 0 receptors saturated.
The model outlined above is unlikely to be analytically tractable, except in the simplest cases (see Paszek 2007). However, efficient simulation algorithms exists, which can model systems involving large-scale stochastic and deterministic components (such as Bertolusso and Kimmel 2009 and references therein). Examples of biological systems, where the modeling paradigm outlined above will be useful include the already mentioned Roeder (2006) model. Another one is the granulopoiesis process, which depends on granulocyte colony stimulating factor (G-CSF) signaling. For the granulocyte lineage, the most essential growth factor is G-CSF. Mice devoid of either G-CSF or its cognate receptor have profound neutropenia and subsequent susceptibility to infection (Lieschke et al., 1994; Liu et al., 1997). G-CSF initiates within the cell a signaling cascade, which leads to activation of a transcription factor called the G-CSF Receptor (G-CSFR), which is a member of the hematopoietin cytokine receptor superfamily Mutations in the G-CSF Receptor have been found in patients with Severe Congenital Neutropenia, Myelodysplastic Syndrome, and Acute Myeloid Leukemia (Dong et al., 1994, 1995a, 1995b). G-CSF is widely used to mobilize and harvest stem cells into the peripheral circulation for both stem cell transplantation (Vose and Armitage, 1995). At least five types of changes in the G-CSF Receptor have so far been found in patients with Myelodysplastic Syndrome/ Acute Myeloid Leukemia (Sultana et al., 2003).
One of the recent developments in modeling of chemotherapy, which take into account the internal signaling within cells, is a series of papers by Clairambault, Perthame, Levi and co-authors, exploring the importance of periodic change in cells’ status. Periodic oscillations in proliferating cells are caused by the internal oscillator governing the cell cycle. In a recent paper (Brikci et al., 2008), a nonlinear model of the dynamics of a cell population divided into proliferative and quiescent compartments has been presented. The proliferative phase represents the complete cell cycle (G1–S–G2–M) of a population committed to divide at its end. The model is structured by the time spent by a cell in the proliferative phase, and by the amount of Cyclin D/(CDK4 or 6) complexes. Cells can transit from one compartment to the other, following transition rules which differ according to the tissue state: healthy or cancerous. The asymptotic behavior of solutions of the nonlinear model is analyzed in two cases, exhibiting tissue homeostasis or tumor exponential growth. The model is simulated and its analytic predictions are confirmed numerically.
Another source of oscillations is the diurnal cycle of night and day, which is imposed on cells on higher organisms. This cycle should be taken into account when designing chemotherapy. (Clairambault et al., 2006). In another paper concerning this topic (Basdevant et al., 2005), the chronotherapy is defined as a regimen, which takes advantage of the circadian rhythm of cells physiology in maximizing a treatment efficacy on its target while minimizing its toxicity on healthy organs. In the paper, a mathematical model describing the time evolution of efficiency and toxicity of an oxaliplatin anti-tumor treatment has been derived. The authors then applied an optimal control technique to search for the best drug infusion laws. We will omit further details here, only noting that the authors prove that disregarding diurnal rhythms leads to less efficient treatments.
Recently, there has been a veritable explosion of models of tumor growth in 3 dimensions, involving mechano-chemical and mechano-biological considerations. Some of these models involve treatment components. These models, in principle, hold a great promise, although not much has transpired as yet in the terms of applications. We will very briefly and idiosyncratically review some of them including contributions, which refer to older but informative papers.
One of the earliest contributions was that of Chaplain et al. (2001) who examined the spatio-temporal pattern formation in reaction-diffusion systems on the surface of the unit sphere in 3D. A related mathematical problem has been considered by Bertuzzi et al. (2004), who considered the growth of the tumor cord, a cylindrical structure representing a prevascular stage of tumor growth. These models have been improved upon and developed in much detail, including chemotherapy, in further works (see papers by Sanga et al., 2006; Frieboes et al., 2009) . Many prevous publications are cited in these two papers. Modeling involving large number of factors interacting in spatial domain, is more realistic than previous, ideal mixing-based, computations. However, it is also clear that there is a large number of parameters involved, some of which can be deduced from chemistry considerations, but most of which are as yet unknown. Much work remains to be done to understand complicated patterns of tumor and metastasis, particularly under treatment. One interesting topic is the relationship between microscopic models, which refer to single cells, cell membrane receptors and growth factor molecules, and macroscopic models, which include continuous densities of these species. Marciniak-Czochra and Ptaschnyk (2008) devised a rigorous procedure based on homogenization, which helps understand this relationship.
Another interesting topic involves the consequences of the cancerization field theory. Some human tumors involve tubular or linear structures. In the lung, early lesions, the so-called Ground Glass Opacities, frequently have the form of intertwined linear shadows visible on Computed Tomography scans (as reviewed by Marciniak-Czochra and Kimmel (2007)). It seems that this type of growth may be present in the prevascular stage of lung cancer development. Since this type of growth does not produce nodular structures in Xray or CT images, its interpretation is ambiguous. It is an interesting question how such structures arise and what are their further invasion paths. Computational and analytical results (Marciniak-Czochra and Kimmel, 2008), show that reaction-diffusion systems are capable of producing quasi-chaotic patterns of tumor growth in space.
We are witnessing a revolution in the biological sciences, in which the multiscale systems biology approach begins to penetrate thinking at all levels, from molecules, to cells to organisms. This approach is an opposite of the typically reductionist approach, in which the whole is a simple sum of its parts. In systems biology, emergent properties, due to complex connections between the components, play a major role. This new thinking is due to a new appreciation of complexity and stochasticity in biology (Wilkinson, 2006). This in turn, is a logical consequence of new techniques such as DNA sequencing and measurements of gene and protein expression.
In this situation, mathematical modeling, which was transiently eclipsed by a deluge of new facts discovered by molecular biology and bioinformatics, starts experiencing a true renaissance. Models, which were developed two or three decades ago, such as some models of chemotherapy discussed in this paper, and then largely forgotten, now return to the fore. This happens simply because their structure can be informed by the new measurements and their parameters can be much better estimated.
Optimized cancer treatment protocols are likely to become in the future a standard element of personalized medicine (Faratian et al., 2009). When this happens, theoretical considerations such as those reviewed in this paper, will yield to engineering-type predictions.
This work has been partially supported by Polish Ministry for Science and Higher Education grants N N514 415334 in years 2007–2010 (JS), BK-218/RAu1/2009 (AS) and N N514 411936 in years 2009–2011 (MK) and by the NIH grant R01GM086885 (MK). The support has been used for the writing of the report and the decision to submit the paper for publication.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.