Some Systems of Ordinary Di erential Equations from Cancer Modeling: Qualitative Analysis and Optimal Treatment by Xiang Wan A thesis submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Master of Science Auburn, Alabama Aug 6, 2011 Keywords: cancer model, immunotherapy, radiotherapy, chemotherapy, bang-bang Copyright 2011 by Xiang Wan Approved by Georg Hetzer, Chair, Professor of Mathematics Wenxian Shen, Professor of Mathematics Xiaoying Han, Assistant Professor of Mathematics Abstract In this thesis, we study four systems of ordinary di erential equations, which model the interrelationships between di erent cell populations while tumor cells exist, and treatments, such as immunotherapy, chemotherapy, and radiotherapy, are applied. For the rst two models, we only consider the case with radiation treatment. For the rst model, we consider a single general cell population with its corresponding radiated cell population. Meanwhile, two di erent kinds of radiation are studies separately: constant and decay; for the second model, we consider the host and tumor cell populations together with their corresponding radiated cell populations, which behave in the same way as that in the rst model. For the third and fourth models, we consider both the immunotherapy and chemother- apy. The third model includes three cell populations: host cells, tumor cells, and immune cells, as well as the drug concentration. We study the properties of its null-surfaces, equilib- ria, and the stability of them; in the fourth model, we not only extend the previous model into one with six populations, with the immune cells in the third model being speci ed into three di erent ones: CD8+T cells, circulating lymphocytes, and IL-2. but also focus on the situation when controls are added in a linear manner. We investigate the existence of controls and nd the characterization of optimal bang-bang control. ii Acknowledgments This thesis would not have been possible without the guidance and the help of several individuals who in one way or another contributed and extended their valuable assistance in the preparation and completion of this study. First and foremost, my greatest thanks extend to Dr. Georg Hetzer for his patience and guidance throughout this endeavor at Auburn University. I got lots of advisable suggestions from him back to the time when I applied to this department, and Dr. Hetzer took consid- erable time and energy to further the progress in my studies of Mathematics. I will never forget the sincerity and encouragement that he showed to me. I would also like to thank Dr. Wenxian Shen on the committee, who is such a great teacher of mine for many courses. I have been learning so much from her. I am also thankful for the excellent example she has provided as a successful professor and a good person. My sincere thanks also go to Dr. Xiaoying Han for serving as the committee members. I wish to express my greatest appreciation for my parents, my sister, and her lovely daughter, who have given their love and constant support throughout my life. Thanks for always being with me. I o er my regards and blessings to all other people who supported me in any respect during the completion of the project. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 Background and Preliminary Knowledge . . . . . . . . . . . . . . . . . . . . . . 6 3 Model Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 General assumptions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2 Model I: Single population with Radiotherapy only . . . . . . . . . . . . . . 22 3.3 Model II: Double populations with Radiotherapy only . . . . . . . . . . . . . 23 3.4 Model III: Tumor model with Immune Resistance and Chemotherapy . . . . 24 3.5 Model IV: Immuno-Chemotherapy with controls . . . . . . . . . . . . . . . . 26 4 Analytic Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1 Model I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.1.1 Constant radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.1.2 Decay radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Model II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2.1 Equilibria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.3 Model III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 4.3.1 Null-surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.2 Equilibria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.3.3 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 4.4 Model IV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.4.1 Objective Functional . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 iv 4.4.2 Existence of Optimal Control . . . . . . . . . . . . . . . . . . . . . . 47 4.4.3 Characterization of the Optimal Bang-bang Control . . . . . . . . . . 50 5 Conclusion and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 A One-point compacti cation and associated results . . . . . . . . . . . . . . . . . 60 B Fundamental knowledge of optimal control theory . . . . . . . . . . . . . . . . . 61 v Chapter 1 Introduction Cancer has been known as a deadly disease of mankind from prehistoric time, and re- search into cancer tumor treatment began only recently from a historic point of view [15]. The growth of cancerous tumor is a complicated process involving multiple biological in- teractions. Also, the response of tumors to active treatments such as chemotherapy and radiotherapy is complex, but important to understand. A tumor?s response to treatment depends on many factors, including the severity of the disease, the application of the treat- ment, and the strength of patient?s own immune response. Mathematical modeling of this process is viewed as a potentially powerful tool in the development of improved treatment regimens. The mathematical modeling of tumor growth and treatment has been approached by a number of researchers using a variety of models over the past decades [33]. The main types of cancer treatments involve surgery, chemotherapy, radiotherapy, and immunotherapy, either in isolation, or in combination of two or more of these. A given speci c type of cancer will have a preferred treatment depending on, among other things, where the cancer is located and its stage of development. We introduce three therapies that are widely used in practice in the following context: Radiotherapy, Chemotherapy, and Immunotherapy. Radiotherapy Radiotherapy is a treatment procedure that uses radiation to kill malignant tumor cells. This treatment targets rapidly growing and dividing cells such as those in cancer[1]. Radiation destroys cells by causing one or more chromosomes to break. When this hap- pens, the cells cannot reproduce and eventually die o [2] [3] [4]. Hence the question of 1 persistence or extinction of a community of cells exposed to radiation is of paramount interest. Moreover, it is sometimes possible for the broken chromosomes to recombine. This may lead to the original con guration of the chromosomes, it may lead to muta- tion, or the recombination may be completely ine ective [4] [5]. Here we view the rst of these possibilities as a probability of broken cells becoming whole again. Finally, the radiation protocol may be one of two modes: constant dosage, which is similar to the case of long term radiation after a nuclear accident; decaying radiation, such as radioactive material implanted to ght lung cancer [6]. Chemotherapy Another of the most common and fundamental forms of treatment is chemotherapy, which involves injecting into the body a type of drug designed to attack the cancer cells. This type of drug also attacks normal (host) cells causing common side e ects such as hair loss (Radiotherapy has the same side e ects sometimes). Much research into chemotherapy is involved with designing the drug so as to maximize the e ect on cancer and minimize the side e ects [7]. Therefore, when developing e ective treat- ment strategies, understanding the e ects of chemotherapeutic drugs on tumors is of primary importance. Several approaches to modeling chemotherapeutic induced cell- kill (killing of tumor cells) have been developed. One of the early approaches was by citeskipper1964experimental, which propose that cell-kill due to a chemotherapeutic drug was proportional to the tumor population. This hypothesis is based on in vitro studies in the murine leukemia cell-line L1210. It states that for a xed dose, the reduction of large tumors occurred more rapidly than for smaller tumors. The concept in [8] is referred to as the log-kill mechanism. [9] [10] nd this model to be inconsistent with clinical observations of Hodgkin?s disease and acute lymphoblastic leukemia which showed that, in some cases, reduction in large tumors was slower than in histologically similar smaller tumors. Therefore, [9] [10] hypothesize that the cell-kill is proportional 2 to the growth rate (e.g., exponential, logistic, or Gompertz) of the tumor. A third hy- pothesis notes that some chemotherapeutic drugs must be metabolized by an enzyme before being activated. This reaction is saturable due to the xed amount of enzyme. Thus, [11] develops the Emax model which describes cell-kill in terms of a saturable function of Michaelis-Menton form. Immunotherapy Meanwhile, immunotherapy refers to the use of natural and synthetic substances to stimulate the immune response. This involves stimulating the immune system to work harder or using an outside source of cells, such as synthesized immune system proteins. Immunological therapies include the use of antigen and non-antigen speci c agents such as cytokines. Cytokines are hormones produced in the immune system that regulate the growth and activity of other immune system cells and blood cells. Cytokines alone can give the immune system a boost or given with other immunotherapies they can be used as adjuvants [12]. Cytokines have been used to treat melanoma, leukemia, lymphoma, neuroblastoma, Kaposis sarcoma, mesothelioma, brain cancer, cancer of the kidney, and cancer of the cervix. It is therefore important that we begin to develop mathematical models of tumor growth that include an immune system response, and ultimately a response to immunotherapy. Interleukin-2 (IL-2) is a cytokine that was approved by the FDA in 1992 for treatment of metastatic renal cell (kidney) cancer. IL-2 became the rst cytokine approved for use alone in treating advanced cancer [21]. Since that time, it has also been approved to treat people with metastatic melanoma. IL-2 can be used as a single-drug treatment for these cancers, or it may be combined with other forms of immunotherapy, such as vaccines. IL-2 helps immune system cells reproduce more rapidly once they are in the patient. The use of IL-2 together with chemotherapy or with other cytokines (such as 3 interferon-alpha) may increase their e ectiveness against some cancers, but the side e ects of the combined treatment are also increased [34]. Optimal Control Once an adequate model of interacting cell populations is constructed, we focus on the design of an improved treatment protocol. To this end, we employ the tools of optimal control theory. This theory originated in economics, where it was used to optimize cost or pro t. It was subsequently applied to engineering problems and nally to biological models [21]. The goal of di erent therapies is to destroy the tumor cells, while maintaining adequate amounts of healthy tissue. From a mathematical point of view, adequate destruction of tumor cells might mean forcing the system out of the basin of an unhealthy spiral node, or out of a limit cycle, and into the basin of attraction of a stable, tumor-free equilibrium. Alternatively, if the therapy pushes the system into a limit cycle in which the size of the tumor is small for a long period of time (as long as the life of the patient, for example), this could also be considered as ?cure? [20]. Optimality in treatment might be de ned in a variety of ways. The general goal is to keep the patient healthy while killing the tumor. In this thesis, we choose to minimize the tumor population, while constraining the normal cells to stay above some minimal level. Therefore, the development of therapies protocol can be phrased as an optimal control problem with constraints: for a xed time interval, nd the points within that interval at which the drug should be administered so that the number of tumor cells has been minimized, while the number of healthy cells has been kept above a prescribed threshold [20]. There have been many models, which tried to focus on simulating single or multiple important elements of the multifaceted process of tumor growth and response to therapy, with or without considering the optimal control. Based on these studies, people are trying 4 to design or improve treatment protocol by employing the tools of optimal control theory. See [13], [14], [15], [1] [16] [17]. This thesis is based on the models developed in [18] [19] [20] [21]. We build four di erent systems by stating the biological behaviors of each population, specifying most of the func- tional terms of them, and summarizing the mathematical form of these systems. We analyze the qualitative properties of them, especially for the fourth one, to which we consider the existence and characterization of the optimal bang-bang control. The main mathematical tools that we employ here are from Ordinary Di erential Equa- tions (ODE) and Optimal Control Theory. The organization of this thesis is as follows: We rst give a wide range of background, including some fundamental de nitions and impor- tant theorems from ODE. Afterward, we formulate the four systems of ordinary di erential equations. Next, we investigate the qualitative details of these systems one by one, where we frequently use the theorems we stated before. In the end, we give two appendices, one of which gives the de nitions and theorems concerning the One-point Compacti cation, the other of which contains a brief summary of the optimal control theory, including the problem, formulation, hypothesis, and solution outline. 5 Chapter 2 Background and Preliminary Knowledge De nition 2.1. (Hyperspace topology [22]) Suppose that X is a topological space. Then the hyperspace of X, denoted by 2X, is the space of compact subsets of X. Suppose that U1; ;Un is a nite collection of open subset of X, then R(U1; ;Un) =fK22XjK [ni=1Ui and for all 1 i n;k\Ui6=;g is a basic open set for the hyperspace topology. Theorem 2.2. [22] If X is compact, so is 2X. Theorem 2.3. [22] If X is a metric space with the bounded metric d, then the associated Hausdor metric generated the topology of 2X. De nition 2.4. (Asymptotically autonomous process with limit semi ow [23]) Assume (X;d) is a metric space, t0 2 R, and = f(t;s)jt0 s t < 1g. A continuous mapping : X!X is called a nonautonomous process if it satis es: (i) (s;s;x) = x;s t0. (ii) (t;s; (s;r;x)) = (t;r;x);t s r t0. The process is called autonomous if, additionally, (iii) (t+r;s+r;x) = (t;s;r);t s t0;r> 0. In this case, we set (t;x) = (t+t0;t0;x), and call an autonomous semi ow. 6 Moreover, the nonautonomous process is called asymptotically autonomous with limit semi- ow , if is an autonomous semi ow on X, and (ti +si;si;xi)! (t;x) as i!1 in the metric space X for any sequences ti!t;si!1;xi!x. De nition 2.5. (Chain Recurrent [23]) Let X and be as in de nition 2.4, A X be a nonempty positively invariant set under , and x;y2A. For " > 0;t > 0, an (";t)-chain from x to y (in A) is a sequence fx = x1;x2;:::;xn+1 = y;t1;t2;:::;tng of points xi in A and time ti > t such that dist( (ti;xi);xi+1) < ";i = 1;2;:::;n. A point x2A is called chain recurrent in A if for every "> 0;t> 0; there exists an (";t)-chain from x to x in A. The set A is said to be chain recurrent if every point in A is chain recurrent in A. Remark 2.6. We are interested in the case that A is compact, connected, and invariant. De nition 2.7. (!-limit set [23] [24]) If Phi is a nonautonomous semi ow on X and (s;x)2[t0;1) X, then the forward orbit of through (s;x) is de ned to be O (s;x) =f (t;s;x) : t sg X: IfO (s;x) has compact closure in X, then the !-limit set of (s;x) (or ofO (s;x)) is de ned by ! (s;x) = \ s f (t;s;x) : t g; where, for a subset A of X, A denotes the closure of A in X. In other words, y2! (s;x) if there is a sequence tj!1;tj >s, such that (tj;s;x)!y;j!1. In the case of an autonomous semi ow , the !-limit set is independent of s and hence we denote it by ! (x): ! (x) = \ 0 f (t;x) : t g: 7 An equivalent de nition for !-limit set for an autonomous semi ow is as follows: sup- pose that t is a autonomous semi ow on Rn and p 2 Rn. A point x in Rn is called an !-limit point of the orbit through p if there is a sequence of numbers t1 t2 t3 such that ti !1 and limi!1 ti(p) = x. The collection of all such !-limit points is called the !-limit set of p. Lemma 2.8. [25] Suppose ffng;n = 1;2;:::; is a sequence of functions de ned and con- tinuous on an open set D Rn+1 with limn!1fn = f0 uniformly on compact subsets of D. Suppose (tn;xn) is a sequence of points in D convergent to (t0;x0) 2 D as n !1, and let n(t);n = 0;1;:::; be a solution of the equation _x = fn(t;x) passing through the point (tn;xn). If 0(t) is de ned on [a;b] and is unique, then there is an integer n0 such that each n(t);n n0 can be de ned on [a;b] and converges to 0(t) uniformly on [a;b]. Proof. Please refer to [25] for the detailed proof. De nition 2.9. Again, (X;d) is a metric space, and t0 2R. Consider the following systems _x = f(t;x) (2.1) _y = g(y) (2.2) where x;y 2 Rn;f : R+ Rn ! Rn and g : Rn ! Rn are continuous vector functions. Assume the initial value problems for each system have unique solutions de ned for all future time after the initial time. De ne = (t;s;x0) to be the solution of system (2.1) with the initial value being x(s) = x0, and = (t;x0) to be the solution of system (2.2) with the initial value being y(0) = x0. Proposition 2.10. We use the same assumptions and notations as in the de nition 2.9. Then is asymptotically autonomous with limit semi ow if one of the two following conditions is satis ed: 8 (A) : f(t;x)!g(x) as t!1 uniformly on compact subsets of Rn. (B) : g is locally Lipschitz, and, for each compact subset K Rn, there is a function K : [0;1)![0;1) satisfying K(t)!0 as t!1 and Z t+ t [f(s;x) g(x)]ds K(t) for every (x; )2K [0;1] and t> 0. Proof. By de nition, we want to show (tj + sj;sj;xj) ! (t;x) as j!1 for any three sequences tj!t;sj!1;xj!x as j!1. (A) : Let T > 0 be such that tj 0 such that K is compact and [ t2[0;T] j 0 fy :jy (t;xj)j<"g K: Since g is locally Lipschitz on K, let L be the associated Lipschitz constant, and choose > 0 such that eLT <". Then by our assumption, there exists a function K(T) on [0;1) satisfying K(t)!0 as t!1 and Z t+ t [f(s;x) g(x)]ds K(t) 9 for x2K and 2 [0;T]. Since K(t) ! 0, choose S > 0 such that K(s) < for s > S. Denote uj = (t+sj;sj;xj);vj(t) = (t;xj). Then juj(t) vj(t)j = Z t 0 f( +sj;uj( )) d +xj Z t 0 g(vj( )) d +xj Z sj+t sj [f( ;uj( sj)) g(vj( sj))] d + Z t 0 jg(uj( )) g(vj( ))jd K(sj) +L Z t 0 juj( ) vj( )jd : By the Gronwall?s inequality, juj(t) vj(t)j k(sj)eLt eLT " for 0 t T. By the de nition of K, it follows that uj(t)2K for all T 2[0;T] if sj >S and uj!u;j!1 uniformly. As in part (A), we nish the proof. Lemma 2.11. We use the same assumptions and notations as in the de nition 2.5. Assume A is connected and chain recurrent with representing the associated autonomous semi ow. Then for any x;y2A;"> 0;T > 0, there exists a (";T)-chain from x to y. Proof. Let = "3. Since A is connected, there exists a1; ;an 2A with x = a1;y = an, and d(ai;ai+1) < ; ;i = 1; ;n 1: Now, since A is chain recurrent, for each 1 i n, there exist a ( ;T)-chain from ai to ai: fai = bi1; ;bin = ai;ti1; ;tin 1g;bj2A;tj T: 10 To get a chain from x to y, we can ?connect? these chains together, namely, we get a new sequence fx = a1 = b11; ;b1n;b21; ;bnn = an = y;t11; ;tnn 1g;bij 2A;tij T: To verify that this is a (";T)-chain, it su ces to check whether or not those points where the chain is connected in the above way satisfy the de nition of (";T)-chain. For 1 i n 1, d( (tin;bin);bi+11) = d( (tin;ai);ai+1) d( (tin;ai);ai) +d(ai;ai+1) < + < ": Lemma 2.12. We use the same assumptions and notations as in the de nition 2.4, where is an autonomous semi ow. Let T > 0 and OT (x) =f (t;x) : t Tg. Given y2! (x) and "> 0;t0 > 0. There exists a (";t0)-chain fy = y1; ;yl;yl+1 = y;t1; ;tlg such that yi2OT (x) for i = 1;2; ;l, ti = t0 for i = 1;2;:::;l 1 and t0 tl < 2t0. Proof. Since y2! (x), there exists n !1 such that ( n;x) !y as n!1. Choose N > 0 such that d( ( N +t;x); (t;y)) <" for t2[0;t0]: Similarly, we can choose M >N > 0 such that M > N + 2t0 and d( ( M;x);y) < : Let l be such that M N = lt0 +r for some r2[0;t0). 11 Now let y1 = y;yi = ( N + (i 1)t0;x) for i = 2;:::;l, yl+1 = y, ti = t0 for i = 1;:::;l 1, tl = t0 + r. We can check that d( (t1;y1);y2) = d( ( N + t0;x); (t0;y)) <", d( (ti;yi);yi+1) = 0 for i = 2;:::;l 1, and d( (tl;yl);yl+1) = d( ( N +lt0 +r;x);y) = d( ( M;x);t) <": Proposition 2.13. We use the same assumptions and notations as in the de nition 2.5, where is an autonomous semi ow, and suppose thatO (x) =f (t;x) : t> 0ghas compact closure in X. Then ! (x) has the following properties: (a) ! (x) is nonempty, compact, and connected. (b) ! (x) is invariant. (c) ! (x) attracts (t;x): dX( (t;x);! (x))!0;t!1: (d) ! (x) is is chain recurrent. Proof. (a), (b), and (c) are classical results. Check [24] To prove (d), we need to show that for any y2! (x), " > 0, and t0 > 0, there exists an (";t0)-chain lying in ! (x) which connects y to itself. From lemma 2.12, for each n2N, there exists a ( 1n;t0)-chain: fy = yn1;:::;ynln+1 = y;tn1;:::;tnlng with yi2OT (x) for i = 1;:::;ln + 1. Furthermore, We denote the sets fynigln+1i=1 by Yn. First we notice that Yn are all nite therefore compact and Yn OT (x) which is also compact. Apple Theorem 2.2 and consider OT (x) as the whole space, then, W.L.O.G., 12 Yn!eY in the hyperspace topology with some eY OT (x) being compact. (Otherwise we take a subsequence of Yn.) Since y 2 Yn for all n, then y 2 eY ! (x). Next, we will construct our chain by points from eY. Since eY is compact, there exists a > 0 such that if ey1;ey2 2eY and d(ey1;ey2) < , then d( (t;ey1); (t;ey1)) < "3 for t2[0;2t0]: Furthermore, let N > 0 be such that 1N < "3 and D(YN;eY) < "3, where D is the Hausdor metric. We nish the proof by carrying out the following steps. 1. Let ey1 = y, pick up ey2 2eY such that d(ey2;yN2 ) < , and t1 = t0. Then d( (t0;ey1);ey2) = d( (t0;yN1 );ey2) < d( (t0;yN1 );yN2 ) +d(yN2 ;ey2) < 1N + <": Also, since d(ey2;yN2 ) < , by how we chose , d( (t0;ey2); (t0;yN2 )) < "3. 2. Pick up ey3 2eY such that d(ey3;yN3 ) < , and t2 = t0. Then d( (t0;ey2);ey3) < d( (t0;ey2); (t0;yN2 )) +d( (t0;yN2 );yN3 ) +d(yN3 ;ey3) < "3 + 1N + <": Also, d(ey3;yN3 ) < induces d( (t0;ey3); (t0;yN3 )) < "3. ... (lN): Set eylN+1 = y and tlN 2[t0;2t0), we have d( (tlN;fylN);eylN+1) < d( (tlN;fylN); (tlN;yNlN)) +d( (tlN;yNlN);yNlN+1) < "3 + 1N <": 13 Therefore, fy = ey1;:::;eylN+1;t1;:::;tlNg is the desired chain. Remark 2.14. If is a nonautonomous process but asymptotically autonomous with limit semi ow , then and can be embedded in a single autonomous semi ow on a larger metric space Z = [t0;1] X, where [t0;1] is the one-point compacti cation of [t0;1) in the usual sense (see Appendix A.), and a metric on Z is de ned by ((s;x);(t;y)) = jh(s) h(t)j+d(x;y), where h : [t0;1]![0;1] is the map de ned as follows: h(t) = 8 >>< >>: t t0 1 +t t0; t<1; 1; t =1: The embedding is: : [0;1) Z!Z, (t;(s;x)) = 8 >>< >>: (t+s; (t+s;s;x)); t0 s<1 (1; (t;x)); s =1 (2.3) Clearly, is continuous and a semi ow on Z. Lemma 2.15. We use the same assumptions and notations as in the de nition 2.4, where is an asymptotically autonomous process with limit semi ow , and assume O (s;x) has compact closure in X. Then O (s;x) has compact closure in Z and f1g ! (s;x) = ! (s;x): (2.4) Proof. Since is the limit semi ow, by the de nition, sj!1, therefore we take the second embedding of (2.3). (2.4) follows. The topology on Z is the product topology on f1g X. Since O (s;x) has compact closure in X, by the de nition of compactness, it has a nite open cover, say C, in X. Then 14 we can construct an open cover C0 of O (s;x) by setting C0 =f[t0;1] EjE2Cg: It follows immediately that O (s;x) has compact closure in Z. Theorem 2.16. We use the same assumptions and notations as in the de nition 2.4, where is an asymptotically autonomous process with limit semi ow , and assume the forward orbit O+ (s;x) has compact closure in X. The !-limit set ! (s;x) has the following proper- ties: (a) ! (s;x) is nonempty, compact, and connected. (b) ! (s;x) is invariant under the semi ow : (t;! (s;x)) = ! (s;x) for each t 0: (c) ! (s;x) attracts (t;s;x): dX( (t;s;x);!)!0;t!1: (d) ! (s;x) is is chain recurrent for . Proof. By lemma 2.15,O (s;x) has compact closure in Z, then by proposition 2.13, ! (s;x) is nonempty, compact, and connected. (a) then follows. Similarly, ! (s;x) is invariant under , therefore we proved (b). For (c), we observe that dZ( (t;(s;x));! (s;x)) = t+s t0 1 +t+s t0 1 +dX( (t;s;x);! (s;x))!0 as t!1: This induces dX( (t;s;x);! (s;x)) ! 0 as t!1. Finally, f1g ! (s;x) being chain recurrent in Z implied that ! (s;x) is also chain recurrent in X. 15 Theorem 2.17. (Poincar e-Bendixson Theorem in R2 [26]) Suppose f2C1(E) where E is an open subset of R2, and is the solution ow of the system _x = f(x). If is a nonempty compact !-limit set of , and does not contain a rest point, then is a periodic orbit. Corollary 2.18. [26] We use the same notations and assumptions as in the Theorem 2.17. If E contains a periodic orbit of the system _x = f(x) and its interior U, then U contains at least one rest point of the system. Theorem 2.19. (Dulac?s Criterion [24]) Consider a smooth di erential equation system _x = f(x;y); _y = g(x;y): If there is a smooth function B(x;y) de ned on a simply connected region Rn such that @ @x(B f) + @ @y(B g) is not identically zero and of a xed sign on , then the system has no periodic solution on . Theorem 2.20. (Lyapunov?s Stability Theorem [24]) Consider autonomous system _x = f(x);x2Rn: Let x0 be a rest point of this system and U Rn be an open set containing x0. A continuous function h : U!R is called a Lyapunov Function of the above system at x0 if it satis es: 1. h(x0) = 0. 2. h(x) > 0 for x2Unfx0g. 3. h is continuously di erentiable on the set Unfx0g, and, on this set, _h(x) = gradh(x) f(x) 0: h(x) is called a strict Lyapunov Function if, additionally, 16 (4) _h(x) < 0 on Unfx0g. For the above system, if there exists a Lyapunov Function de ned on an open neighborhood of a rest point of the system, then this rest point is stable. Moreover, if the Lyapunov Function is a strict Lyapunov Function, then this rest point is asymptotically stable. De nition 2.21. (Subsolution and Supersolution [27]) Let f(t;x) be de ned in D R2. Consider the following system with initial value _y = f(t;x) for t2J = [t0;t0 +a]; y(t0) = K; (2.5) where K is a constant, and a> 0. Suppose y(t) is a solution of the system. y /y+ is called a subsolution/supersolution of the above system if it is di erentiable on J and satis es _y f(t;x) on J; y (t0) K: subsolution _y+ f(t;x) on J; y+(t0) K: supersolution Theorem 2.22. (Kamke?s Comparison Theorem, version of unique solution [27]) Under the previous de nition, where f is a continuous function, and D is an open subset of R2, if the initial value problem (2.5) has a unique solution on J, then y (t) y(t) y+(t) on J: Theorem 2.23. (Carath eodory?s Existence Theorem [28]) Consider the initial value problem _x = f(t;x(t)); x( ) = ; where ( ; ) 2 D, with D a nonempty open subset of R Rn, and f : D ! Rn. It has a solution if for some open set Ra;b D centered at ; , the restriction of f to Ra;b is 17 continuous in x for xed t, measurable in t for xed x, and satis es jf(t;x)j m(t); (t;x)2Ra;b; for some m integrable over the interval [ a; +a]. Next we will introduce several de nitions and theorems from optimal control theory. Please refer to Appendix B for the fundamental background, such as basic terminologies and general problem formulation, and an outline of solving optimal (bang-bang) control problem. De nition 2.24. (Optimal control [29]) Consider the following optimal control problem: min J = Z T 0 F(x(t);u(t);t) dt+S(x(T);T) (2.6) _x(t) = f(x(t);u(t);t); x(0) = x0 (2.7) g(x(t);u(t);t) 0 (2.8) h(x(t);t) 0 (2.9) a(x(T);T) 0 (2.10) b(x(T);T) = 0 (2.11) where T is free on [0;tf], F : Rn Rm R!R;S : Rn R!R;f : Rn Rm R!Rn;g : Rn Rm R!Rs;h : Rn R!Rq;a : Rn R!Rl, and b : Rn R!Rl0. Assume that F, S, f, g, h, a, and b are continuous over their domain, respectively. De ne the (state-dependent) set of admissible values (x;t) =fu2Rmjg(x;u;t) 0g Rm 18 and the set N(x;t) =f(F(x;u;t) + ;f(x;u;t))j 0;u2 (x;t)g Rn+1 where m and n are the number of control and state variables, respectively. We are looking for a measurable function u( ) mapping from [0;T] into Rm and a cor- responding function x( ) mapping from [0;T] into Rn which is absolutely continuous, such that the constraint (2.7) | (2.11) are satis ed and the objective functional (2.6) takes its minimum value. Theorem 2.25. (Filippov-Cesari Theorem [29]) Under the previous de nition, we assume that F;S;f;g;h;a; and b are continuous in all their arguments at all points (x;u;t). Fur- thermore, suppose that the following conditions hold: 1. There exists an admissible solution pair. 2. N(x;t) is convex for all (x;t)2Rn [0;tf]. 3. There exists > 0 such that kx(t)k< for all admissible fx(t);u(t)g and t. 4. There exists 1 > 0 such that kuk< 1 for all u2 (x;t) with kxk< . Then there exists an optimal triple fT ;x ;u g with u ( ) measurable. Theorem 2.26. (Pontryagin?s Maximum/Minimum Principle, or PMP [21], [30]) Use the same notations and assumptions as in the de nition 2.24, where we assume u(t) = (u1(t);:::;um(t)) is a piecewise continuous control function and x(t) = (x1(t);:::;xn(t)) be the corresponding continuous and piecewise di erentiable state function de ned on the xed interval [t0;t1] that minimizes Z t1 t0 f(t;x(t);u(t)) dt 19 subject to the di erential equations _xi(t) = gi(t;x(t);u(t)); i = 1;:::;n; initial conditions xi(t0) = xi0; xi0 xed; i = 1;:::;n; terminal conditions xi(t1) = xi1; xi0 xed; i = 1;:::;p; xi(t1) xit; xi0 xed; i = p+ 1;:::;q xi(t1) is free; i = q + 1;:::;n; and control variable restriction u(t)2U; U is a given set in Rm: We assume that f;g;@f=@xj, and @gi=@xj are continuous functions of all their arguments for all i = 1;:::;n and j = 1;:::;n. Then there exist continuous functions (adjoint functions) (t) = ( 1(t);:::; n(t)) : R!Rn, where for all t0 t t1 we have (t)6= 0 such that for every t0 t t1,, H(t;x (t);u(t); (t)) H(t;x (t);u (t); (t)); where the Hamiltonian function H is de ned by H(t;x;u; ) = f(t;x;u) + nX i=1 igi(t;x;u): Moreover, except at points of discontinuity of u (t), _ (t) = @H(t;x (t);u (t); (t))=@xi; i = 1;:::;n: 20 Finally, the following transversality conditions are satis ed: i(t1) no conditions; i = 1::::;p; (2.12) i(t1) 0; the equality holds if x i(t1) >xi1; ;i = p+ 1;:::;q; (2.13) i(t1) = 0; ;i = q + 1;:::;n: (2.14) 21 Chapter 3 Model Description 3.1 General assumptions We will study three types of therapies in this thesis: Immunotherapy, Chemotherapy, and Radiotherapy. What we list below are several assumptions about these therapies that some models in this thesis share, each of which represents a certain well-accepted behavior of a population or between populations. Immune response: The intrinsic growth rates of tumor and immune cell populations are assumed to obey logistic laws. Also, the immune cells will kill the tumor cells at some rate, which results in the death of both populations. Furthermore, in the absence of any tumor cell, the immune cells will die o at a constant rate. Chemotherapy: Medicines will kill not only tumor cells, but also host and immune cells, of which the rate is modeled by the term k(1 e M) with k denoting the response coe cient. The parameters k and can be adjusted according to clinical data. Radiation: The radiation a ects all types of cell. Furthermore, it is sometimes possible for broken chromosomes to recombine, so that ?broken? tumor or host cells can become viable cells again. This is here modeled to occur at a constant rate. 3.2 Model I: Single population with Radiotherapy only We start with the case with radiotherapy only. More speci cally, we deal with one cell population (either host or tumor cells) and establish some general mathematical properties and dynamics of the system. We will study the double-population case in the next model. 22 A(t) and Ar(t) denote the population densities at time t of the cell populations under consideration and the associated radiated cell populations, respectively. We assume logistic growth for the cell population with r and K denoting the intrinsic growth rate and carrying capacity, respectively. Radiated cells with broken chromosomes are represented by D(t)u, where D(t) is the rate of radiation protocol, a nonnegative function. It?s reasonable to ask furthermore that D(t)6 0. We consider two cases for D(t) in this model: 1. D(t) = D0 > 0, constant. 2. D(t) = D0e t, decay. p is the rate at which the radiated cells recombine into normal cells, and is the washout rate of radiated cells. After the discussion above, the system turns out to be of the following form: 8 >>>> >< >>> >>: _A = rA(1 A K) D(t)A+pAr _Ar = D(t)A pAr Ar A(0) 0;Ar(0) 0 (3.1) 3.3 Model II: Double populations with Radiotherapy only In this model, we inherit the same assumptions as in Model I and consider the case with both host and tumor cells, which includes additional competition between these two populations. H(t) and T(t) denote the population densities at time t of the host and tumor cell populations, respectively, with Hr(t) and Tr(t) for radiated host and tumor cell popu- lations, respectively; ri andKi denote the intrinsic growth rates and carrying capacities, respectively. 23 The radiotherapy is designed in such a way that full radiation concentration a ects the cancer cells, and a small proportion of the radiation a ects the host cells. We only consider the case with constant radiation rate here, i.e., D(t) = D, where D is a constant. We assume the same washout rate, , but di erent recombining rates for these two populations, which are denoted by p1 and p2 for the radiated host and tumor cells, respectively. The system, in this way, becomes 8 >>>> >>> >>>> >>< >>>> >>> >>> >>>: _H = r1H(1 H K1 ) "DH +p1Hr c1HT _Hr = "DH p1Hr Hr _T = r2T(1 T K2 ) DT +p2Tr c2HT _Tr = DT p2Tr Tr H(0) 0;Hr(0) 0;T(0) 0;Tr(0) 0 (3.2) 3.4 Model III: Tumor model with Immune Resistance and Chemotherapy Next, we will look at some more detailed cells population behaviors. In this model, we will introduce three populations: tumor cells, immune cells, and normal (host) cells, of which the population densities at time t are denoted by T(t);I(t), and H(t), respectively. We assume that the immune cells have a source with a constant in ux rate s. Moreover, in the absence of any tumor, they will die o at a per capita rate d, so that the immune cells population won?t blow up, but has an upper bound sd. The presence of tumor cells stimulate the immune response, represented by the Michaelis- Menten form I(t)T(t) +T(t) , which is wildly used in enzyme kinetics modeling (the same as 24 the terms used in [31],[32]). It?s easy to see that this term is positive and monotonically increasing with respect to T. We will use similar terms again in the next model. The reaction between immune and tumor cells would result in the death for both, which leads to two competition terms: c1I(t)T(t) and c2I(t)T(t): Both the tumor and normal cell populations obey logistic growth laws, with ri and Ki representing the intrinsic growth rates and carrying capacities for these two types of cells, respectively. Thus, we have the following growth terms: r2T(1 TK 2 ) and r3H(1 HK 3 ): In addition, there are two terms representing the competition between tumor and host cells: c02T(t)H(t) and c3T(t)H(t): We add the e ect of drugs into the system. M(t) denote the amount of drugs at time t. By section 3.1, the response rates for all three types of cells are given by the terms: ki(1 e iM); i = 1;2;3 where for mathematical convenience and unknown of details of pharmacokinetics, we let i = 1;i = 1;2;3 in the preliminary studies. The amount of drug is determined by the given dose, v(t), and a per capita decay rate of the drug once it?s injected, at the rate d4. 25 Thus, we model this system as 8> >>> >>> >>>< >>> >>>> >>> : _I = s d1I + IT +T c1IT k1(1 e M)I _T = r2T(1 T K2 ) c2IT c 0 2TH k2(1 e M)T _H = r3H(1 H K3 ) c3TH k3(1 e M)H _M = v(t) d4M (3.3) 3.5 Model IV: Immuno-Chemotherapy with controls In this model, we introduce a mixed-therapy method to treat cancer. More speci cally, we combine two biological therapies, Tumor Speci c T Cell (CD8+T CTL) and cytokine IL-2, together with a type of chemotherapy. We also consider controls of these therapies in the system. There are six populations in this model. The population densities at time t are denoted by: T(t): Tumor cells N(t): Natural Killer cells C(t): Number of circulating lymphocytes S(t): Tumor Speci c T cells (CD8+T cells) I(t): Immunotherapy (cytokine IL-2) concentration M(t): Chemotherapy medicine concentration We assume the similar functional terms to the ones in [14]. The functions of ac- tion/reaction between these populations can be categorized into several types. Growth and death terms: The growth of tumor cells are assumed to obey logistic laws. Both Natural Killer cells and CD8+T cells have constant death rates, while the growth 26 of the former comes partially from circulating lymphocytes. The circulating lympho- cytes have a constant production rate and a constant death rate. Both chemotherapy and immunotherapy are assumed to decay proportionally to their concentrations. Killing terms: The tumor cells will be killed by both the Natural Killer cells and the CD8+T cells. Also, the chemotherapy medicine will a ect all types of cells. The form of inhibition is the same as in the system (3.1). Stimulation and recruitment: The IL-2 stimulates the recruitment of both the Natural Killer cells and the CD8+T cells. The latter will also be stimulated by the debris of tumor cells. On the other hand, the CD8+T cells will also stimulate the produc- tion of IL-2. These stimulations shall be modeled by the Michaelis-Menten form [33]. The Natural Killer cells killing the tumor cells may enhance the production of the CD8+T cells. The same e ect will be produced by the encounters between circulating lymphocytes and tumor cells. Control terms: We can control the population concentration of CTLs, IL-2, and medicine. Before building the whole model, let?s discuss the detailed mathematical terms rst in the following text. T(t) : The same as before, simple logistic growth term of the tumor cell population, rT(1 T K), is assumed in the absence of medicine and immune interactions. Death of tumor cells due to the Natural Killer cells and CD8+T cells is modeled by a mass action term, cNT, and a ratio dependent term, D(S;T)T = d (S=T) l s=nl + (S=T)lT [33], respectively. Death due to chemotherapy medicine is given in the form that we used to assume, kT(1 e TM)T. N(t) : The constant rate source term of Natural Killer cells due to circulating lymphocytes is eC, while the linear natural death term is fN. By killing the tumor cells, there 27 is an extra death term, pNT. Stimulation due to IL-2 is modeled in the Michaelis- Menten form, pNNIg N +I . The same as before, death due to chemotherapy medicine is given by kN(1 e NM)N. C(t) : The circulating lymphocytes has a constant source term, , and a linear death term, C, together with the death due to medicine, kC(1 e CM)C. S(t) : For the CD8+T cells, we assume there exists a modi ed (in the presence of IL-2) linear death term, mS +I, as well as a modi ed quadratic death term, u0S 2CI +I , due to the inhibition by circulating lymphocytes. CTLs will also die through interaction with tumor cells, modeled by qST. Also, the stimulations on the CD8+T cells by Natural Killer cells and circulating lymphocytes interacting with tumor cells are represented by a1NT and a2CT, respectively. The stimulatory e ect on the CD8+T cells by the debris of tumor cells and the IL-2 are represented by pTSTg T +T and pISIg I +I , respectively. Similarly, death due to medicine is modeled by kS(1 e SM)S. I(t) : The IL-2 is assumed to be with a linear death rate, II, and a constant source from Circulating Lymphocytes, C. We use pSSIg S +I to represent the stimulatory production from the CD8+T cells. M(t) : Once injected, medicine (chemotherapy) is assumed to have a linear decay rate, M. All of the controls are in the term of XvX(t), where X denotes corresponding populations. After the discussion above, we can summarize the system as 28 8 >>>> >>> >>> >>>> >>> >>>> >< >>> >>> >>>> >>> >>>> >>> >>: _T = rT(1 T K) cNT D(S;T)T kT(1 e TM)T _N = eC fN pNT + pNNI gN +I kN(1 e NM)N _C = C kC(1 e CM)C _S = mS +I u0S2CI +I qST +a1NT +a2CT + pISIg I +I + pTSTg T +T kS(1 e SM)S + SvS(t) _I = II + C + pSSI gS +I + IvI(t) _M = M + MvM(t) (3.4) where D(S;T) = d (S=T) l s=nl + (S=T)l. Also, T(0) 0;N(0) 0;C(0) 0;S(0) 0;I(0) 0;M(0) 0. 29 Chapter 4 Analytic Results 4.1 Model I Before starting to analyze the three cases individually, we rst take a look at some preliminary properties of the system (3.1). Proposition 4.1. Assume that (A(t);Ar(t)) is the solution of the initial value problem (3.1), and D(t) is continuous, then A(t) 0, Ar(t) 0. Proof. (0;0) is an equilibrium. It su ces to prove that A(t) and Ar(t) cannot cross the A-axis and Ar-axis. By the continuity of D(t), A(t) and Ar(t) are at least continuous over a small time interval [0; ] since they are the solutions of the initial value problem (3.1). Assume A(0) = 0;Ar(0) > 0, then _A(0) = pAr(0) > 0, hence there exists a 0 < " < such that _A(t) > 0;Ar(t) > 0 for t2[0;"]. Thus, A(t) = Z t 0 _A(t) dt> 0 for t2(0;"]; Similar argument can imply the same conclusion for A(0) > 0;Ar(0) = 0, as _Ar(0) = D(0)Ar(0) > 0. For the case that A(0) > 0;Ar(0) > 0, if the solution orbit intersects either A-axis or Ar-axis at some time point, for instance, t1 >", we use the same arguments above to show that the orbit will stay in the rst quadrant, with A(0) and Ar(0) being substituted by A(t1) and Ar(t1), respectively. 30 Proposition 4.2. IVP System (3.1) is dissipative within the set R= R2+\ (A;Ar)jA+Ar max A(0) +Ar(0);K(r + ) 2 4r : where R2+ =f(A;Ar)jA 0;Ar 0g. Proof. Let X = A + Ar, then we combine the two equations in the initial value problem (3.1) together: _X = rA(1 A K) Ar = X rKA2 + (r + )A = X + K(r + ) 2 4r r K(A K(r + ) 2r ) 2 X + K(r + ) 2 4r Solving the corresponding equality, we get eX(t) = X(0) K(r + ) 2 4r e t + K(r + ) 2 4r : Therefore, by the Kamke comparison theorem (Theorem 2.22) , X max X(0);K(r + ) 2 4r where X(0) = A(0) +Ar(0). Next, we will discuss cases with di erent types of D(t) and see some analytical properties for each of them. 4.1.1 Constant radiation In this case, D(t) = D0 > 0, where D0 is a constant. The system (3.1) becomes 31 _A = rA(1 A K) D0A+pAr _Ar = D0A pAr Ar A(0) 0;Ar(0) 0 (4.1) Then, by solving the corresponding equations 8> < >: rA(1 AK) D0A+pAr = 0 D0A pAr Ar = 0 we can get at most two reasonable distinct equilibria of the initial value problem (4.1): E0 = (0;0) and E1 = Kr(p+ ) D0r(p+ ) ;KD0r(p+ ) D0r(p+ )2 = (u1;v1)6= E0; u1 0;v1 0: Theorem 4.3. E1 exists if and only if E0 is unstable. Proof. For one direction, E1 existing means that u 0;v 0, and the equalities cannot hold at the same time. Since D0 6= 0, and all other parameters are positive in this model, it is also equivalent to say that r(p+ ) D0 > 0, which is 0 0; which indicates + has positive real part. Therefore, E0 is unstable. Similarly arguments yield the other direction. Theorem 4.4. If E1 does not exist, then E0 is asymptotically stable. Proof. If E1 does not exist, then E0 is the only equilibrium. By the dissipativity (proposition 4.2), it su ces to prove that there is no periodic solution. If there is a periodic orbit, then by the corollary 2.18, it must surround an equilibrium, which is E0 in this case. Use the dissipativity again, it?s clear that this type of periodic orbit does not exist, because if it does, then part of the orbit must lie in the dissipative region in proposition 4.2, therefore cannot keep being periodic but approaches E0. Hence, E0 is globally asymptotically stable. Theorem 4.5. If E1 exists, then it is globally asymptotically stable on R2nE0. Proof. From proposition 4.3, E0 is unstable. Since the system is dissipative, we only need to show that there is no periodic solution. We use the Dulac?s Criterion (Theorem 2.19) to prove this. Choose B(A;Ar) = 1A A r in the Dulac?s Criterion and apply it to the initial value problem (4.1): @ @A 1 A Ar rA 1 AK D0A+pAr + @@A r 1 A Ar(D0A pAr Ar) = rK A r pA2 D0A2 r < 0 33 This holds for any (A;Ar) 2R2+. Then by the Dulac?s Criterion, the initial value problem does not have any periodic solution in R2+. Theorem 4.6. The IVP system (4.1) persists if and only if 0 < D0 < r(p+ ) , otherwise both populations die o . Proof. Combine the above 3 theorems. The system persists if and only if E1 exists, which by Theorem 4.3 is equivalent with 0 < D0 < r(p+ ) . Furthermore, Theorem 4.5 guarantees it is asymptotically stable. If not, then by Theorem 4.4, any orbit will approach the equilibrium E1, the extinction. 4.1.2 Decay radiation In this case, we assume D(t) = D0e t. The system (3.1) becomes _A = rA(1 A K) D0e tA+pAr _Ar = D0e tA pAr Ar A(0) 0;Ar(0) 0 (4.4) Theorem 4.7. Let (A(t);Ar(t)) be the solution of the initial value problem (4.4), then limt!1(A(t);Ar(t)) = (0;0) or limt!1(A(t);Ar(t)) = (K;0): Proof. We consider the asymptotic system by taking the limit of t to 1. _A = rA(1 A K) +pAr _Ar = pAr Ar (4.5) 34 By Theorem 2.16, the !-limit sets of system (4.4) are contained within the chain recur- rent sets of system (4.5). But the only chain recurrent sets of (4.5) are (0;0) and (K;0). This follows from two observations: (1) For any orbit starting from a point of which the Ar-coordinate is not 0, by the second equation of the initial value problem (4.5), it will approach some point on the A-axis. Therefore, for any of these points, there would not be any chain connecting it to itself. (2) For any orbit starting from a point lying on the A-axis other than (0;0) and (K;0), the Ar-coordinate will remain the same; but by the rst equation of the initial value prob- lem (4.5), the evolution of A obeys the logistic law, which will approach K. Therefore, as the same reason as last case, there would not be any chain connecting any of these points to itself. 4.2 Model II Similarly as proposition 4.2, we get the following corollary: Corollary 4.8. The initial value problem (3.2) is dissipative within the set R= R4+ \ (H;Hr;T;Tr)jX max X(0);K1(r1 + ) 2 4r1 + K2(r2 + )2 4r2 : where R4+ =f(H;Hr;T;Tr)jH 0;Hr 0;T 0;Tr 0g and X = H +Hr +T +Tr. 35 Proof. The proof is similar as in the Theorem 4.2. Let X = H +Hr +T +Tr, then by (3.2), _X = r1H(1 H K1 ) Hr c1HT +r2T(1 T K2 ) Tr c2HT X + (r1H + H r1H 2 K1 ) c1HT + (r2T + T r2T2 K2 ) c2HT X r1K 1 [H K12r 1 (r1 + )]2 + K1(r1 + ) 2 4r1 c1HT r2K 2 [T K22r 2 (r2 + )]2 + K2(r2 + ) 2 4r2 c2HT X + K1(r1 + ) 2 4r1 + K2(r2 + )2 4r2 Solving the corresponding equality, we can get eX = X(0) K 1(r1 + )2 4r1 + K2(r2 + )2 4r2 e t + K1(r1 + ) 2 4r1 + K2(r2 + )2 4r2 : By using the Kamke comparison theorem (Theorem 2.22), we can get the desired result. Also, in absence of radiation, the system (3.2) reduces to a competition system 8 >>> >>> < >>> >>> : _H = r1H(1 H K1 ) c1HT _T = r2T(1 T K2 ) c2HT H(0) 0;T(0) 0: (4.6) A well-known result shows that cancer cells will win in the competition under the hypothesis c2K1 >>> >< >>> >>: r1H(1 HK 1 ) "DH +p1Hr c1HT = 0 "DH p1Hr Hr = 0 H(0) 0;T(0) 0: has a nonnegative solution. Since T = 0, by (4.2), we get the conditions 0 0;d > 0, in which the host cells die o . Similarly as above, the condition for the existence of this type of equilibrium is 0 >>> >>> >>>< >>> >>>> >>> : 2 64 r1 K1 c1 c2 r2K 2 3 75 2 64H T 3 75 = 2 64A B 3 75 Hr = "Dp 1 + H Tr = Dp 2 + T (4.12) where A = r1 "D + p1"Dp 1 + ;B = r2 D + p2Dp 2 + . Therefore, E3 exists if and only if the matrix 2 64 r1 K1 c1 c2 r2K 2 3 75 is invertible, namely, r1r2 6= c1c2K1K2: and 0 BB BB BB B@ H Hr T Tr 1 CC CC CC CA = 0 BB BB BB BB B@ r2K1A c1K1K2B r1r2 c1c2K1K2"D(r 2K1A c1K1K2B) (p1 + )(r1r2 c1c2K1K2) c2K1K2A r1K2B c1c2K1K2 r1r2D(c 2K1K2A r1K2B) (p2 + )(c1c2K1K2 r1r2) 1 CC CC CC CC CA 38 4.2.2 Stability From (4.3), we calculate the variational matrix being 2 66 66 66 66 4 r1 2r1HK 1 "D c1T p1 c1H 0 "D p1 0 0 c2T 0 r2 2r2TK 2 D c2H p2 0 0 D p2 3 77 77 77 77 5 = 2 64M1 M2 M3 M4 3 75 (4.13) where Mis are 2 2 matrices. If (^a;^b;^c; ^d) is a general equilibrium, by solving the characteristic equation, the eigenval- ues are given as follows. For mathematical convenience, denote = r1 2r1HK 1 "D c1T; = r2 2r2TK 2 D c2H. Since det(M2) = 0, det(M3) = 0, it?s easy to get 1 = 12 h ( p1 ) p ( +p1 + )2 + 4"Dp1 i 2 = 12 h ( p2 ) p ( +p2 + )2 + 4Dp2 i (4.14) Notice all of the eigenvalues are real. Null state. For E0, 01 = 12 h (r1 "D p1 ) p (r1 "D +p1 + )2 + 4"Dp1 i 02 = 12 h (r2 D p2 ) p (r2 D +p2 + )2 + 4Dp2 i By Theorem 4.3, E0 is unstable if and only if (4.8) and (4.10) are satis ed, i.e. E0 is stable if and only if D> max r 1(p1 + ) " ; r2(p2 + ) . 39 Tumor free. For E1, 11 = 12 (r1 2r1aK 1 "D p1 ) r (r1 2r1aK 1 "D +p1 + )2 + 4"Dp1 # (4.15) 12 = 12 h (r2 D c2a p2 ) p (r2 D c2a+p2 + )2 + 4Dp2 i (4.16) Theorem 4.11. E1 is stable if and only if D> max r 1(p1 + )(r2 c2K1 p2 ) r1(p1 + ) " c2K1 ; r1(p1 + )(p2 + )(r2 c2K1) [r1(p1 + ) "c2K1(p2 + )] Proof. Here, by (4.8), we assume su ciently small " to guarantee E1 exists. By The- orem 4.5, it su ces to make 12 < 0. By (4.16), 12 are real, and 12 < 0 ) 8 >< >: r2 D c2a p2 < 0 p(r 2 D c2a+p2 + )2 + 4Dp2 < >: [r1(p1 + ) " c2K1]D>r1(p1 + )(r2 c2K1 p2 ) [r1(p1 + ) "c2K1(p2 + )]D>r1(p1 + )(p2 + )(r2 c2K1) ) 8> >< >>: D> r1(p1 + )(r2 c2K1 p2 )r 1(p1 + ) " c2K1 D> r1(p1 + )(p2 + )(r2 c2K1) [r 1(p1 + ) "c2K1(p2 + )] The third step deduction is guaranteed by " being su ciently small. Theorem 4.12. (Global stability of the tumor free equilibrium) Recall from (4.9) that E1 = (a;b;0;0). Then E1 is globally asymptotically stable if 1. 2 2 2 < 0 40 2. (c1 +c2) 2 2c2 K2 < 0 3. r2 c2a< 0 Proof. We prove this by constructing a Lyapunov function. Let V(H;Hr;T;Tr) = H a aln Ha +Hr b bln Hrb +T +Tr (4.17) This function satis es the rst two conditions in Theorem 2.20. Now compute _V = (H a) r1(1 HK 1 ) "D +p1HrH c1T +(Hr b) "DHH r p1 +r2T(1 TK 2 ) c2HT Tr Since "D = r1(1 aK 1 ) + p1ba and ab = p1 + "D = (H a) r1(1 HK 1 ) r1(1 aK 1 ) p1ba + p1HrH +(Hr b) "DHH r "Dab +r2T(1 TK 2 ) c1T(H a) c2HT Tr = (H a) r1K 1 (H a) p1baH(H a) + p1H(Hr b) +(Hr b) "Da bHr (Hr b) + "D Hr (H a) c2K 2 T2 c1T(H a) c2T(H a) +r2T c2aT Tr = ( r1K 1 + p1baH)(H a)2 "DabH r (Hr b)2 + (p 1H + "DH r )(H a)(Hr b) c2K 2 T2 (c1 +c2)(H a)T + (r2 c2a)T Tr let = r1K 1 + p1baH; = "DabH r ; = p 1H + "DH r = (H a)2 (Hr b)2 + (H a)(Hr b) c2K 2 T2 (c1 +c2)(H a)T + (r2 c2a)T Tr 41 For mathematical convenience, we consider 2 _V = [(H a) (Hr b)]2 + ( 2 2 2 )(Hr b) 2 [(H a) c1 +c2 T]2 + (c 1 +c2)2 2c2 K2 T2 +2(r2 c2a)T 2 Tr Consequently, if 2 2 2 < 0; (c1 +c2)2 2c2 K2 < 0;and r2 c2a< 0 then _V < 0, which makes V satisfy the third condition of Theorem 2.20, therefore E1 is globally asymptotically stable. Dead. Corollary 4.13. The dead equilibrium E2 is stable if and only if D> max r 2(p2 + )(r1 c1K2 p1 ) "r2(p2 + ) c1K2 ; r2(p1 + )(p2 + )(r1 c1K2) ["r2(p2 + ) c1K2(p1 + )] 4.3 Model III To understand the dynamics of this system, we analyze the null-surfaces and equilibria of the drug-free system which is the simpli cation of system (3.3) without any term involving the medicine M(t). 42 8 >>> >>>< >>> >>>: _I = s d1I + IT +T c1IT _T = r2T(1 T K2 ) c2IT c 0 2TH _H = r3H(1 H K3 ) c3TH (4.18) 4.3.1 Null-surfaces The three sets of null-surfaces of the system (4.18) are as follows: N1: s d1I + IT +T c1IT = 0 )I = s( +T)(d 1 +c1T)( +T) + T = f(T): If (d1 + c1T)( + T) + T 6= 0. This is a curved surface parallel to the N-axis in the I-T-N surface. N2: r2T(1 TK 2 ) c2IT c02TH = 0 )T = 0 or T = K2 + c2K2r 2 I + c 0 2 r2H = g(I;H): N2 is a plane. N3: r3H(1 HK 3 ) c3TH = 0 )H = 0 or H = K3 + c3K3r 3 T = h(T): N3 is a plane which parallel to I-axis. 4.3.2 Equilibria There are three types of equilibria: Tumor free In this category, we consider the tumor population to be zero and the host cells population to be nonzero. The equilibrium is of the form ( sd 1 ;0;K3). 43 Dead We classify an equilibrium as \dead" if the host cells population is zero. There are two types of this category of equilibria. ( sd 1 ;0;0) in which both tumor and normal cells die o . (a;b;0) in which a;b satisfy f(b) = 0;g(a) = 0. The solution, if it exists, is unique upon xed parameters. In this case, only the normal cells died o and the tumor cells survived. Coexisting (a;b;c) in which a;b;c satisfying a = f(b);b = g(a;c);c = h(b); (4.19) which induces a = s( +b)(d 1 +c1b)( +b) + b b = K2 + c2K2r 2 a+ c 0 2 r2c c = c3K3r 3 b By solving the second and third linear equations and then plugging into the rst one, we can get a cubic equation of a, which means, depending on the parameters, there could be zero, one, two, or three di erent coexisting equilibria, which may not necessarily locate in the rst quadrant, but the whole R3, theoretically. 44 4.3.3 Stability We use linearization around the equilibria to analyze the stability. The linearized system is as follows 2 66 66 4 _u _v _w 3 77 77 5 = 2 66 66 64 d1 c1T I ( +T)2 c1I 0 c2T r2 2r2TK 2 c2I c02H c02T 0 c3H r3 2r3HK 3 c3T 3 77 77 75 2 66 66 4 u v w 3 77 77 5 (4.20) Tumor free equilibrium. In principle, we want the tumor free equilibrium to be stable so that the system will move toward the tumor free state starting at least locally. By linearization around this equilibrium, we can get the system 2 66 66 4 _u _v _w 3 77 77 5 = 2 66 66 64 d1 sd 1 c1sd 1 0 0 r2 c2sd 1 c02K3 0 0 c3K3 r3 3 77 77 75 2 66 66 4 u v w 3 77 77 5 (4.21) with eigenvalues 1 = d1 < 0 2 = r2 c2sd 1 c02K3 3 = r3 < 0 So the equilibrium is stable when 2 < 0, i.e. r2 < c2sd 1 +c02K3: (4.22) 45 As r2 is the growth rate of the tumor cells, which is normally large compared with c2 and c02, (4.22) tells us there is a big chance that the tumor free equilibrium is unstable. This indicates the importance or the therapies, for example chemotherapy. Dead equilibria. Similarly, we use (4.20) to analyze the stability of them. Of the rst type of dead equilibrium ( sd 1 ;0;0), the eigenvalues are 1 = d1 < 0 2 = r2 c2sd 1 3 = r3 > 0 which indicates this equilibrium is always unstable. Of the second type of dead equilibrium (a;b;0), 1 and 2 are the solutions of the equation 2 (d1 r2 + 2r2bK 2 +c2a) d1r2 + 2r2bd1K 2 +c2ad1 + c2b a ( +b)2 +c1c2ab = 0; and 3 = r3 c3b. So this equilibrium could be either stable or unstable, depending on the parameters. 4.4 Model IV For this model, we focus on investigating the qualitative properties of bang-bang controls in the system, therefore determining the optimal therapies treatments schedule for patients. To achieve this goal, we need to look at the existence and characterization of the optimal bang-bang control. For the background, please, refer to Appendix B. 46 4.4.1 Objective Functional In particular, we wish to minimize the objective functional J(vS;vI;vM) = Z tf t0 T(t) +"SvS +"IvI +"MvM dt (4.23) where "S;"I and "M are weight factors, t0 = 0, and tf is the terminal time which is xed. This control problem is under two constraints. The rst one requires the tumor cell population to be bounded throughout the time interval [0;tf]. It is formulated by T(tf) ; is a constant: (4.24) The other constraint con nes the total amount of Chemotherapy medicine used through the therapy. It is given in a integration form. Z tf t0 vm(t)dt ; is a constant: (4.25) 4.4.2 Existence of Optimal Control We establish the existence of an optimal control by the Filippov-Cesari existence theo- rem (Theorem 2.25). Theorem 4.14. (Existence of Optimal Control) Given the IVP system (3.4) under the constraints (4.24) and (4.25), with the set of all admissible controls being U =fu = (vs(t);vI(t);vM(t))jvS(t);vI(t);vM(t)2[0;1]piecewise continuousg; then there exists an optimal control bu = (bvS;bvI;cvM) for the objective functional (4.23), namely, J(bu) = min u2[0;1]3 J(u): 47 Proof. To apply the Filippov-Cesari existence theorem (Theorem 2.25), we need to verify its four conditions. Use the notation of Theorem 2.25 for our system, then x = 0 BB BB BB BB BB BB BB @ T N C S I M 1 CC CC CC CC CC CC CC A where x2R6, and N(x;t) : R6 R+ !R7: N(x;t) = 0 BB BB BB BB BB BB BB BB BB B@ T(t) +"SvS +"IvI +"MvM + rT(1 TK) cNT D(S;T)T kT(1 e TM)T eC fN pNT + pNNIg N +I kN(1 e NM)N C kC(1 e CM)C mS +I u0S 2CI +I qST +a1NT +a2CT + pISI gI +I + pTST gT +T kS(1 e SM)S + SvS(t) II + C + pSSIg S +I + IvI(t) M + MvM(t) 1 CC CC CC CC CC CC CC CC CC CA where 0. The rst one asks the existence of an admissible solution pair for the state and controls, which was showed in [34] and [35]. For the second condition, notice all of the control terms and are linear, therefore for any A1;A2 2N(x;t), A1 + (1 )A2 2N(x;t) for 2[0;1]. This shows N(x;t) is convex. 48 For the third condition, consider the rst equation of the IVP system (3.4): _T = rT(1 T K) cNT D(S;T)T kT(1 e TM)T rT: Solve the corresponding equality, we get eT = T0ert. By the Kamke Comparison Theorem, T(t) T0ert. Therefore, T0ertf is an upper bound of T(t). We denote it as T . Now, consider the system of x+ = T+ N+ C+ S+ I+ M+ T : (T+)0 = rT+ (4.26) (N+)0 = eC+ +pNN+ (4.27) (C+)0 = (4.28) (S+)0 = a1N+T +a2C+T + (pI +pT)S+ + SvS(t) (4.29) (I+)0 = C+ +pSS+ + IvI(t) (4.30) (M+)0 = MvM(t) (4.31) Compare with the system (3.4), it?s easy to see that (x+)0 x0. Therefore, x+ is a superso- lution of the system (3.4). We only need to prove it is bounded, which can be shown by the following steps. Notice that the time interval is [t0;tf], which is nite. 1. By (4.28), C+ is bounded above. Say an upper bound is C . Also, T+ is bounded above by T . M+ is also bounded above because of (4.25). Let M denote an upper bound of it. 2. Combine the above results with (4.27), we get (N+)0 Constant + pnN+, therefore N+ is bounded above by, say, N . 49 3. Combine the results from above two steps with (4.29), together with the fact that vS(t) is nite, we get (S+)0 Constant+ (pI +pT)S+. Then S+ could be bounded by, say, S . Similarly, I+ could be bounded by I . Thus, we nd a constant > 0 such that kxk>> >>< >>>> >: ("S; S); for vS(t); ("I; I); for vI(t); ("M; M); for vM(t): Thus, we get (4.34){(4.36). The characterizations are given by the de nition. 54 Chapter 5 Conclusion and Discussion We have studied four di erent ODE systems about cancer therapies. For the rst three models, we focus on the qualitative analysis of their dynamical properties; for the last one, we investigate the optimal control for the whole system for the purpose of improving the treatment protocol. There are several further subjects that we can look at in the future: 1. For the optimal control, we can further look at higher order representation of the control. Also, the competition between di erent populations may need some further, more detailed assumptions, so as to deal with di erent kinds of tumor cells. 2. All these four systems are built by ODEs. To study and understand the interrelation- ship more speci cally, we can employ the tool of partial di erential equations (PDE). A summary of some work that has been done in this area can be found in [37]. 55 Bibliography [1] G. Belostotski and H.I. Freedman. A control theory model for cancer treatment by radiotherapy. Int. J. Pure Appl. Math, 25:447{480, 2005. [2] L.G. Hanin, L.V. Pavlova, and A.Y. Yakovlev. Biomathematical problems in optimiza- tion of cancer radiotherapy. CRC, 1993. [3] R.K. Sachs, P.L. Chen, P.J. Hahnfeldt, and L.R. Hlatky. Dna damage caused by ionizing radiation. Mathematical biosciences, 112(2):271{303, 1992. [4] H. Schollnberger. Adaptive response and dose-response plateaus for initiation in a state- vector model of carcinogenesis. International journal of radiation biology, 75(3):351{364, 1999. [5] J.A. Nickolo and M.F. Hoekstra. DNA Damage and Repair: Advances from phage to humans. Contemporary cancer research. Humana Press, 2001. [6] H.D. Thames and J.H. Hendry. Fractionation in radiotherapy. Taylor & Francis, 1987. [7] R.T. Skeel. Handbook of Cancer Chemotherapy. Number v. 236 in Handbook of Cancer Chemotherapy. Lippincott Williams & Wilkins, 2007. [8] H.E. Skipper, F.M. Schabel Jr, and W.S. Wilcox. Experimental evaluation of poten- tial anticancer agents. xiii. on the criteria and kinetics associated with" curability" of experimental leukemia. Cancer chemotherapy reports. Part 1, 35:1, 1964. [9] L. Norton and R. Simon. Tumor size, sensitivity to therapy, and design of treatment schedules. Cancer treatment reports, 61(7):1307, 1977. [10] L. Norton and R. Simon. The norton-simon hypothesis revisited. Cancer treatment reports, 70(1):163, 1986. [11] N.H. Holford and L.B. Sheiner. Pharmacokinetic and pharmacodynamic modeling in vivo. Critical reviews in bioengineering, 5(4):273, 1981. [12] K.R. Fister and J.H. Donnelly. Immunotherapy: an optimal control theory approach. Mathematical biosciences and engineering: MBE, 2(3):499, 2005. [13] C. Chan, A.J.T. George, and J. Stark. T cell sensitivity and speci city-kinetic proof- reading revisited. Discrete and Continuous Dynamical Systems Series B, 3(3):343{360, 2003. 56 [14] L.G. de Pillis, W. Gu, and A. Radunskaya. Mixed immunotherapy and chemotherapy of tumors: modeling, applications and biological interpretations. Journal of theoretical biology, 238(4):841{862, 2006. [15] S. T. R. Pinho, H. I. Freedman, and F. Nani. A chemotherapy model for the treatment of cancer with metastasis. Mathematical and Computer Modelling, 36(7-8):773 { 803, 2002. [16] L.G. de Pillis and A. Radunskaya. A mathematical model of immune response to tumor invasion. Computational Fluid and Solid Mechanics, 2003. [17] K.R. Fister and J.C. Panetta. Optimal control applied to competing chemotherapeutic cell-kill strategies. SIAM Journal on Applied Mathematics, 63(6):1954{1971, 2003. [18] H.I. Freedman and S.T.R. Pinho. Persistence and extinction in a mathematical model of cell populations a ected by radiation. Periodica Mathematica Hungarica, 56(1):25{35, 2008. [19] H.I. Freedman and S.T.R. Pinho. Stability criteria for the cure state in a cancer model with radiation treatment. Nonlinear Analysis: Real World Applications, 10(5):2709{ 2715, 2009. [20] L. G. de Pillis and A. Radunskaya. A mathematical tumor model with immune resistance and drug therapy: an optimal control approach. Journal of Theoretical Medicine, 3:79{ 100, 2000. [21] L.G. de Pillis, K.R. Fister, W. Gu, C. Collins, M. Daub, D. Gross, J. Moore, and B. Preskill. Seeking bang-bang solutions of mixed immuno-chemotherapy of tumors. Electronic Journal of Di erential Equations, 2007(171):1{24, 2007. [22] A. Illanes and S.B. Nadler. Hyperspaces: fundamentals and recent advances. Mono- graphs and textbooks in pure and applied mathematics. M. Dekker, 1999. [23] Konstantin Mischaikow, Hal Smith, and Horst R. Thieme. Asymptotically autonomous semi ows: chain recurrence and Lyapunov functions. Trans. Amer. Math. Soc., 347(5):1669{1685, 1995. [24] C.C. Chicone. Ordinary di erential equations with applications. Texts in applied math- ematics. Springer, 2006. [25] J.K. Hale. Ordinary Di erential Equations. Dover Publications, 2009. [26] L. Perko. Di erential equations and dynamical systems. Texts in applied mathematics. Springer, 2001. [27] W. Walter. Ordinary di erential equations. Graduate texts in mathematics. [28] D.L. Lukes. Di erential equations: classical to controlled. Mathematics in science and engineering. 57 [29] R.F. Hartl, S.P. Sethi, and R.G. Vickson. A survey of the maximum principles for optimal control problems with state constraints. SIAM review, 37(2):181{218, 1995. [30] M.I. Kamien and N.L. Schwartz. Dynamic optimization: the calculus of variations and optimal control in economics and management. Advanced textbooks in economics. [31] Vladimir Kuznetsov, Iliya Makalkin, Mark Taylor, and Alan Perelson. Nonlinear dy- namics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bulletin of Mathematical Biology, 56:295{321, 1994. [32] Denise Kirschner and John Carl Panetta. Modeling immunotherapy of the tumor im- mune interaction. Journal of Mathematical Biology, 37:235{252, 1998. [33] L.G. de Pillis and A. Radunskaya. A mathematical model of immune response to tumor invasion. Computational Fluid and Solid Mechanics, 2:1661{1668, 2003. [34] K.R. Fister and J.H. Donnelly. Immunotherapy: an optimal control theory approach. Mathematical biosciences and engineering: MBE, 2(3):499, 2005. [35] T.N. Burden, J. Ernstberger, and K.R. Fister. Optimal control applied to immunother- apy. Discrete and Continuous Dynamical Systems Series B, 4(1):135{146, 2004. [36] D.E. Kirk. Optimal control theory: an introduction. Dover books on engineering. Dover Publications, 2004. [37] P.S. Kim, D. Levy, and P.P. Lee. Modeling and simulation of the immune system as a self-regulating network. Methods in enzymology, 467:79{109, 2009. [38] J.R. Munkres. Topology. Prentice Hall, 2000. [39] L.M. Hocking. Optimal control: an introduction to the theory with applications. Oxford applied mathematics and computing science series. Clarendon Press, 1991. [40] R.V. Gamkrelidze. Principles of optimal control theory. Mathematical concepts and methods in science and engineering. Plenum Press, 1978. [41] G. Leitmann. The calculus of variations and optimal control: an introduction. Mathe- matical concepts and methods in science and engineering. Plenum Press, 1981. [42] J. Zabczyk. Mathematical control theory: an introduction. Modern Birkhauser Classics. Birkhauser, 2007. [43] S. Lenhart and J.T. Workman. Optimal control applied to biological models. Chapman and Hall/CRC mathematical & computational biology series. Chapman & Hall/CRC, 2007. 58 Appendices 59 Appendix A One-point compacti cation and associated results De nitions: A topological space X is said to be locally compact at x if there is some compact subspace C of X that contains a neighborhood of x. If X is locally compact at each of its points, X is said to be locally compact. Theorem [38]: Let X be a topological space. Then X is locally compact Hausdor if and only if there exists a space Y satisfying the following conditions: 1. X is a subspace of Y. 2. The set Y X consists of a single point. 3. Y is compact Hausdor space. De nitions: If Y is a compact Hausdor space and X is a proper subspace of Y whose closure equals Y, then Y is said to be a compacti cation of X. If Y X equals a single point, then Y is called the one-point compacti cation of X. 60 Appendix B Fundamental knowledge of optimal control theory 1. General formulation [39] [40] We consider the following di erential equation in Rn, _x = f(t;x;u); (B.1) The point x = (x1; ;xn) 2Rn will be called the state vector (variables), and the parameter u = (u1; ;ur)2Rr will be called the control vector (variables). We call (B.1) the state equation, where f = f(t;x;u) maps from R1+n+r to Rn. Any absolutely continuous function x(t);t2 [t1;t2] is called a solution of this equation on the interval [t1;t2] if it satis es (B.1) for almost all t2[t1;t2]. 2. Hypotheses [39] [40] In order to assure that the state equation possesses a unique solution for a given initial value, we introduce certain hypotheses: For the state equation, we assume that f is continuous on its domain and has a continuous derivative with respect to x: @f @x = @f @x1; ; @f @xn : Now we introduce conditions for control vectors. We call a function x(t) : [t0;t1]!Rr a control vector function, where t0;t1 2R+[f0g;t0 < >: a if (t) < 0; ? if (t) = 0; b if (t) > 0: If = 0 cannot be sustained over a time interval, but occurs only at nitely many points, then the control u is bang-bang. If 0 on some interval of time, we say u is singular on that interval. A character- ization of u on this interval must be found using other information. 64