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

Drexler, Dániel András; Sápi, Johanna; Kovács, Levente

doi:https://doi.org/10.1155/2017/5985031

Complexity

On this page

Abstract Introduction Discussion Conclusions Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2017 | Article ID 5985031 | https://doi.org/10.1155/2017/5985031

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

Dániel András Drexler,¹Johanna Sápi,¹and Levente Kovács¹

Academic Editor: David Basanta

Received25 Aug 2017

Revised25 Oct 2017

Accepted01 Nov 2017

Published26 Nov 2017

Abstract

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 extended model is validated using measurements with mice as hosts, colon adenocarcinoma as tumor, and bevacizumab as the drug used for treatment. The results show that the extended model can describe the important physiological phenomena and shows a good fit to the average of the measurements.

1. 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–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 closed-loop 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 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.

2. 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 represents tumor volume, the species 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:(i) that defines that the tumor cells proliferate (divide) with a tumor growth rate ;(ii) that defines that there is an outflow of the inhibitor with a reaction rate coefficient , that is, the clearance of the inhibitor;(iii) that defines that if there is inhibitor () present then the growth of tumor () is inhibited.

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–20] using mass-action kinetics. The resulting differential equations and the output of the system arewhere denotes the function of tumor volume in , denotes the function of inhibitor level in , is the output function (that we can measure) in , 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 , is the inhibition rate in kg/(mg·day), and is the clearance of the inhibitor in .

The clearance of intravenous bevacizumab injection is = /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 first-order 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.

3. 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 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 represents the proliferating tumor volume, the species represents the necrotic tumor volume, and the species represents the inhibitor serum level. The equations of the model are as follows:(i) that defines that the tumor cells proliferate (divide) with a tumor growth rate . Using mass-action kinetics, this equation results in the term ;(ii) that defines the necrosis (death) of tumor cells with necrosis rate . 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 , ;(iii) 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 , where the parameter is the Michaelis-Menten constant of the inhibitor;(iv) 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 (called the median effective dose [22]) resulting in the velocity term . 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 and . Since these terms have the dimension mm³/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³) to define the term . 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 is the time function of proliferating tumor volume in , is the time function of necrotic tumor volume in , 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³ that is the sum of the proliferating () and necrotic () tumor volumes; that is,

The dynamics of the output is described by the differential equationthat 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 ] 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(i); thus the first component of the vector field is nonnegative if ;(ii), if ;(iii), if .

Thus the system with dynamics defined by the differential equations (4)–(6) is nonnegative. Note that here we have supposed that 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.

4. 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)).

(a)

(b)

(ii) Therapy 2. Nine mice (mouse E1–E9) got mg/ml injection of bevacizumab each day for 18 days (Figure 1(b)).

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 () of the inhibitor was not identified; this value is given for bevacizumab in [22]. The initial proliferating tumor volume 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 ().

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 , with 1529 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 with 1931 necrotic part.

(a)

(b)

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 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.

(a)

(b)

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].

5. 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.226 1/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 , 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 mixed-order model, resulting in the term . 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 cyclin-dependent 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, . In this case the proliferating tumor volume dynamics is described by and the solution to that differential equation for any is 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 . Due to the pharmacodynamics, the value of the term in (4) describing the effect of the drug is and thus the proliferating tumor volume dynamics is described by the differential equation and the solution of this equation for is Using the parameter values from Table 1, the proliferating tumor volume for positive times with infinite serum level of the drug is given byThis implies that, for any positive time and positive serum level , the proliferating tumor volume is bounded asThe solution to (4) can be approximated asdue to the small value of (thus in our approximation we have neglected the effect of on ). 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 inequalityThe 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

(1)increasing the value of , which could be done by finding an angiogenic drug that is more effective;(2)increasing the value of , which could be done by combining antiangiogenic therapy with other therapies attacking the proliferating tumor, that is, chemotherapy or radiotherapy.

The model behavior shows the known fact that bevacizumab can not stop tumor growth. Due to the best of the author’s knowledge, this is the only tumor growth model incorporating the effect of angiogenic inhibition that possesses this property. This model can be extended later to incorporate the effect of other drugs applied in combined therapy with bevacizumab, which is one of the main motivations behind the research presented here.

Since the structure of the extended model is different from the original minimal model presented in [13], the comparison of the parameters resulting after identification can not be done directly, only after some transformations discussed below. Note moreover that identification was done in both cases by minimizing the least squared error from the measurements; in [13] only the measurements from mice C1–C4 were used; however, here we have used measurements from mice C1–C5 and E1–E9; thus the resulting parameters are necessarily different.

