Modeling of Tumor Growth Incorporating the Effects of Necrosis and the Effect of Bevacizumab

Tumor growth models are important to create an engineering background for cancer treatment either by using the models for simulations and evaluation of treatment protocols or, if combined with control engineering, by designing treatment protocols. A well-defined tumor growth model must describe the physiological processes and the measurements as well. Growing tumors are composed of dead tumor cells (forming the necrotic part) and living, proliferating tumor cells (forming the proliferating part); when tumor volume is measured, these parts are measured together. Most of the known tumor growth models do not consider the modeling of the necrotic part. Starting from a minimal model of the tumor growth under bevacizumab treatment, the aim of the current research is to extend it incorporating the volume and dynamics of the necrotic part and the pharmacodynamics and mixed-order pharmacokinetics of the applied drug.The extendedmodel is validated usingmeasurements with mice as hosts, colon adenocarcinoma as tumor, and bevacizumab as the drug used for treatment.The results show that the extendedmodel can describe the important physiological phenomena and shows a good fit to the average of the measurements.


Introduction
Angiogenesis, the formation of new blood vessels, is a fundamental process required for the growth of primary tumors [1].Inhibition of angiogenesis is thus a promising way of fighting against cancer [2,3].However, the optimal protocol for antiangiogenic treatment is still under research in clinical practice [4].
Application of control theory has been considered in the literature to give a solution to the dosage problem of antiangiogenic drugs by many authors; see, for example, [5][6][7][8][9][10].Nevertheless, control engineering methods can only be applied sufficiently if there is a reliable model of tumor growth that incorporates the effect of the drug.Most of the closedloop control approaches in the literature are based on the Hahnfeldt model [11] or some of its modified versions, for example, [12].This model defines the dynamics of the tumor and the vasculature using a second-order model (without pharmacokinetics), described by nonlinear differential equations.
In [13], it was shown that tumor dynamics can be described by using a simple, first-order model (without pharmacokinetics) containing a linear and a bilinear term, where the bilinear term defines the effect of the drug on the tumor dynamics and ensures positivity of the system.This model has the ability to explain the experiments in which mice received one injection and almost every fundamental physiological process behind tumor dynamics.This minimal model, which is detailed in Section 2, seems to be a promising model for the applications, since it has a relatively simple structure (as opposed to the Hahnfeldt model, its differential equations contain linear terms and only one bilinear term) and can be used for controller design purposes as well [14].
However, the model does not incorporate the phenomenon of tumor cell necrosis which is a considerable process happening in growing tumors.Tumors have a necrotic part, composed of dead tumor cells, but the effect of necrotic cells on tumor growth dynamics is still under investigation.On the one hand, necrotic cells usually evoke inflammatory 2 Complexity response which may lead to tumor regression and hence it could be used for cancer therapy.On the other hand, necrotic cells stimulate proliferation and angiogenesis (via tumor necrosis factor) in most cancer cells, resulting in tumor promotion.As a consequence, a necrotic region could have either pro-or antitumor effect [15,16].These phenomena are not known exactly yet, but we can clearly state that necrotic tumor cells consume space inside the tumor, so when tumor measurement is done, we measure the necrotic and proliferating tumor cells as well.However, it is important to separate necrotic tumor volume from the living tumor volume.Thus, in order to have a more reliable tumor model, we need to incorporate the effect of tumor necrosis as well.
Moreover, most tumor growth models do not consider the pharmacodynamics of the drug, that is, the phenomenon increasing the drug dosage does not necessarily result in proportional increase in the effect of the drug.The effect of the drug is considered with Michaelis-Menten kinetics, and we use mixed-order pharmacokinetic model to describe the depletion of the drug.In [17] it was shown that the therapy is much more effective if we give small doses of the drug each day as opposed to giving a much larger amount at the beginning of the therapy, regardless of the fact that the drug depletes slowly.Using the mixed-order pharmacokinetics and the pharmacodynamics in our model, we are able to describe this phenomenon that could not be described by the minimal model in [13].
In Section 3, we extend the minimal model to incorporate the effect of tumor necrosis and give a third-order model to describe the dynamics of proliferating and necrotic tumor cells and use mixed-order pharmacokinetics and standard pharmacodynamics to describe the dynamics of the applied inhibitor.The output of the new model is the sum of the proliferating and necrotic tumor volumes; we use this output for parametric identification in Section 4.
We use measurements from experiments with C57Bl/6 mice, using C38 adenocarcinoma as tumor and bevacizumab as inhibitor [17].C38 colon adenocarcinoma is a well-known and widely used mouse tumor, which is originated from columnar epithelium of colon's mucosa and has the following advantages: (i) it grows fast in mice (after 2-3 weeks it reaches a lethal size); (ii) due to its specificity, there is no need to use immunosuppressed mice; (iii) a piece of tumor can be implanted subcutaneously into the mouse; (iv) it has large relative vascular area; (v) it typically does not metastasize; (vi) tumor cells inflict strong hypoxial reaction.
These properties of the tumor used in the experiments were favorable for our intentions of creating a simple but descriptive model of tumor growth under the effect of angiogenic inhibition.
The results show that the proposed model is capable of explaining the measurements (shows a good fit for the average of the measurements) while modeling the most important physiological properties of tumor dynamics as well.

