Mathematical Modelling of the Dynamics of Tumor Growth and its Optimal Control

The dynamics of tumor cell growth and its treatment process is discussed in this paper. We analyze some simple mathematical models and generalized models to understand the growth of tumor cells. A couple of diffusion models are discussed to explain how the tumor cells spread and become more dangerous as well as the treatment process of cancer and how cancer cell behaves in the presence of different therapy and drugs. The optimal control of chemotherapy has been discussed. It has also been explained how much the model is effective in reducing tumor cells over time.


Introduction
Mathematical modelling has been playing a great role in exploring a variety of information not only in different branches of physical and engineering sciences, but also in biological and medical sciences. A number of theoretical studies have been performed by researchers on different topics of life sciences. Among others, during the last five decades Misra and his co-workers have carried out extensive studies on a variety of topics in Biomedical Mathematics and Physiological Fluid Dynamics (e.g. (Maiti & Misra, 2011;Mallick et al. , 2019;Misra & Chakravarty, 1982;Misra & Dravid, 2006;Misra & Mitra, 2006, 2008Misra & Pandey, 2001;Misra & Singh, 1985;Misra et al. , 2004Misra et al. , , 2010Misra et al. , , 2011Misra et al. , , 2020). A mathematical model was developed (Misra et al. , 2004) to explore the enzymatic action of DNA Knots and links. Another mathematical model was formulated and analyzed (Misra & Dravid, 2006) with aim to identify transcription factor binding sites of genes. In (Misra & Mitra, 2006) performed a mathematical analysis for investigating in single-species and host-parasite system, with special reference to period-doubling bifurcations and chaos. The method of mathematical modelling was employed in (Misra & Chakravarty, 1982) to study the dynamic of arterial wall tissues and in (Misra & Singh, 1985) to determine the left ventricular wall stresses of a human-sized heart. A mathematical model was developed by Misra and Pandey (Misra & Pandey, 2001) in order to explore different aspects of the oesophageal swallowing of food bolus. The peristaltic flow of bile in a patheological state was explored in (Maiti & Misra, 2011), while non-Newtonian characteristics of blood flow in contracting/ expanding arteries were explained in  by developing different mathematical models. Recent theoretical studies reported in (Mallick et al. , 2019) and (Misra et al. , 2020) on nanofluid flow, temperature distribution and entropy generation in porous microchannel/microfluidic tubes bear promises of important applications in human microcirculatory systems. Although the mathematical models mentioned above are considered as landmark contributions in mathematical biosciences, none of them refers to growth or treatment of cancer.
In cancer research, several mathematical models have been developed. The derived pieces of information have enriched the subject oncology significantly. Knowledge of tumor growth is very important in the treatment of cancers. Studies of cancer growth by way of mathematical modelling are very useful to have a better understanding of the dynamics of tumor or cancer growth. Such studies also help evaluate different approaches. The information derived from different theoretical models enable the clinicians to decide upon the optimal cancer therapy and use it to treat a particular cancer patient (Wang, 2018). Mathematical modelling being a well-built tool to test different hypotheses to validate experimental observations and to study dynamical behaviour of complex systems can play a vital role in cancer research.
A mathematical model was formulated and analyzed to investigate synchronization among tumor-like cell aggregations coupled by quorum sensing (Misra & Mitra, 2008). In another study carried out by developing a mathematical model an attempt was made by Misra et al. (Misra et al. , 2010) to find theoretical estimate of arterial blood flow during cancer treatment by the method of electromagnetic hyperthermia. But cancer growth has not been studied in any of these communications.
Considering all the above, it has been our endeavour in this paper to discuss the following matters thoroughly: The growth dynamics of tumor cells are discussed through mathematical models.
Diffusion models are formulated to study the spread of cancer cells.
A mathematical model is also developed to investigate the optimal control in the case of chemotherapy.
The effectiveness of the model in reducing tumor cell volume has also been discussed.
It is believed that the observations reported here will be beneficial to the clinicians involved in the treatment of cancer patients and also to researchers of oncology.

Exponential Models
The growth of tumor can be modeled through first-order ordinary differential equations by examining the change of tumor volume V over time t with an initial volume V (0) = V 0 . This initial condition has been used in most tumor models. Some exponential models that can be conveniently used to describe tumor growth are discussed below.

Malthusian Model
This is the simplest model to describe tumor growth. The growth is proportional to the population of tumor cells. This model is often used to describe a single species population's growth rate. The model has the following form (Murray, 2007) : where V is the volume of the tumor, r is the growth parameter and t represents time. The solution of equation 1 is given by : The behavior of Malthusian model can be seen from the following figure: The volume of the tumor over time grows exponentially if the growth parameter r is positive. If r is negative, it decreases exponentially and if r is zero, the cell volume remains constant.

Power Law Model
The generalization of Malthusian model can be described by the power law model. It was first introduced by Mendelsohn in 1963. This model has the following form (Wang, 2018) : where r is the intrinsic growth rate. If b = 1, we get Malthusian model which we have already discussed. This equation has the solution of the following form : b−1 . Now we can draw the graph from this solution using different value of b. For b = 1, we get the Malthusian model. We consider here two cases for b > 1 and b < 1 and analyze the effect of the exponent in the model. In Figure 2, we can see that the graphs are different for b = 1.05 and b = 0.9, therefore the model behaves differently. For b = 1.05 we can observe that the volume increases faster for r > 0, and it decays at faster rate for r < 0. For b = 0.9, we see that the growth rate of tumor volume is slower than the case of the other graph. For r > 0 the volume still increases exponentially but at a slower rate and for r < 0 the volume reduces slowly. So, we can conclude that the growth rate gets faster with increase in the value b.

Migration Model
Migration model can be described by a population that obeys an exponential law of growth with migration. Migration can affect the cell volume in two different ways. It can increase the volume by adding newly affected cells or it can reduce the cell volume, which can happen when the tumor cells die. This can be described by following equation (Murray, 2007) : where r is the growth rate and the migration rate is K. When K is positive, more and more normal cells transform into tumor cells and when K is negative the tumor cells die. Now to solve the equation 3, let V = u − K r where r = 0. Then we have the solution: Plots of the solution for different migration rates are shown below : By taking the growth rate r = 10 −7 m 3 day −1 and initial population V 0 = 5 × 10 −4 m 3 we can generate the graphs in Figure 3 over time t = 0 to t = 10 days. Here we can see that if the migration rate is positive, the volume increases at a faster rate than in the case when the migration rate is negative.
For different migration rates, we can explain the behavior of the model. From the above graphs we observe that for growth rate r < 0 the cell volume V is decreasing over time and V (t) → − K r as t → ∞. Also we see that, if we take the initial volume V 0 = − K r , the system gets stuck. Although the point, V = − K r is unstable, since the system moves away when we take initial volume near this point. Again for r > 0 we observe that, the cell volume explodes and tends to ∞.

Gompertz Model
The Gompertz model exhibits an exponential decay of the growth rate. It has been successfully used to model breast and lung cancer growth. It is a sigmoid function (function that has S shaped curve) which describes growth as being slowest at start and end. It has been modified suitably for use in biology, with regard to detailing populations. The model can be described by the following form (Tsoularis & Wallace, 2002;Wang, 2018) : where a is the intrinsic growth parameter and β is the parameter of growth deceleration. Integrating both sides by taking limits from V 0 to V and t 0 to t, we get: The volume of tumor increases exponentially over time when the growth rate a > 0. Also we can see that the volume decreases exponentially over time if a < 0 and volume is constant if the growth rate is 0.

Logistic Models
The exponential models describe the tumor growth satisfactorily for a certain amount of time. But the volume tends to infinity if growth rate is positive over time which is not realistic. The volume can only increase to a certain level since there are limited resources that are necessary for cell growth and then it gets stable which can be fairly described by the logistic models. Some logistic models and their generalization are discussed here.

Logistic Model for Population Growth
The model assumes a linear decrease of the relative growth rate with population size. The maximum size is limited by a carrying capacity K. The model is described by (Murray, 2007) : where a is the coefficient of cell proliferation. We used the carrying capacity, where, b is the population deceleration rate. The solution of equation 5 is: The equation is called the logistic law of growth. We can analyze the behavior of the model from the solution graph presented below.  Figure 6: Graph of the logistic model with a = 5 × 10 −7 m 3 day −1 and K = 5 × 10 −4 .
If we take the initial volume between 0 and K/2, we can see that the volume increases exponentially over time and it gets stable at K. But if we consider the initial volume between K/2 and K the population attains stability faster. If the initial volume is greater than the carrying capacity, then cell volume decreases due to lack of nutrition and therefore the stable position K is attained quickly.

Von Bertalanffy Model
Assume that growth is proportional to surface area, since nutrient enters through the surface, and that death is proportional to the tumor size. The model is also known as the surface rule model. It has been successfully applied to describe human tumor growth (Tsoularis & Wallace, 2002).
where a is the growth parameter and b is the growth deceleration parameter. The solution of 6 is given by : We can show the behavior of the Von Bertalanffy model in the following figure. It may be observed that by taking the initial cell volume between 0 and ( a b ) 3 the cell volume increases over time and gets stable at V = ( a b ) 3 . Hence, V = ( a b ) 3 is the saturation level of the model. If we take initial volume V 0 > ( a b ) 3 the volume reduces exponentially over time and is stable at V = ( a b ) 3 . Also if we take V 0 = ( a b ) 3 the volume remains constant.

Richards' Model
This model was initially developed for modelling population growth. It can be used to discuss tumor growth also. Mathematical form of this model is given by (Tsoularis & Wallace, 2002) : Here α is the positive exponent. We can analytically derive the solution of the equation 7 by letting V α−1 K α = U and by using the partial fraction Using these in the equation 7 and integrating both sides we obtain the solution: The graphical presentation of this function is shown in Figure 8 : It can be observed that this model behaves just like the logistic model. However in this model, the system can be made stable more quickly by increasing the value of α suitably.

Generalised Logistic Model
The generalised form of the logistic equation can be written as (Tsoularis & Wallace, 2002) dV where α, β and γ are non-negative exponents and a is the growth rate parameter. From this generalised form, it is possible to derive all the models described above by choosing the values of α, β and γ suitably.
For α = 1 and γ = 0, one gets the exponential growth equation : The graph of its solution behaves exactly like the exponential growth curve that increases in an unbounded manner, as t → ∞ for a > 0.
For α = 1, β = 1 and γ = 1, equation 8 reduces to the logistic equation : The logistic curve can also be generated by setting all the exponents equal to 1. The graph is stable at V = K 2 for different initial volumes. If we take α = 2/3, β = 1/3 and γ = 1 we get from 8, which is the equation of Von Bertalanffy model. We can clearly observe that Von Bertalanffy's curve gets stable at a slower rate than the logistic curve.
Lastly, we can consider α = 1, and γ = 1. Then we get the Richard's equation The graph of its solution becomes stable at a faster speed than the normal logistic equation with the smallest increment of β.
One can now analyse the behavior of the generalised logistic model by considering different values of the parameters and observing the nature of the graphs of the solution. We notice that the model behaves differently with different parameter values. By increasing the value of α slightly one can observe that the volume of the tumor reaches the carrying capacity K quickly. Similarly, an increase in the value of β causes the volume to increase faster. But for γ the observation is to the contrary. The higher the exponent, the slower is the growth rate becomes.
Moreover differentiating 8 with respect to t and considering d 2 V dt 2 = 0, one can obtain the inflection volume Tsoularis & Wallace (2002), It may further be observed that for all three exponents α, β and γ, V inf → K as α, β → ∞ and γ → 0, that is, lim The observation presented above can be schematically shown by using the following diagram :

Generalised
Logistic Equation

Diffusion Model
Diffusion model consists a set of mathematical equations that attempts to estimate the spread of information (idea or rumor) or a contagious disease or any other phenomena through a corresponding environment. It is possible to use a diffusion model to explain how the tumor cells spread through the other cells. Let us consider the diffusion model in the following form (Bankman, 2008) : subject to the boundary condition, ∇u(t, x).n(x) = 0 where u is the tumor cell density, ∂/∂t is the partial differential operator, d the diffusion tensor for tumor cells and f (t, x, u) is the reaction term. This equation can explore two characteristics of tumor growth, viz. diffusion and proliferation. The term d∆u in 9 describes the invasion of the diffusion tensor d, while the term f (t, x, u) describes the proliferation of tumor cells. The equation 10 represents the no flux boundary condition, which means tumor cells do not diffuse towards these structures.

Diffusion Model Using Malthusian Equation
Some of the tumor cell population models have been discussed above. Using these models as the proliferation models in 9 we can discuss some diffusion models. Let us consider the Malthusian model. Then diffusion model takes the form: where r is the growth parameter. Using 10, in this equation we have : Solution of the initial boundary problem is given by: Now if we plot the solution we get the diffusion graph as below and at different values of t we can show the diffusion behavior in the following graphs: The graph shows that, over the time the diffusion of the tumor cells towards x = 5(the middle), is the highest. This implies that the cells near the middle of the tissue medium get more oxygen and nutrients and so can grow very easily. So the cells tend to shift to the middle at a high rate. When most of the tumor cells have shifted at the middle, we can observe that the overpopulation of the cells there, create lack of nutrients and energy sources and eventually they start to move away from that place. The diffusion u(t, x), however, is an increasing function of time t.

Diffusion Model Based on Gompertz Equation
Using Gompertz model as the proliferation equation 9, we can formulate another diffusion model. Then the diffusion model takes the form : where r is the growth parameter. This equation with the no flux boundary condition can be written in the following form: Solving the non-homogeneous partial differential equation 14 along with the initial condition 15 and boundary condition 16, we get the solution: The solution graphs of the diffusion behavior for different values of time t are presented below : Analysis of these graphs puts forward observations similar to those derived from the diffusion model using Malthusian growth equation for proliferating cells.

Diffusion Model Using Logistic Equation
We can describe the cell proliferation process more efficiently, if we use the logistic model in equation 9. The model can then be formulated as follows : u(0, x) = u 0 (x) (18) u(t, 0) = u(t, L) = 0 Thus we have a boundary value problem consisting of a nonlinear partial differential equation with initial and boundary conditions. We get the solution Bokhari et al. (2008) of 17, where C 1 and C 2 are arbitrary constants and I 0 and K 0 are the modified Bessel function of 1st kind of order zero and 2nd kind of order zero respectively. Now by applying the initial and boundary conditions we get the nonzero constants, and the Bessel functions we used are defined as, Now if we observe the nature of the graphs of the model for different values of t, we can examine the diffusion behavior in the following graphs: Analyzing the graphs of the model, we find almost similar result as the previous models. The cells behavior is to diffuse towards the middle. So the diffusion rate is the highest towards that. But after a while they start migrating away from the middle because of too much volume of the tumor and lack of nutrients, oxygen level and other things that are necessary for the cells to grow. The difference between the previous models and this model is that in previous models that we discussed above, the diffusion of tumor cells increased unboundedly over time. But here the diffusion u(t, x) reduces, as time progresses. This is more realistic since there are limitations of resources and they reduce over time.

Modeling Treatment
In this section, we study the controlling of tumor growth by modeling the treatments. For this purpose, we consider some cell populations, viz. tumor cells, natural killer cells, dendritic cells and cytotoxic CD8 + T cells. To model treatment, we consider all these types of cells follow the same type of equation of the form: where f (·) and g(·) are proliferation, competition and inhibition functions of the cells and K [·] , z(M )[·] are the killing cells term arising out of the effect of chemotherapy drug thus if K [·] = 0, there is no effect of drug. The term z(M ) = 1 − e −M represents the effectiveness of chemotherapy at some specific phases of cell cycles (Chang et al. , 2003;Unni & Seshaiyer, 2019).

Dynamics of Tumor Cell
Let us denote the tumor cells by T , natural killer cells by N , dendritic cells by D, cytotoxic CD8 + T cells by L and considering the logistic growth model of the tumor cells. The dynamics of tumor cells is governed by the partial differential equation, In equation 21, a denotes the intrinsic growth rate and b the growth declaration rate. So, aT (1 − bT ) is the proliferation term. Also j denotes the parameter of interaction between tumor cells and dendritic cells, c 1 denotes the parameter of interaction between tumor cells and natural killer cells and k stands for the parameter of interaction between tumor cells and CD8 + T cells.

Dynamics of Natural Killer (NK) Cells
The dynamics of NK cells is governed by the equation : where s 1 is a constant source of NK cells and the term g 1 N T 2 h 1 +T 2 represents proliferation of NK cells, g 1 being the maximum proliferation rate and h 1 the steepness coefficient. The parameter c 2 is interaction rate between tumor and NK cells, d 1 is the rate at which dendritic cells kill NK cells and e is the natural death rate of the NK cells.

Dynamics of Dendritic Cells
To examine the dynamics of dendritic cells, we may consider s 2 as a constant source of dendritic cells, f 1 the interaction rate between dendritic and CD8 + T cells, d 2 the killing parameter of dendritic cells by NK cells, d 3 the proliferation parameter of dendritic cells from tumor cells and g natural death parameter of dendritic cells. Then the dendritic cells dynamics takes place according to the equation : Dynamics of CD8 + T Cells CD8 + T cells are immune cells that kill target cells. The dynamics of this type of cells can be described by the equation : where f 2 is the activation parameter of CD8 + T cells which arises due to the interaction between tumor and dendritic cells, h is competition rate between CD8 + T and tumor cells, u is the rate at which the CD8 + T cells are eliminated due to inactiveness towards drug, r 1 is the rate at which CD8 + T cells are activated and i is the natural death rate of CD8 + T cells. The term P I LI g I +I represents activeness of CD8 + T cells due to immunotherapy, where I is the concentration of immunotherapy drug (Bankman, 2008).

Drug and Vaccine Intervention
In order to study the effects of chemotherapy and immunotherapy drugs, one can consider drug intervention in immunotherapy approach when CD8 + T cells are promoted through antigen-specific cytolytic immune cells. Introducing the term v L = v L (t) in equation 24, we have, In order to study the chemotherapy and immunotherapy drugs, we describe the dynamics of the respective concentrations in the blood stream by considering the equations : and The drug intervention terms d 4 M and d 5 I in 26 and 27 represent the amount of chemotherapy and immunotherapy drugs in some interval of time. It ay be noted that the chemotherapy and immunotherapy drugs are eliminated from the body over time at a rate proportional to its concentration (Unni & Seshaiyer, 2019).

Model Equations
As a whole the model consists of the following set of partial differential equations :

Analyzing Model's Behavior
The system of equations 28 can be solved by Runge Kutta 4 th order method. The discussion that follows is based on computational results. All the parameters in this section are estimated in cells −1 day −1 unit. As in (Unni & Seshaiyer, 2019), the effects of the existence of dendritic cells, and chemotherapy drugs, immunotherapy drugs etc.
I. We assume that values of additional activation parameters are zero for CD8 + T cells and NK cells i.e, r 1 = 0, g 1 = 0, h 1 = 0. Also let P I = g I = u = 0, so that the influence of drug killing terms (K T = K N = K D = K L = 0), no influence of drug and vaccine interventions (v L = v M = v I = 0) along with the corresponding death rates (d 4 = d 5 = 0). We will also assume d 3 = 0 which corresponds to the new term that has been added to the model to indicate the growth of dendritic cells for being impacted by tumor cells. We will also assume for simplicity and illustration purposes of a weak immune system that the dynamics start with 100 tumor cells with one natural killer, one dendritic, and one CD8 + T cell.  III. Along with d 3 , we also wanted to study the influence of the source term s 2 on the dynamics of tumor, NK and CD8 + T cells. When the source term of dendritic cells is considered, NK cells and CD8 + T cells are found to increase. It may be noted that tumor growth is effected by both NK and CD8 + T cells. Reduction of cells is quite clear from the same figure. This suggests that, an external source term of dendritic cells has the potential to reduce tumor growth. It may also be noted that for tumor growth, dendritic cells play a significant role in recruiting CD8 + T cells. IV. Next, to study the effect of drug intervention term only for the CD8 + T cell population as an immunotherapy where the immune cell levels are boosted by the addition of antigen-specific cytolytic immune cells, one can increase the value of v L from 1 to 10 6 . The effect of the chemotherapy drug can be studied by considering the influence of v M . We set v M = 1 and study the influence of increasing K T in the dynamics of tumor cells. The results are shown in Figure 20, and we notice from the first graph that the effect of adding the drug in small doses has no significant impact on tumor growth whereas the other graph shows that tumor volume can be reduced through chemotherapy.
V. Next, we turn our attention to the effect of the nonlinear term introduced as an inactivation term, that accounts for the regulation and suppression of CD8 + T cell activity. Figure 21 shows the effect of the nonlinear term in system 28 when u = 0 and u = 3 × 10 10 . Here the dynamics shows that there is a drastic drop in the number of cytotoxic CD8 + T cells, when there is a small increase in tumor cell growth.

Optimal Control
In this section, a modelis developed and analyzed with an aim to discuss different strategies to cure cancer. Here we study the optimal control for chemotherapy. This optimal control helps reduce the tumor volume and the side effects of the drugs over a given period of time. Fister and Panetta (Panetta & Fister, 2000) used an optimal technique. There are several models describing chemotheraputic killing of cancer cells. The analysis presented here is based on the skipper's log-kill hypothesis, which states that killing of cancer cells by injecting chemotheraputic drugs is proportional to the tumor at time t.
The results for this hypothesis using the Gompertzian growth model (Lenhart & Workman, 2007) have already be presented in section 1.4. If V (t) is the tumor volume at time t, the model can be described by : where r is the natural growth rate of the tumor cell, δ is the magnitude of the dose and u(t) is the effect and strength of the drug. If u(t) = 0, there is no drug effect and if u(t) > 0, the strength of the drug effect is considerable.
The objective function can be used minimize the cost of the control, the possible side-effects and the tumor volume V over a time interval. Finally, we require u(t) ≥ 0 for all t. Considering a quadratic objective function, the minimization problem can be defined as follows : where a is a positive weight parameter. If we solve this problem we can generate different graphs using different values of the parameters. The effects of drugs over time are illustrated in Figure 22. Using the values of the parameters as r = 0.3, δ = 0.45, V 0 = 0.975 and T = 20, we have compared them by changing the value of a in the same optimization problem. It is observed that volume of the tumor reduces with increase in a.
Also we can observe the results by changing values of the initial density and time periods. We compare the graphs for V 0 = 0.975 and V 0 = 0.6, also taking a longer period of time, T = 40. One can see that after a short period of time, the result do not change significantly. The reduction rate of the tumor volume is faster for first few days, then it slows down and after sometime, it becomes uniform. We can also see from the above figure that long time treatment results are similar to those for short time.
If we change the dosage δ of the drug by taking the parameters, r = 0.3, a = 3, V 0 = 0.8, T = 20 and δ = 0.3 and 0.5, we notice that for higher dosage of drugs, reduction of tumor volume is faster. For δ = 0.5 we notice much more reduction in tumor volume than for δ = 0.3. The observations of the study reveal that the time period of the treatment does not effect the reduction process and also the initial volume of tumor. The reduction process is very slow during the first few days. But the dosage of medicine plays a vital role in the reduction process. The higher the dose the faster and better is the result. In this case, tumor size is considerably reduced in the first few days and then the process slows down. Therefore we can say that the optimal control is a high drug dosage strategy to reduce cancer, which is common in cancer treatment these days.

Discussion
Modeling cancer growth and its treatment has drawn serious attention of researchers. We have discussed the behavior of tumor growth mathematically through some well known models and their generalization. We discussed the generalization of the logistic model and the diffusion models with different growth equation of proliferative cells. We have also discussed the effect of chemotherapy and immunotherapy and observed the study shows that chemotherapy is very effective in reducing the tumor cells. From this study, we can conclude that the optimal treatment of chemotherapy as a standard strategy to reduce tumor cell is highly effective in the first few days and then the reduction process slows down. This has been the observation even for longer period of treatment. It can be further concluded that with higher dose of medicine, tumor cells reduce faster.