First consider the tumor dynamics without drug. In the minimal model [13], the tumor growth dynamics without drug is described by (consider (1) with ), while for the extended model, the tumor growth without drug is governed by ((4) with ). The tumor growth rate for the minimal model is 0.27 1/day, while the transformed “net” tumor growth rate for the extended model is .

Next, consider the effect of drug on the tumor growth dynamics in the time instant when the first big dose injection (10 mg/kg for the minimal model and the equivalent 0.171 mg/ml serum level for the extended model) is given. Denote the inhibition rate for the minimal model as , while denote the inhibition rate for the extended model as . Note that the dimension of is kg/(mg·day), while the dimension of is . When the injection appears, the effect of the drug on the tumor growth dynamics is given by mg/kg for the minimal model and for the extended model (since ). These effects are equivalent if we equate these terms; thus . This implies that the new inhibition rate should be compared to 10 times the old inhibition rate, that is, . The value of the new inhibition rate is 1/day, while mg/kg = 0.074 1/day.

The values of the “net” tumor growth rate () and the inhibition rate for the extended model are higher than the values for the minimal model in [13]. In the identification process for the extended model, measurements from mice E1–E9 were also used, while for the minimal model, only measurements from mice C1–C4 were used. Since for mice E1–E9 the inhibitor input and the resulting tumor volume are lower than for mice C1–C4; the inhibiting effect of the drug is higher according to the extended model.

The values acquired after the identification describe the average of the measurements and so the average of the population, not the tumor growth dynamics in a specific mice. However, with small perturbations in tumor growth rate, necrotic rate, and inhibition rate, all the individual mice measurements can be described with the proposed model. In our future works, we will create a mixed-effect model [41]; that is, we will determine the variance of the parameters based on the measurements.

6. 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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant Agreement no. 679681).

References