Minimal Tumor Growth Model
A minimal tumor growth model under bevacizumab treatment was proposed in [13] and is further modified to create a more realistic tumor growth model in the following sections.The minimal tumor growth model can be explained with an analogy to chemical reactions.Suppose that the species  1 represents tumor volume, the species  2 represents the inhibitor level, while  denotes a compartment outside of the model (we use it to denote inflows and outflows).Then the model can be regarded as a fictive chemical reaction with the following equations: The last equation expresses the inhibitory effect of the applied drug on tumor growth.Note that this equation defines direct effect of the drug on the tumor volume; however in the physiological process, the effect is indirect.This will be detailed in Section 5. Without proper treatment, the tumor creates angiogenic signaling molecules, and as a result, the host body creates dense vasculature around the tumor, and the tumor grows into that vasculature; thus it has the required sources of nutrients for proliferation [1].If the angiogenesis is inhibited, for example, with the usage of angiogenic inhibitors (like bevacizumab), the vasculature will be less dense; thus the tumor growth rate will decrease.
The differential equations of the model can be created based on the chemical reaction equations with techniques, for example, from [18][19][20] using mass-action kinetics.The resulting differential equations and the output of the system are where  1 denotes the function of tumor volume in mm 3 ,  2 denotes the function of inhibitor level in mg/kg,  is the output function (that we can measure) in mm 3 , and  is the input function, the rate of inhibitor injection measured in mg/(kg⋅day).The model has three parameters,  is the tumor growth rate in 1/day,  is the inhibition rate in kg/(mg⋅day), and  is the clearance of the inhibitor in 1/day.The clearance of intravenous bevacizumab injection is  = ln 2/3.9 1/day acquired from [21].This parameter was not tuned for the measurements in [13,17]; however the other parameters were identified from mice experiments in [13], and their values are  = 0.27 1/day and  = 0.0074 kg/(mg⋅day).
This model is a very simple one (containing only one bilinear term, while all the other terms are linear) that describes tumor growth under treatment with angiogenic inhibitor.The novelty of the model is the differential equation of the tumor growth that contains a linear and a bilinear term.The differential equation of the inhibitor is based on the firstorder pharmacokinetic equation of the drug and can be found in standard literature; see, for example, [11].The differential equations define a positive system, that is, if all the initial conditions are positive, then the solutions of the differential equations are positive as well.

Extended Tumor Growth Model
The tumor growth model discussed in the previous section incorporates the fundamental physiological processes involved in the evolution of tumor volume, that is, the proliferation of the tumor and the inhibitory effect of the drug.It was shown in [13] using parametric identification and validation based on measurements from mice experiments that the model is also capable of explaining some of the experimental results as well with some limitations.
There is a fundamental physiological process that is not modeled in the minimal model: the necrosis of tumor cells.Necrosis appears in every tumor [1]; thus the measured tumor volume is the total volume of the proliferating (living) and necrotic (dead) tumor cells.Since necrotic cells have important effect on tumor growth; this implies that modeling the necrotic tumor volume is crucial if we want to create a valid tumor growth model.
Moreover, the minimal model in [13] is not able to explain an important phenomenon recognized in [17].In [17], it was shown statistically that frequent, small doses of drug (1/18 mg/kg each day, equivalent to 9.5 ⋅ 10 −4 mg/ml serum level each day), have significantly larger effect on tumor growth than one large dose (10 mg/kg, equivalent to 0.171 mg/ml serum level) at the beginning of the treatment.The model in [13] was not able to reproduce this result; thus the model here is modified such that the pharmacodynamics of the drug is incorporated into the model.
We extend the previously defined model by adding the dynamics of the necrotic tumor volume and the pharmacodynamics of the drug.Similar to Section 2, we give the chemical equivalents of the differential equations with the following notations: the species  1 represents the proliferating tumor volume, the species  2 represents the necrotic tumor volume, and the species  3 represents the inhibitor serum level.The equations of the model are as follows: that defines that there is an outflow of the inhibitor with a reaction rate coefficient , that is, the clearance of the inhibitor.We use Michaelis-Menten kinetics in order to have a mixed-order model for the pharmacokinetics, so this equation results in the term ẋ 3 = − 3 /(  +  3 ), where the parameter   is the Michaelis-Menten constant of the inhibitor; (iv)  1 +  3   →  2 that defines that the inhibitor binds to the angiogenic signaling molecules whose concentration is proportional to the proliferating tumor volume; as a result, formation of blood vessels is inhibited; thus growth of the proliferating tumor cells is inhibited, and necrosis takes place.The effect of the inhibition is considered with Michaelis-Menten kinetics with Michaelis-Menten constant ED 50 (called the median effective dose [22]) resulting in the velocity term  1  3 /(ED 50 +  3 ).This inhibitory effect on the volumes is considered with reaction rate coefficient .The effect of this equation on the dynamics of the proliferating and necrotic tumor volumes is expressed by the terms ẋ 1 = − 1  3 /(ED 50 +  3 ) and ẋ 2 =  1  3 /(ED 50 +  3 ).Since these terms have the dimension mm 3 /day, these terms can not be directly used to modify the dynamics of the inhibitor serum level, since that has the dimension mg/(ml⋅day).Thus, we use the constant  with dimension mg/(ml⋅mm 3 ) to define the term ẋ 3 = − 1  3 /(ED 50 +  3 ).Instead of , we will use the constant   =  in the remaining of the paper.
The combination of these terms give the differential equation of the extended tumor growth model: where  1 is the time function of proliferating tumor volume in mm 3 ,  2 is the time function of necrotic tumor volume in mm 3 ,  3 is the time function of inhibitor serum level in mg/ml, and  is the input that is the time function of inhibitor injection rate in mg/(ml⋅day).
The parameters of the model are listed in Table 1 where their notations, names, and dimensions are given.The table also contains the values of the parameters that resulted after identification which will be detailed in Section 4.
The output  of the system is the measured tumor volume in mm 3 that is the sum of the proliferating ( 1 ) and necrotic ( 2 ) tumor volumes; that is, The dynamics of the output is described by the differential equation that is the sum of ( 4) and ( 5); thus the change of the measured tumor volume depends only on the tumor growth rate constant  and the actual volume of the proliferating tumor volume.
We will use the method discussed in [23] to show that this model is nonnegative; that is, for nonnegative initial conditions and nonnegative input, the solution of the differential equation is always nonnegative.It was shown in [23, Proposition 2.1] that the dynamical system is nonnegative if and only if the vector field corresponding to the right-hand side of the differential equation of the dynamical system is essentially nonnegative.A vector field is said to be essentially nonnegative if its th component is nonnegative whenever the th state variable is zero and all other states are nonnegative.The vector field defining the dynamics of the system given by the right-hand sides of (4)-( 6) is essentially nonnegative, since Thus the system with dynamics defined by the differential equations ( 4)-( 6) is nonnegative.Note that here we have supposed that  ≥ 0 which is based on the physiological consideration that we can not take out drug from the patient, so the input can not be negative in practice; however it could be negative theoretically.