B. Döme, S. Paku, B. Somlai, and J. Tímár, “Vascularization of cutaneous melanoma involves vessel co-option and has clinical significance,” Journal of Pathology, vol. 197, pp. 355–362, 2002.
View at: Google Scholar
A. L. Harris, “Angiogenesis as a new target for cancer control,” European Journal of Cancer Supplements, vol. 1, no. 2, pp. 1–12, 2003.
View at: Publisher Site | Google Scholar
I. Ilic, S. Jankovic, and M. Ilic, “Bevacizumab combined with chemotherapy improves survival for patients with metastatic colorectal cancer: Evidence from meta analysis,” PLoS ONE, vol. 11, no. 8, Article ID e0161912, 2016.
View at: Publisher Site | Google Scholar
J. H. Distler, A. Hirth, M. Kurowska-Stolarska, R. E. Gay, S. Gay, and O. Distler, “Angiogenic and angiostatic factors in the molecular control of angiogenesis,” The Quarterly Journal of Nuclear Medicine and Molecular Imaging, vol. 47, no. 3, pp. 149–161, 2003.
View at: Google Scholar
J. Klamka, H. Maurer, and A. Swierniak, “Local controllability and optimal control for a model of combined anticancer therapy with control delays,” Mathematical Biosciences and Engineering, vol. 14, no. 1, pp. 195–216, 2017.
View at: Publisher Site | Google Scholar
F. S. Lobato, V. S. Machado, and V. Steffen, “Determination of an optimal control strategy for drug administration in tumor treatment using multi-objective optimization differential evolution,” Computer Methods and Programs in Biomedicine, vol. 131, pp. 51–61, 2016.
View at: Publisher Site | Google Scholar
N. Nath, T. Burg, D. M. Dawson, and E. Iyasere, “Optimizing antiangiogenic therapy for tumor minimization,” in Proceedings of the 2010 American Control Conference, (ACC '10), pp. 1242–1247, IEEE, Baltimore, MD, USA, July 2010.
View at: Google Scholar
U. Ledzewicz and H. Schättler, “Analysis of optimal controls for a mathematical model of tumour anti-angiogenesis,” Optimal Control Applications and Methods, vol. 29, no. 1, pp. 41–57, 2008.
View at: Publisher Site | Google Scholar
J. Sápi, D. A. Drexler, I. Harmati, Z. Sápi, and L. Kovács, “Qualitative analysis of tumor growth model under antiangiogenic therapychoosing the effective operating point and design parameters for controller design,” Optimal Control Applications and Methods, vol. 37, no. 5, pp. 848–866, 2016.
View at: Google Scholar
J. Sápi, D. A. Drexler, and L. Kovács, “Potential benefits of discrete-time controller-based treatments over protocol-based cancer therapies,” Acta Polytechnica Hungarica, vol. 14, no. 1, pp. 11–23, 2017.
View at: Google Scholar
P. Hahnfeldt, D. Panigrahy, J. Folkman, and L. Hlatky, “Tumor development under angiogenic signaling: a dynamical theory of tumor growth, treatment response, and postvascular dormancy,” Cancer Research, vol. 59, no. 19, pp. 4770–4775, 1999.
View at: Google Scholar
A. d'Onofrio and A. Gandolfi, “Tumour eradication by antiangiogenic therapy: analysis and extensions of the model by Hahnfeldt et al. (1999),” Mathematical Biosciences, vol. 191, no. 2, pp. 159–184, 2004.
View at: Publisher Site | Google Scholar | MathSciNet
D. A. Drexler, J. Sápi, and L. Kovács, “A minimal model of tumor growth with angiogenic inhibition using bevacizumab,” in Proceedings of the 2017 IEEE 15th International Symposium on Applied Machine Intelligence and Informatics, pp. 185–190, IEEE, Herl'any, Slovakia, 2017.
View at: Google Scholar
D. A. Drexler, J. Sápi, and L. Kovács, “Positive control of a minimal model of tumor growth with bevacizumab treatment,” in Proceedings of the 12th IEEE Conference on Industrial Electronics and Applications, pp. 2081–2084, IEEE, Siem Reap, Cambodia, 2017.
View at: Google Scholar
X. Wang and Y. Lin, “Tumor necrosis factor and cancer, buddies or foes?” Acta Pharmacologica Sinica, vol. 29, no. 11, pp. 1275–1288, 2008.
View at: Publisher Site | Google Scholar
S. Y. Proskuryakov and V. L. Gabai, “Mechanisms of tumor cell necrosis,” Current Pharmaceutical Design, vol. 16, no. 1, pp. 56–68, 2010.
View at: Publisher Site | Google Scholar
J. Sápi, L. Kovács, D. A. Drexler, P. Kocsis, D. Gajári, and Z. Sápi, “Tumor volume estimation and quasi- continuous administration for most effective bevacizumab therapy,” PLoS ONE, vol. 10, no. 11, Article ID e0142190, pp. 1–20, 2015.
View at: Publisher Site | Google Scholar
P. Érdi and J. Tóth, Mathematical models of chemical reactions. theory and applications of deterministic and stochastic models, Princeton University Press, NJ, USA, 1989.
D. A. Drexler and J. Tóth, “Global controllability of chemical reactions,” Journal of Mathematical Chemistry, vol. 54, pp. 1327–1350, 2016.
View at: Google Scholar
D. A. Drexler, E. Virágh, and J. Tóth, “Controllability and reachability of reactions with temperature and inflow control,” Fuel, vol. 211, pp. 906–911, 2018, http://www.sciencedirect.com/science/article/pii/S0016236117312036.
View at: Publisher Site | Google Scholar
F. Wu, M. Tamhane, and M. E. Morris, “Pharmacokinetics, lymph node uptake, and mechanistic PK model of near-infrared dye-labeled bevacizumab after IV and SC administration in mice,” The AAPS Journal, vol. 14, no. 2, pp. 252–261, 2012.
View at: Publisher Site | Google Scholar
H.-P. Gerber and N. Ferrara, “Pharmacology and pharmacodynamics of bevacizumab as monotherapy or in combination with cytotoxic therapy in preclinical studies,” Cancer Research, vol. 65, no. 3, pp. 671–680, 2005.
View at: Google Scholar
W. M. Haddad, V. S. Chellaboina, and Q. Hui, Nonnegative and Compartmental Dynamical Systems, Princeton University Press, NJ, USA, 2010.
View at: Publisher Site | MathSciNet
L. Goldshaid, E. Rubinstein, A. Brandis et al., “Novel design principles enable specific targeting of imaging and therapeutic agents to necrotic domains in breast tumors,” Breast Cancer Research, vol. 12, no. 3, article no. R29, 2010.
View at: Publisher Site | Google Scholar
“Genentech, Inc., ‘Prescribing information of Avastin’ (bevacizumab)”,” http://www.gene.com/download/pdf/avastin_prescribing.pdf.
View at: Google Scholar
J.-F. Lu, R. Bruno, S. Eppler, W. Novotny, B. Lum, and J. Gaudreault, “Clinical pharmacokinetics of bevacizumab in patients with solid tumors,” Cancer Chemotherapy and Pharmacology, vol. 62, no. 5, pp. 779–786, 2008.
View at: Publisher Site | Google Scholar
F. Kazazi-Hyseni, J. H. Beijnen, and J. H. M. Schellens, “Bevacizumab,” The Oncologist, vol. 15, no. 8, pp. 819–825, 2010.
View at: Publisher Site | Google Scholar
J. Gaudreault, G. Lieberman, F. Kabbinavar, and., and etal., “Pharmacokinetics (PK) of bevacizumab (BV) in colorectal cancer (CRC),” Clinical Pharmacology & Therapeutics, vol. 69, 25 pages, 2001.
View at: Google Scholar
M. H. Cohen, J. Gootenberg, P. Keegan, and R. Pazdur, “FDA drug approval summary: Bevacizumab (Avastin®) plus carboplatin and paclitaxel as first-line treatment of advanced/metastatic recurrent nonsquamous non-small cell lung cancer,” The Oncologist, vol. 12, no. 6, pp. 713–718, 2007.
View at: Publisher Site | Google Scholar
E. J. Israel, D. F. Wilsker, K. C. Hayes, D. Schoenfeld, and N. E. Simister, “Increased clearance of IgG in mice that lack β2-microglobulin: possible protective role of FcRn,” The Journal of Immunology, vol. 89, no. 4, pp. 573–578, 1996.
View at: Publisher Site | Google Scholar
V. Ghetie and E. S. Ward, “Multiple roles for the major histocompatibility complex class I-related receptor FcRn,” Annual Review of Immunology, vol. 18, pp. 739–766, 2000.
View at: Publisher Site | Google Scholar
“European Medicines Agency, “Scientific discussion of Avastin”,” http://www.ema.europa.eu/docs/en_GB/document_library/EPAR_-_Scientific_Discussion/human/000582/WC500029262.pdf.
View at: Google Scholar
L. Yan, K. Hsu, and R. A. Beckman, “Antibody-based therapy for solid tumors,” Cancer Journal, vol. 14, no. 3, pp. 178–183, 2008.
View at: Publisher Site | Google Scholar
T. H. Cheung and T. A. Rando, “Molecular regulation of stem cell quiescence,” Nature Reviews Molecular Cell Biology, vol. 14, no. 6, pp. 329–340, 2013.
View at: Publisher Site | Google Scholar
W. Chen, J. Dong, J. Haiech, M.-C. Kilhoffer, and M. Zeniou, “Cancer stem cell quiescence and plasticity as major challenges in cancer therapy,” Stem Cells International, vol. 2016, Article ID 1740936, 16 pages, 2016.
View at: Publisher Site | Google Scholar
Y. Cai, J. Zhang, and Z. Li, “Multi-scale mathematical modelling of tumour growth and microenvironments in anti-angiogenic therapy,” Biomedical Engineering Online, vol. 15, article no. 155, 2016.
View at: Publisher Site | Google Scholar
J. Sápi, D. A. Drexler, I. Harmati et al., “Tumor growth model identification and analysis in case of C38 colon adenocarcinoma and B16 melanoma,” in Proceedings of the 8th IEEE International Symposium on Applied Computational Intelligence and Informatics, pp. 303–308, IEEE, Timisoara, Romania, 2013.
View at: Google Scholar
J. Sápi, D. A. Drexler, Z. Sápi, and L. Kovács, “Identification of C38 colon adenocarcinoma growth under bevacizumab therapy and without therapy,” in Proceedings of the 15th IEEE International Symposium on Computational Intelligence and Informatics, pp. 443–448, IEEE, Budapest, Hungary, 2014.
View at: Google Scholar
L. M. Ellis and D. J. Hicklin, “VEGF-targeted therapy: mechanisms of anti-tumour activity,” Nature Reviews Cancer, vol. 8, no. 8, pp. 579–591, 2008.
View at: Publisher Site | Google Scholar
F. Fan, J. S. Wey, M. F. McCarty et al., “Expression and function of vascular endothelial growth factor receptor-1 on human colorectal cancer cells,” Oncogene, vol. 24, no. 16, pp. 2647–2653, 2005.
View at: Publisher Site | Google Scholar
T. Ferenci, J. Sápi, and L. Kovács, “Modelling tumor growth under angiogenesis inhibition with mixed-effects models,” Acta Polytechnica Hungarica, vol. 14, no. 1, pp. 221–234, 2017.
View at: Google Scholar
J.-B. Tylcz, K. El Alaoui-Lasmaili, E.-H. Djermoune, N. Thomas, B. Faivre, and T. Bastogne, “Data-driven modeling and characterization of anti-angiogenic molecule effects on tumoral vascular density,” Biomedical Signal Processing and Control, vol. 20, pp. 52–60, 2015.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2017 Dániel András Drexler et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

4737

Downloads

1437

Citations