Results of the Parametric Identification
The parametric identification has been carried out using experiments made with mice, for the details of the experiment see [17].The C57Bl/6 mice with C38 colon adenocarcinoma were treated with bevacizumab using two different therapies.
(i) Therapy 1. Five mice (mouse C1-C5) got an injection of 0.171 mg/ml bevacizumab at the beginning of the treatment (Figure 1(a)).
The parametric identification was done using a combination of random and gradient-based search algorithm with more than 100000 runs that aimed to minimize the distance of the simulated total tumor volume and the measured tumor volumes in the least squares sense.The parameters that resulted in the smallest least squares error are presented in Table 1.Note that the median effective dose (ED 50 ) of the inhibitor was not identified; this value is given for bevacizumab in [22].The initial proliferating tumor volume  1 (0) was identified as a parameter; this value is the initial value for the proliferating tumor volume in the simulations, while we suppose that initially the necrotic part has zero volume.We use this assumption since in the experiments morphologically homogeneous cancer cells were implanted subcutaneously that have no necrotic part [24].
The simulated total tumor volumes for Therapy 1 and for Therapy 2 are shown along with the measurements and their average in Figures 1(a) and 1(b), respectively.The simulated therapies using the model with the identified parameters show a good fit for the measurements in the least squares sense; they are close to the average of the measurements.Since there were nine mice for Therapy 2 and only five for Therapy 1, but the measurements were used with the same weight in the identification process, the resulting outcome is closer to the average of the measurements from mice E1-E9 (that received Therapy 2).Note that the average is only valid at the time of measurements (denoted by thick dots); the dotted curve is just a linear interpolation and is used for visual clarity.It should be emphasized that the quantities shown in Figures 1(a) and 1(b) are the total tumor volumes, that is, the sum of the proliferating and necrotic tumor volumes ( 1 +  2 ). Figure 3(a) shows the inhibitor serum level as a result of the simulation of Therapy 1.The inhibitor serum level is initially 0.171 mg/ml and decreases and depletes in about two weeks.Figure 3(b) shows the inhibitor serum level as a result of the simulation of Therapy 2. The inhibitor serum level is initially 9.5 ⋅ 10 −4 mg/ml, decreases on a certain time interval, but gets refilled each day due to the daily injections.Initially, the mean of the inhibitor serum level increases, but it starts to decrease as the proliferating tumor volume grows.The simulated final tumor volumes and total amount of used drugs for the two therapies are shown in Table 2.The simulation results show that the application of Therapy 2 results in much lower tumor volume than the application of Therapy 1 using ten times less drug, as it was shown in [17].

Discussion
A key factor of angiogenesis is VEGF (vascular endothelial growth factor).Bevacizumab inhibits VEGF-induced proliferation of endothelial cells [25].The clearance of bevacizumab is a complex process which depends on several factors, namely on body weight, gender, serum albumin, alkaline phosphatase (ALP), and serum aspartate aminotransferase (AST) [26,27].In our previous mice experiments [17] all mice were 20 gram male C57Bl/6 mice.Since our aim was to create the most compact model possible, it is out of our scope to incorporate effects like serum albumin, ALP, and AST changes.However, the clearance of bevacizumab also depends on the tumor burden.Patients with a lower tumor burden have a lower clearance rate [28,29] which was important to be modeled.The estimated clearance of bevacizumab is 0.207 1/day (95% CI, 0.188-0.2261/day) [26].
The clearance parameter was identified as 0.1825 1/day in the extended tumor growth model.
Another notable phenomenon is the metabolic profile of bevacizumab that was found to be similar to a native IgG (Immunoglobulin G antibody) molecule which does not bind VEGF [30].The major histocompatibility complex (MHC) class I-related receptor FcRn (neonatal Fc Receptor), which can be found in the endothelium of small arterioles and capillaries, modulates IgG trafficking across tissues [31] and hence plays a major role in the clearance of bevacizumab as well.Serum bevacizumab has two depletion ways (Figure 4).The evident process is the binding to cell surface receptors; this is the drug's mechanism of action, namely, inhibition of angiogenesis, leading to necrotic tumor cells.This process was modeled with Michaelis-Menten kinetics using median effective dose ED 50 , resulting in the term The reaction rate of binding to cell surface receptors is  (  = ).The other depletion way is the bulk phase nonspecific endocytosis [32] which is the uptake of extracellular bevacizumab via pinocytosis into endosomes of catabolic cells where it binds to FcRn [33].This binding protects bevacizumab from degradation and systemic elimination and indicates recyclation to the serum, accounting for the longer half-life of the angiogenic inhibitor.This process was modeled with Michaelis-Menten mixedorder model, resulting in the term ẋ 3 = − 3 /(  +  3 ).The reaction rate of internalization into endosome is .
Cancers are not composed of homogeneous cell groups but rather contain different subpopulations of tumor cells.There are active, proliferating cells (typically these cells are located on the surface of the tumor or near living/nascent vessels) but inactive cells are present as well.Most of the inactive cells have died, for example, due to hypoxia; these cells form necrotic regions; however, there are special types of tumor cells which are only temporarily inactive and have propagating potential.These quiescent cells are also known as cancer stem or stem-like cells.Cancer stem cells are in a reversible G0 phase from which these cells may escape on account of specific signaling pathways and molecules like tumor protein p53, retinoblastoma protein (RB), and cyclindependent kinases inhibitors (CDKIs) [34].The quiescent cells can disseminate to distant locations causing tumor dormancy even for long time periods.In addition, stem cells also exist within the tumor leading to the ability of resistance to conventional chemo-and radiotherapies [35].In our experiment, C38 colon adenocarcinoma was used as tumorous cell line which typically does not metastasize; hence modeling disseminated quiescent cells would be inexpedient.Quiescent cells within the primer tumor have critical role in conventional therapies (like chemotherapy), but their role in antiangiogenic monotherapy is not clear.Cai et al. [36] modeled the tumor growth and its microenvironment in antiangiogenic therapy assuming four different tumor cell phenotypes: proliferative cells, quiescent cells, necrotic cells, and migrating cells.Their criterion for a quiescent cell is the case when the tumor cell "satisfies the survival condition but there is no neighboring space for it to proliferate" [36].Upon this unfavorable space condition ceases (and the oxygen supply is sufficient), the quiescent cell will turn back into a proliferating cell.In our experiment, due to the subcutaneous localization of the tumors, these limiting space conditions are not presumptive; therefore our extended tumor model does not incorporate the effect of quiescent cells; we take into account these cells as proliferating cells.
In the extended tumor growth model, the effect of angiogenic inhibitor is taken into account as a direct effect on tumor cells, instead of an indirect effect through the tumor vasculature.The aim of this modeling approach is manifold.On the one hand, our aim was to model tumor growth taking into account the most important and separatable effects.Previously we have found that vasculature volume changes and tumor volume changes are hardly separatable in time [37,38], which implies that identification of the parameters corresponding to vasculature dynamics is very hard if measurements of tumor vasculature are not available.This would limit the future application of the tumor growth model in model-based individual therapy, since the identification of the parameters specific to a patient would require measurement of the vasculature of the patient's tumor.
Here we have proposed a descriptive model that does not require vasculature measurements for identification.On the other hand, it was found that VEGF inhibition may have direct effects on tumor cells that impair tumor growth and metastasis [39].Inhibiting vascular endothelial growth factor receptor 1 (VEGFR1) blocked tumor cell migration, invasion, and colony formation in human colon cancer cells [40].
Due to the pharmacodynamics of the drug, the effect of the inhibitor on the tumor growth is limited.In order to analyze this effect, we examine two extreme cases.
(1) There is no drug present, that is,  3 ≡ 0. In this case the proliferating tumor volume dynamics is described by and the solution to that differential equation for any Using the parameter values from Table 1, the proliferating tumor volume for positive times without treatment is given by (2) The serum level of the drug is high; consider the case  3 → ∞.Due to the pharmacodynamics, the value of the term in (4) describing the effect of the drug is lim and thus the proliferating tumor volume dynamics is described by the differential equation and the solution of this equation for  ≥ 0 is Using the parameter values from Table 1, the proliferating tumor volume for positive times with infinite serum level of the drug is given by This implies that, for any positive time  and positive serum level  3 (), the proliferating tumor volume  1 () is bounded as The solution to (4) can be approximated as due to the small value of   (thus in our approximation we have neglected the effect of  1 on  3 ).This shows that the tumor growth is described by an exponential curve and the inhibitor modifies the rate of the growth.As a consequence of ( 16), the applied drug can not stop the growth of the tumor; however, it can decrease the tumor growth rate resulting in longer survival.The tumor growth can be stopped if and only if the parameters satisfy the inequality The values of the parameters in Table 1 do not satisfy this inequality; thus the tumor growth can not be stopped.This result implicates that in theory the tumor growth can be stopped by

Conclusions and Future Works
We have shown that the extended model is able to describe the most relevant physiological phenomena and explains the measurements as well.The model structure was kept as simple as possible; most of the terms in the differential equations are linear, except for the pharmacodynamics and pharmacokinetics; however these nonlinear terms were necessary to describe an important phenomenon recognized in [17] that is, that much smaller, but more frequent doses are more effective than a great initial dose, regardless of the long depletion time of the drug.
The simplicity however has a price: the dynamics of the vasculature is not modeled.For a good model of vasculature dynamics and the effect of inhibition on the tumor vasculature, see [42].Nevertheless, the results have shown that the extended model can describe the tumor growth dynamics efficiently even without modeling the vasculature; thus we have neglected the vasculature dynamics (note that modeling vasculature dynamics would increase the order of the system; thus the model would be more complicated).Moreover, experiments done with untreated tumor in [38] have shown that tumor dynamics can be described with a first-order model, implicating that tumor growth dynamics dominates vasculature dynamics, which also implies that identification of the parameters corresponding to vasculature dynamics is a hard task or nearly impossible without having measurements of the vasculature.Nevertheless, in our experiments, these measurements were not available.
The structure of the model and the identified parameters imply that the tumor growth can not be stopped using bevacizumab as monotherapy; however, it can slow down the tumor growth, as it turned out in clinical practice as well [22].Based on the model, tumor growth can be stopped by either finding a more effective drug or combining bevacizumab treatment with other treatments.This is recognized in clinical practice as well, since bevacizumab in combination therapy is widely used recently [3].
The main goal of the presented work is to model the effect of bevacizumab on tumor growth; however the long-term goal is to model the effect of therapies where bevacizumab is combined with other drugs.Since modeling the effect of combined therapies is challenging in many ways (identification, experiment design, etc.), analyzing and modeling the effect of bevacizumab only are beneficial; in this way its effect is separated from the other drugs.In our future works, we will extend the model presented here to be able to describe mixed therapies.The presented model structure gives a good basis for this extension.

Figure 2 (
a) shows the total tumor volume and necrotic tumor volume as a result of the simulation of Therapy 1.The final total tumor volume is 4741 mm 3 , with 1529 mm 3 necrotic part.Figure 2(b) shows the total tumor volume and necrotic tumor volume as a result of the simulation of Therapy 2. The final tumor volume in this case is 3713 mm 3 with 1931 mm 3 necrotic part.

Figure 1 :Figure 2 :Figure 3 : 6 ComplexityTable 2 : 2 Figure 4 :
Figure 1: (a) The measured tumor volumes for mice C1-C5 that got Therapy 1 and their average, and the simulated tumor volume using Therapy 1 with the identified model (solid curve).(b) The measured tumor volumes for mice E1-E9 that got Therapy 2 and their average, and the simulated tumor volume using Therapy 2 with the identified model (solid curve).
Note that this necrosis is independent of the treatment.Using mass-action kinetics, this equation modifies the dynamics of the proliferating and necrotic tumor volumes with the terms ẋ 1 1 that defines that the tumor cells proliferate (divide) with a tumor growth rate .Using massaction kinetics, this equation results in the term ẋ 1 =  1 ; (ii)  1   →  2 that defines the necrosis (death) of tumor cells with necrosis rate .

Table 1 :
The parameters of the extended tumor growth model and their values after the identification.