Multiscale modeling of tumor response to vascular endothelial growth factor (VEGF) inhibitor
Abstract
Vascular endothelial growth factor (VEGF) has been known as a key mediator of angiogenesis in cancer. Bevacizumab is antiVEGF monoclonal antibody that has been approved by the FDA as a firstline treatment in many types of cancer. In this paper, we extend a previously validated multiscale tumor model to comprehensively include the multiple roles of VEGF during the course of angiogenesis and its binding mechanism with bevacizumab. We use the model to simulate tumor system response under various bevacizumab concentrations, both in standalone treatment and in combination with chemotherapy. Our simulation indicates that periodic administration of bevacizumab with lower concentration can achieve greater efficacy than a single treatment with higher concentration. The simulation of the combined therapy also shows that the continuous administration of bevacizumab during the maintenance phase can lead to antitumor activity which further suppresses its growth. Agreement with experimental results indicates the potential of the model in predicting the efficacy of antiVEGF therapies and could therefore contribute to developing prospective clinical trials.
1Introduction
Angiogenesis, the formation of new blood vessels from an existing vasculature, is a fundamental process that enables a solid tumor to sustain its growth and metastasize to distant organs. The driving force behind tumorinduced angiogenesis is the vascular endothelial growth factor (VEGF), that is released by starving tumor cells when the supply of oxygen or nutrient becomes limited. Inhibition of VEGF to suppress angiogenesis has become one major focus in cancer therapy in recent decades.
Experiments have provided vast evidences on multiple roles of VEGF during the course of angiogenesis. Upon released by the tumors, VEGF diffuses into the surrounding tissue and reaches the endothelial cells (EC) that form the lining of the blood vessel. The binding of VEGF with the receptor VEGFR located on the EC membrane triggers the sprouting of the new vessels. The tip of these sprouts is made of specialized ECs, called the tip cells. The ECs behind the tip cells are the stalk cells that proliferate and elongate the sprout, forming the trunk of the capillaries. The tip cells migrate and lead the growth direction of the sprout towards the tumor [1–5]. This migratory movement of EC tip cells is directed primarily by positive chemotaxis up the VEGF gradient and also by haptotaxis that occurs along the positive gradients of cellular adhesion sites in the extracellular matrix [6–8].
As in other living cells, ECs also undergo apoptosis, a programmed cell death that occurs naturally as part of its growth and development. In normal condition, apoptosis occurs as a balance between pro and antiapoptosis proteins. This balance can be disrupted by VEGFVEGFR binding, which triggers the activation of several other multiple downstream signaling pathways, in particular the PI3K/AKT pathway that regulates cell survival [9, 10]. Experiment by Nör [11] has shown the correlation between higher VEGF concentration and the upregulation of antiapoptosis protein BCl2. This result indicates the ability of VEGF to enhance EC survival by inhibiting cellular apoptosis, and thus promotes angiogenesis.
Due to multiple roles of VEGF in facilitating angiogenesis, various cancer therapies targetting VEGF have been developed in preclinical and clinical settings. Bevacizumab, also known by its brand name Avastin, is a cancer drug designed to prevent VEGF from binding to VEGFR located on the EC membrane. This blockage ultimately prevents the blood vessels from reaching the tumor. Bevacizumab was the first antiangiogenic drug approved by FDA for treatment of certain types of cancer, such as colorectal, kidney and lung cancers. It is administered intravenously and can be used alone or in combination with other methods, such as chemotherapy.
Understanding the complex dynamics of VEGF and its inhibitor is important in developing an effective treatment strategy. Although many preclinical and clinical trials have been performed extensively, assessment of the drug efficacy is still difficult to conduct, and even more so for combined therapies. Computational models are useful in providing a framework to study the detailed mechanism of VEGF inhibitor, and to investigate tumor response to different treatment strategies.
For a long time, mathematical models have been developed to study various antiangiogenesis therapies [12–20]. Finley et al. [21] provides a summary of developments in this area. Most models studied one or more aspects of antiangiogenesis strategies using either discrete or continuous approach. Continuous models consider the interactions between cell density and drug concentration that influence the viability of tumor and ECs. Although computationally inexpensive, these models can only capture largescale features of angiogenesis, without maintaining cellspecific properties and individual cell interactions. On the other hand, discrete models are able to capture singlecell variability and finer processes involved in angiogenesis at the expense of higher computational cost. Hybrid discretecontinuous models aim to maximize the advantages of both implementations within the same simulation. Many hybrid models for angiogenesis are also multiscale models that encompass several biological scales from the intracellular signaling up to vasculature formation at the tissue level.
In this study, we extend our previously tested multiscale hybrid discretecontinuous model [22, 23] to include detailed mechanisms of the role of VEGF in angiogenesis and individual EC response to VEGF and its inhibitor. We specifically aim to model: 1) VEGFVEGFR binding and EC tip activation, 2) EC apoptosis signaling and cell survival, and 3) tumor system response to bevacizumab and chemotherapy drug. The model extension covers several biological scales. At the molecular scale, we employ the pharmacodynamics Emax model to quantify the binding effect of VEGF to VEGFR on EC membrane, which serves as a trigger for angiogenesis initiation. We further combine the VEGF signaling with the apoptosis signaling pathway to determine the viability state of each EC. This molecular level model is connected to the tissue level through reactiondiffusion partial differential equations that govern the changes in VEGF, bevacizumab, and chemotherapy drug concentrations at the extracellular scale. The algorithm at the cellular scale is also extended to incorporate tumor response to chemotherapy.
We perform model simulation to analyze the efficacy of bevacizumab in suppressing angiogenesis and tumor growth. As a model system, we chose A549 lung adenocarcinoma. In year 2020, lung cancer is the leading cause of cancerrelated deaths in the world. A549 is a nonsmall cell lung cancer (NSCLC) cell line that has been used extensively as a testing ground for therapeutic drugs development, such as bevacizumab, paclitaxel and docetaxel, both in vitro and in vivo [15, 24, 25]. We compare the effect of different bevacizumab concentration and treatment strategies on the VEGF expression in the tumor tissue, on the EC apoptosis and the resulting capillary structure, and on the overall tumor volume. In addition, using the same model, we explore a combined treatment strategy with chemoterapy drug paclitaxel. Our simulation shows that continuous suppression of angiogenesis by bevacizumab can effectively reduce the overall growth of the tumor in the long run.
2Methods
We extend our previous model [22, 23] at several biological scales. In the extracellular scale, additional equations governing the changes in VEGF, bevacizumab and paclitaxel concentrations are added into the system. At the molecular scale, the apoptosis signaling pathway of the EC is now added to control its survival. At the cellular scale, additional feature is added into each tumor cell agent, allowing it to respond to chemotherapy drug, and at the tissue scale, each key process of vasculature formation is designed to be VEGFdependent, in accordance with experimental observations. Figure 1 illustrates how these biological scales are interconnected.
Fig. 1
We choose the hybrid discrete–continuous modeling approach to study the complex interplay between the tumor cells, endothelial cells, and their microenvironmental factors. The simulation developed in this work consists of two types of interacting environments: (a) the discrete cell fields representing tumor cells and ECs, and (b) the continuous fields representing the microenvironmental factors. Components in the continuous fields are modeled through a system of six coupled partial differential equations (PDEs) reflecting the spatiotemporal dynamics and interconnections between the following components: oxygen, VEGF, matrix degradative enzyme (MDE) and fibronectin, along with VEGF inhibitor (bevacizumab) and chemotherapy drug (paclitaxel). The system of coupled PDEs is solved using the finite difference scheme, in which both spatial domain and time interval are discretized [26]. The tumor cells and ECs are modeled as discrete agents residing in and confined to a fixed twodimensional domain.
2.1Microenvironmental factors and drugs dynamics
The extracellular scale is a central part of the model that connects the other components. In this model, we consider six extracellular microenvironmental factors, namely oxygen, VEGF, matrix degradative enzyme (MDE) and fibronectin, along with VEGF inhibitor (bevacizumab) and chemotherapy drug (paclitaxel).
Oxygen is a critical factor that determines tumor viability. The availability of oxygen facilitates tumor growth, while hypoxic tumor cells secrete VEGF that diffuses into the surrounding tissue, initiating EC capillary formation. ECs adhere to the extracellular matrix (ECM), which mainly consists of fibronectin that bounds to the surrounding tissue, and require ECM for their growth and migration during angiogenesis. Viable tumor cells produce MDE that diffuses and degrades the ECM, enabling the tumor cells to migrate within the tissue. In clinical treatment, both bevacizumab and paclitaxel are injected intravenously, delivered by the capillaries, diffuse into the tumor tissue and decay naturally. Bevacizumab binds with VEGF, while paclitaxel kills tumor cells directly.
The reactiondiffusion equations governing the spatial and temporal changes in the oxygen concentration u (x, t), MDE m (x, t), fibronectin f (x, t), VEGF v (x, t), bevacizumab B (x, t), and paclitaxel P (x, t) are as follows:
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
Paclitaxel, like most chemotherapy drugs, targets tumor cells directly and we will discuss the effect of paclitaxel on the proliferation of tumor cells in the following section. The parameters used in the extracellular scale and the values used in the simulation are listed in Table 1 in the Appendix. We use the parameter values that come from experiments as much as possible and only estimate those whose data are not available.
The coupled system of reactiondiffusion equations (1)(6) are defined across a twodimensional rectangular domain representing a tumor tissue. The domain is divided into grid of points with grid sizes Δx and Δy. During each timestep of length Δt, these chemicals diffuse across the grid points and cells consume from the nearby nodes. The central difference scheme is employed to discretize the diffusion term in the equations. The gradient at a grid point (i, j) is calculated by taking the average value of the gradient across the four orthogonal neighboring grid points. We calculate the changes in the chemical concentrations in a fixed time step Δt and write t = nΔt, for n = 0, 1, 2, . . .. As an example, the diffusion term D_{v} ∇ ^{2}v of VEGF in equation (4) at time t = nΔt can be approximated by
(9)
The other equations in the system (1)(6) can be implemented in a similar fashion.
2.2Tumor model
In the cellular scale, the agentbased approach is taken to model tumor phenotypic switch and migration. The objectoriented implementation strongly supports agentbased modeling paradigm. We construct a Tumor class, and each tumor cell is an object with a set of attributes, namely its position, viability status, age, and threshold value T_{c} that determines the cell’s resistance to chemotherapy drug. At each time step, a tumor cell executes a sequence of operations following the flowchart in Fig. 2.
Fig. 2
The operation performed by each tumor cell is primarily dependent on its current state. We determine whether a tumor cell is viable, quiescent, or necrotic based on its local oxygen concentration u (x, t). Two threshold values T_{q} and T_{n} are set and a cell is viable if u (x, t) > T_{q}, necrotic if u (x, t) < T_{n}, and quiescent otherwise.
While the quiescent cells release VEGF, the viable cells continue to proliferation check to determine if either cell division or cell migration should take place. The proliferation check consists of checking whether a certain time interval T has elapsed since the last cell division and whether there is available space in its surrounding for a cell to divide and give rise to a daughter cell. Physically, this waittime T ensures that the cell has completed all stages in its cell cycle. In our simulation, we set T to be the cell’s doubling time.
Paclitaxel, which is also known under the brand name Taxol, targets cells that are undergoing mitosis. It specifically prevents cells to achieve metaphase stage and prolong the activation of mitotic checkpoint. Furthermore, it stabilizes microtubules that leads to either cell cycle arrest or apoptotic cell death [27]. In our simulation, we set an individual threshold value T_{c} for each cell so that it can proceed with mitosis if the local paclitaxel concentration P (x, t) ≤ T_{c}. If P (x, t) > T_{c}, the cell pauses until the next iteration, and then randomly chooses whether it undergoes cell cycle arrest and becomes quiescent, or becomes apoptotic. We can capture cell’s variability in drug resistance by setting a random threshold value T_{c} for individual cell.
Viable cells that are not undergoing mitosis may migrate to nearby available space according to the discrete equation
(10)
(11)
(12)
2.3Angiogenesis model
The tissue scale of the model consists of the algorithm that governs the activation and dynamics of ECs, both the tip and stalk cells. We incorporate the apoptosis signaling pathway, which controls EC survival and apoptosis, in response to the change in VEGF concentration. This approach connects the tissue scale to the molecular scale via the extracellular scale.
Similar to the tumor model described in the previous section, we employ the objectoriented paradigm to implement the angiogenesis model. Each EC is an object with the following attributes: cell type, activation status, position, age, and apoptosis protein concentrations. Cell type is an attribute that takes on a value 1 if the EC is a tip cell and 0 if it is a stalk cell. The activation status has binary value as well. It is 1 if EC has been activated through the binding of VEGF to its receptor VEGFR, and 0 otherwise. As each EC object is equipped with an apoptosis signaling pathway, individual apoptosis protein concentrations are also included as part of the EC’s attributes. This enables us to capture cell variability and resistance to apoptosis. The VEGFdriven processes during the course of angiogenesis are described below.
VEGFVEGFR binding The binding of VEGF to its receptor VEGFR located on the tip of EC has been shown to be the key process during the initial stage of angiogenesis. We quantify the binding effect using the nonlinear Emax model that is often used in doseresponse analysis of drugs [13, 28, 29]. The effect E induced by VEGF upon binding to VEGFR is given by
(13)
Signaling pathway model for EC survival and apoptosis In addition to its role in EC tip activation, the upregulation of VEGF concentration has been shown to inhibit EC apoptosis and promote its survival.
Cellular apoptosis is controlled by two major pathways: the intrinsic and extrinsic pathways. We propose a simplified pathway that still captures the activation and inhibition of the key proteins without increasing the computational complexity as much. The simplified intrinsic and extrinsic pathways are shown in red boxes in Fig. 3. The intrinsic pathway is triggered by intracellular apoptosis stimulus, such as DNA damage, while the extrinsic pathway is triggered by death receptor TNF [30–33]. In response to DNA damage, proapoptosis proteins, such as Bax, are activated. The activation of Bax triggers the release of cytochrome c, which in turn activates Caspase 9 protein. In the extrinsic pathway, the tumor necrosis factor (TNF) activates Caspase 8 protein. Both intrinsic and extrinsic pathways converge to the activation of Caspase 3 protein leading to EC apoptosis.
Fig. 3
On the other end, the P13K/Akt pathway, repesented by blue boxes in Fig. 3, is activated once diffused VEGF reaches and activates the EC tip. VEGF upregulation induces the expression of BCl2, which is the antiapoptosis protein that suppresses the activation of Bax [11, 34]. In addition, BCl2 activates the proangiogenic chemokine CXCL8 that induces EC survival, proliferation and migration [20].
We incorporate this signaling pathway of the EC into the model to determine the fate of EC in response to VEGF upregulation or its inhibition by bevacizumab. The biochemical kinetics involved in the model are given in Table 4 in the Appendix. We use the law of mass action to translate the biochemical kinetics into a system of ODEs (Table 5) which govern the changes in the apoptosis protein concentrations. The reaction rate constants and initial value of each apoptosis protein are listed in Table 6.
Tip cells migration Upon activation by VEGF, the EC tip cells start to migrate while the stalk cells right behind them proliferate and sprout, forming new vasculature. We assume that the growth direction of the capillary sprout is determined by migratory direction of the EC tip cell. The EC stalk cells behind it proliferate and elongate following the tip cell, forming a contiguous sprout wall.
In our simulation setup, the tumor tissue is represented by a twodimensional square lattice, where we employ the agentbased approach with biased random walk to model the movement of EC tip on the lattice points. The EC stalk cell right behind the tip cell grows along the horizontal and vertical lines connecting the lattice points that have been previously occupied by the tip cell, forming a continuous capillary, as shown in Fig. 4(a).
Fig. 4
At any given time step, an EC tip cell occupies a grid point (i, j), indicated as node N_{5} in Fig. 4(b). The four grid points N_{1}, N_{2}, N_{3}, N_{4} orthogonally surrounding the center grid point N_{5} form the von Neumann neighborhood for this EC tip cell and these are the possible locations that the cell can migrate to in the subsequent time step. The EC tip cell can move up to position N_{1} (i, j + 1), move down to N_{2} (i, j  1), move right to N_{3} (i + 1, j), move left to N_{4} (i  1, j), or stay at its current position N_{5} (i, j).
For each k = 1, . . . , 5, we construct the probability function P_{k} that corresponds to each lattice point N_{k} to determine the new position of the EC tip in subsequent time step. We take into account the two biasing factors in EC migratory direction, namely the chemotaxis to VEGF and haptotaxis to the extracellular matrix (fibronectin), as proposed in [8, 16]. Thus,
(14)
(15)
(16)
Capillary sprouting and branching A common phenomenon during tumorinduced angiogenesis that has been experimentally observed is the brush border effect, in which the frequency of branching increases at the edge of the network as capillary sprouts get closer to the tumor [35]. Since spatial distance to the tumor correlates with the level of VEGF concentration, we employ positional information approach to model the branching probability as a function of local VEGF concentration. In [23], the branching probability
(17)
At any time step, a random number r ∈ (0, 1) is generated for each agent (the EC tip cell). If r ≤ P_{b}, then the EC tip will branch out and generates two new EC tips. Two random numbers s_{1}, s_{2} ∈ (0, 1) are generated and the two new EC tips will occupy gridpoint N_{k} for which s_{k} ∈ I_{k}, k = 1, 2. However, if r > P_{b}, then EC tip will migrate to an adjacent location, as discussed previously, followed by stalk cell elongation and sprouting.
Our algorithm would also allow two capillary sprouts to encounter each other and merge (anastomosis). In the simulation this process would occur if an EC tip happens to occupy a grid point that is currently or previously occupied by another EC tip. In this case, the algorithm decides randomly which sprout will continue to grow and deactivate the other one. The list of parameter values used in the tissue scale is available in Table 3 in the Appendix.
3Results and discussion
Tumor angiogenesis can be quantified through various parameters, such as the estimation of VEGF expression, MVD (Mean Vascular Density), Chalky counting, MAGS (Microscopic Angiogenesis Grading System) scoring, and Color Doppler study [36]. We implement the proposed model in two simulation setups. In the first set of simulations, we aim to analyze the effect of bevacizumab on VEGF expression, capillary density and structure, as well as EC survival. For the simulation domain, we consider a unit square D = [0, 1] × [0, 1] with parent blood vessel with 6 EC tip cells located along the left boundary (x = 0). The line tumor as the source of VEGF is placed along the right boundary (x = 1). The initial VEGF concentration is taken to be
(18)
(19)
In the second set of simulations, we analyze the effect of bevacizumab on the overall tumor growth, both when administered alone and when in combination with chemotherapy drug paclitaxel. Tumor algorithm in Section 2.2 is now incorporated into the simulation. We initially place 400 active and quiescent cells in the center of the simulation domain D = [0, 1] × [0, 1], forming a lump of cells with radius 0.1. The parent blood vessel is located along the four perpendicular sides of the domain, that is, along x = 0, x = 1, y = 0, and y = 1, with two randomly placed EC tip cells per side.
As the tumor used in the laboratory experiment has been inoculated and grown to a certain size, we assume that certain level of VEGF has already been present in the tumor tissue when the first dose of bevacizumab is administered. In our simulation, we take the initial VEGF distribution to be
(20)
We implement our model in C++ programming environment and use MATLAB for data analysis and visualization. For model validation, our simulation results are compared with the experimental data on A549 lung adenocarcinoma from Zhang et al [24].
3.1Effect of bevacizumab on VEGF expression
Bevacizumab works by directly binding the circulating VEGF, thereby lowering the expression of VEGF in the tumor tissue. Using the first simulation setup described earlier, we compute the total VEGF concentration in the tumor tissue ∑_{(x,y)∈D}v (x, y, t) under bevacizumab treatment as a percentage relative to the total VEGF concentration of the control group (without bevacizumab treatment).
Figure 5 shows the comparison of total VEGF expressions between the control and the case groups with concentrations of 3.0 and 6.0 mg/kg administered in two different ways. In the first treatment method, bevacizumab is administered once at the beginning of experiment, while in the second case, administration is done periodically every other day during the first ten days.
Fig. 5
We compare our simulation result with the experimental data from [24] on day 7, where bevacizumab is administered every other day. In Fig. 5, the total VEGF expression from the experimental control group on day 7 is taken to be the reference value. By comparing (a) and (b) in Fig. 5, we see that periodic treatment, even with lower concentration, shows higher efficacy in suppressing the VEGF expression. With the periodic treatment, bevacizumab with concentration 3.0 mg/kg can reduce the total VEGF expression by more than 45% at the end of day 10, while onetime treatment with concentration 6.0 mg/kg can only decrease the total VEGF expression by 35%.
3.2Effect of bevacizumab on capillary morphology
Using the same setup, we simulate the VEGFdriven capillary formation to investigate the effect of bevacizumab on the capillary structure and branching frequency.
Since the capillary branching probability is defined as a function of VEGF concentration in equation (17), our model is able to quantitatively measure the effect of different bevacizumab concentration on the branching frequency. Figure 6 shows that during the first four days, there is no significant difference between the control group and the case groups. However, at the end of day 10, bevacizumab concentrations of 3.0 and 6.0 mg/kg give approximately 50% and 65% reductions in the total number of branching when given once on day 1. If given periodically, the reduction is much greater, that is, 66% and 83% for concentration of 3.0 and 6.0 mg/kg, respectively.
Fig. 6
Counting branching frequency directly is difficult to achieve in a laboratory experiment. A more common approach to quantify angiogenesis is by calculating the vascular density, which is defined to be the total length of perfused vessels per unit area of observation field. Zhang et al [24] provides experimental data on the vascular density of the control group and case groups after 7 days of periodic bevacizumab treatment. As branching frequency correlates to vascular density and both can be used to measure tumor vascularization, we take these experimental data to test the model. The comparison between experimental data and simulation results on the branching frequency is shown in Figure 6(b), where the total number of branching of the control group on day 7 is taken to be the reference value.
Figure 7 shows the difference in capillary structure between the control and case groups. Branching and brush border formation are clearly seen in the control group, but these phenomena are less prominent as bevacizumab concentration increases.
Fig. 7
3.3Effect of bevacizumab on EC survival
To analyze the effect of various bevacizumab concentrations on EC survival and apoptosis, we implement the signaling pathway in Table 5. Since apoptosis is determined by the level of activated Caspase 3 protein and the activated CXCL8 level is directly associated with EC proliferation and migration, in our model we determine EC viability by taking the numerical difference between the level of these two protein concentrations.
Figure 8 shows the level of activated CXCL8 and Caspase 3 proteins with various bevacizumab concentrations under two different treatment methods. In (a), the level of activated CXCL8 associated with cell survival is constantly higher than the level of Caspase 3 that drives cellular apoptosis, although the difference decreases as bevacizumab concentration increases. In (b), bevacizumab is administered periodically every other day and the level of Caspase 3 protein starts to exceed the level of activated CXCL8 when the bevacizumab concentration is slightly above 3. This result suggests that periodic administration of bevacizumab is more effective even with a much lower concentration.
Fig. 8
3.4Effect of bevacizumab on tumor growth
We now combine the tumor and angiogenesis models described in Sections 2.2 and 2.3 to simulate the effect of bevacizumab on the overall tumor growth.
Recall that in the setup of the second set of simulations described earlier, the parent blood vessels are located along the boundary of the simulation domain D = [0, 1] × [0, 1]. Thus, we define the bevacizumab concentration at the administration time t_{k} by
(21)
For each treatment dosage, we run 10 simulations and record the average tumor radius, which is taken to be the largest distance from a tumor cell to the center of the lump at (0.5, 0.5). Assuming the spherical shape of the tumor mass, the volume can be estimated according to the formula V = (4/3) πr^{3}. Figure 9 shows the growth of the tumor during the first 10 days of treatment. The volume of the tumor in the control group on day 9 is taken to be the reference volume, where the volume at different days or different treatment groups are calculated as a percentage relative to this reference value. Compared with the control group, tumor growth in the case groups slows down significantly once treatment starts. Lower bevacizumab concentration (3 mg/kg) reduces the tumor volume on day 10 by 30% while bevacizumab with twice higher concentration can reduce the tumor volume up to approximately 55%, suggesting that bevacizumab suppresses tumor growth in concentration dependent manner.
Fig. 9
3.5Exploring bevacizumabpaclitaxel combined treatment
Experiments have shown promising efficacy data in combined therapy that involves chemotherapy drug and antiVEGF. Using the proposed model, we investigate the role of bevacizumab in a combined treatment with chemotherapy drug paclitaxel. Our simulation result indicates the importance of continuous bevacizumab administration during the maintenance phase after chemotherapy.
Incorporating equations (6) and (8) into the extracellular component of the model, we adopt the treatment schedule used in Ishikura et al [37] in which 20 mg/kg paclitaxel and 5mg/kg bevacizumab were administered on inoculated A549 lung cancer cells. Four treatment groups are considered. The first group is the untreated control group. In group 2, paclitaxel is administered on days 1, 8 and 15. In group 3, both paclitaxel and bevacizumab are administered simultaneously on days 1, 8 and 15. Finally in group 4, after both paclitaxel and bevacizumab are administered on days 1, 8, 15, weekly bevacizumab administration continues on during the maintenance phase.
Since both paclitaxel and bevacizumab are injected into the tissue intravenously, we set their concentration on the parent blood vessels located along the domain boundary at the administration time t_{k} to be
(22)
Figure 10 shows the comparison in tumor volume between these treatment groups. The volume of each group is calculated as a percentage relative to the volume of the control group on day 32. While paclitaxel alone can reduce the tumor volume by around 40%, combination with bevacizumab supresses the volume by approximately 55%. We note that in both groups, the tumor can still grow after the last treatment on day 15 before reaching their stable volume at around day 25. In group 4, continuous suppression of angiogenesis by bevacizumab may contribute to higher antitumor activity and smaller tumor volume.
Fig. 10
Figure 11 shows the cross section of the tumor morphology under these treatments on day 32. The proliferative cells are colored in green, while the quiescent and apoptotic cells are shown in blue and red, respectively. Under chemotherapy standalone treatment, proliferative cells still make up a good portion of the lump even though the apoptotic cells are visible around the perimeter. These viable cells can potentially keep growing unless paclitaxel is administered with sufficiently high concentration to eventually stop its growth and achieve stable volume. In groups 3 and 4, quiescent cells make up the biggest portion of the lump as bevacizumab suppresses angiogenesis and oxygen delivery. Although apoptotic cells are found on the perimeter in both groups, more viable cells are seen in group 3 than in group 4, and consequently the tumor in group 3 saturates at a higher volume. In group 4, the majority of the cells have become quiescent. These cells, although not apoptotic, cannot grow or divide.
Fig. 11
4Conclusion
Our previously developed multiscale tumor model has been extended to incorporate detailed mechanism of VEGF during angiogenesis and simulate tumor response to VEGF inhibitor, such as bevacizumab. The strengths of the proposed model come from the use of multiscale, hybrid discretecontinuous approach that provides transparency and flexibility in modeling different parts of the tumor system.
We adopt the Emax model to control the activation of EC upon binding with VEGF. A set of biochemical kinetics equations for the EC intracellular signaling is now incorporated into the molecular scale of the model and it takes into account the VEGF expression as a factor that influences cellular apoptosis and survival. By employing a discrete biased random walk in modeling VEGFdriven capillary formation, we can directly investigate the effect of bevacizumab on the capillary density and structure. Furthermore, it allows us to model the blood flow through the fine capillary network that delivers both the oxygen and antitumor drugs, rather than make a simplified assumption that any cell within a fixed distance from the blood vessel receives sufficient supply of these chemicals. Likewise, by employing the agentbased model for the tumor cell, we can include individual variability in cell viability and resistance to chemotherapy drug, mimicking the actual tumor system.
The inclusion of the EC apoptosis pathway allows us to quantify EC cell survival as the numerical difference between two main proteins that are responsible for cell apoptosis and survival. Although this would be difficult, if not impossible, to test in the experiment, the simulation result is able to show that periodic treatment with lower concentration can achieve greater efficacy than a onetime treatment with higher concentration. Moreoever, our simulation of the combined therapy also shows significant reduction in the tumor volume when continuous administration of bevacizumab is maintained post chemotherapy.
We realize that the proposed model still has shortcomings that can be addressed in future work. The current model considers VEGF as the only switch that can trigger angiogenesis even though there are many angiogenic pathways that may allow angiogenesis to occur even though VEGF pathway has been completely blocked. Additionally, although the current model is able to produce results that agree with the experiments, the simulation is run in twodimensional setting to minimize computational cost. Given the transparency and flexibility of the hybrid agentbased model, one could naturally extend this to three dimensions to achieve better precisions.
Author contributions
M.H. developed computational model, wrote computer codes, performed simulations, analyzed data and simulation results. J. S. formulated biological questions and performed literature study. Both authors designed the project and wrote the paper.
References
[1]  Gerhardt H. , Golding M. , Fruttinger M. , et al., VEGF guidesangiogenic sprouting utilizing endothelial tip cell filopodia, Journal of Cell Biology 161 (2003), 1163–1177. 
[2]  Yoshida A. , AnandApte B. and Zetter B.R. , Differential endothelialmigration and proliferation to basic fibroblast growth factor andvascular endothelial growth factor, Growth Factors 13(1996), 57–64. 
[3]  Paweletz N. and Knierim M. , Tumorrelated angiogenesis, Critical Review in Oncology/Hematology 9 (1989), 197–242. 
[4]  Terranova V.P. , Diorio R. , Lyall R.M. , et al., Human endothelialcells are chemotactic to endothelial cell growth factor and heparin, Journal of Cell Biology 101 (1985), 2330–2334. 
[5]  Suchting S. and Eichmann A. , Jagged gives endothelial tip cells anedge, Cell 137 (2009), 988–990. 
[6]  Cao Y. , Chen H. , Zhou L. , et al., Heterodimers of placenta growthfactor/vascular endothelial growth factor: endothelial activity,tumor cell expression, and high affinity binding to FLK1/KDR, Journal of Biological Chemistry 271 (1996), 3154–3162. 
[7]  StetlerStevenson W.G. , Matrix metalloproteinases in angiogenesis: amoving target for therapeutic intervention, Journal of ClinicalInvestigation 103 (1999), 1237–1241. 
[8]  Stokes C.L. and Lauffenburger D.A. , Analysis of the roles ofmicrovessel endothelial cell random motility and chemotaxis inangiogenesis, Journal of Theoretical Biology 152 (1991), 377–403. 
[9]  Lee S. , Chen T. , Barber C. , et al., Autocrine VEGF signaling isrequired for vascular homeostatis, Cell 130 (2007), 691–703. 
[10]  Zachary I. , VEGF signaling: integration and multitasking inendothelial cell biology, Biochem Soc Trans 31 (2003), 1171–1177. 
[11]  Nör J.E. , Christensen J. , Mooney D.J. and Polverini P.J. , Vascular endothelial growth factor (VEGF)mediated angiogenesis isassociated with enhanced endothelial cell survival and induction ofBcl2 expression, Am J Pathol 154 (1999), 375–384. 
[12]  Drexler D.A. , Sapi J. and Kovacs L. , Modeling of Tumor GrowthIncorporating the Effects of Necrosis and the Effect of Bevacizumab, Complexity 2017, Article ID 5985031, https://doi.org/10.1155/2017/5985031. 
[13]  Billy F. , et al., A pharmacologically based multiscale mathematicalmodel of angiogenesis and its use in investigating the efficacy of anew cancer treatment strategy, Journal of Theoretical Biology 260 (2009), 545–562. 
[14]  Gevertz J.L. , Computational Modeling of Tumor Response toVascularTargeting Therapies –Part I: Validation,Article ID 15, Computational and Mathematical Methods in Medicine 2011, Article ID 830515, doi: 10.1155/2011/830515. 
[15]  Wang J. , et al., Multiscale agentbased modeling on melanoma andits related angiogenesis analysis, Theor Biol Med Model 10 (2013), 41. 
[16]  Sun X. , Zhang L. and Tan H. , et al., Multiscale agentbased braincancer modeling and prediction of TKI treatment response: Incorporating EGFR signaling pathway and angiogenesis, BMC Bioinformatics 13 (2012), 218. 
[17]  Lai X. and Friedman A. , Mathematical modeling in scheduling cancertreatment with combination of VEGF inhibitor and chemotherapy drugs, J Theor Biol 462 (2019), 490–498. 
[18]  Mollard S. , et al., Model driven optimization of antiangiogenics +cytotoxics combination: application to breast cancer mice treatedwith bevacizumab + paclitaxel doublet leads to reduced tumor growthand fewer metastasis, Oncotarget 8 (2017), 23087–23098. 
[19]  Jain H.V. , Nör J.E. and Jackson T.L. , Quantification ofendothelial celltargeted antiBCl2 therapy and its suppression oftumor growth and vascularization, Mol Cancer Ther 8(2009), 2926–2936. 
[20]  Jain H.V. , Nör J.E. and Jackson T.L. , Modeling theVEGFBcl2CXCL8 pathway in intratumoral angiogenesis, Bull Math Bio 70 (2008), 89–117. 
[21]  Finley S.D. , Chu L. and Popel A.S. , Computational systems biologyapproaches to antiangiogenic cancer therapeutics, Drug DiscovToday 20 (2015), 187–197. 
[22]  Hendrata M. and Sudiono J. , A computational model for investigatingtumor apoptosis induced by mesenchymal stem cellderived secretome, Computational and Mathematical Methods in Medicine 2016, Article ID 4910603, https://doi.org/10.1155/2016/4910603. 
[23]  Hendrata M. and Sudiono J. , A hybrid multiscale model forinvestigating tumor angiogenesis and its response to cellbasedtherapy, In Silico Biology 13 (2019), 1–20. 
[24]  Zhang N. , Zhang G. , Zheng Y. , Wang T. and Wang H. , Effect of Avastinon the number and structure of tumor blood vessels of nude mice withA549 lung adenocarcinoma, Experimental and TherapeuticMedicine 8 (2014), 1723–1726. 
[25]  Franklin M. , A549  A Model for NonSmall Cell Lung Cancer, MiBioresearch, Retrieved January 20, 2021. 
[26]  Strikwerda J.C. , Finite Difference Schemes and Partial Differential Equations, SIAM, 2004. 
[27]  Bacus S. , et al., Taxolinduced apoptosis depends on MAP kinasepathways (ERK and p38) and is independent of p53, Oncogene 20 (2001), 147–155. 
[28]  Macdougall J. , Analysis of Dose–Response Studies–Emax Model. Dose Finding in Drug Development, Springer (2006), 127–145. 
[29]  Salahudeen M.S. and Nishtala P.S. , An overview of pharmacodynamicmodelling, ligandbinding approach and its application in clinicalpractice, Saudi Pharm J 25 (2017), 165–175. 
[30]  Fujio Y. and Walsh K. , Akt mediates cytoprotection of endothelialcells by vascular endothelial growth factor in ananchoragedependent manner, J Biol Chem 274 (1999), 16349–16354. 
[31]  Gerber H.P. , et al., Vascular endothelial growth factor regulatesendothelial cell survival through the phosphatidylinositol3kinase/Akt signal transduction pathway: requirement for Flk1/KDRactivation, J Biol Chem 273 (1998), 30336–30343. 
[32]  Gerber H.P. , Dixit V. and Ferrara N. , Vascular endothelial growthfactor induces expression of the antiapoptotic proteins Bcl2 and A1in vascular endothelial cells, J Biol Chem 273 (1998), 13313–13316. 
[33]  Spyridopoulos I. , et al., Vascular endothelial growth factorinhibits endothelial cell apoptosis induced by tumor necrosisfactoralpha: balance between growth and death signals, J MolCell Cardiol 29 (1997), 1321–1330. 
[34]  Pidgeon G.P. , Barr M.P. , Harmey J.H. and Foley D.A. , and D.J.BouchierHayes, Vascular endothelial growth factor (VEGF)upregulates BCl2 and inhibits apoptosis in human and murine mammaryadenocarcinoma cells, British Journal of Cancer 85(2001), 273–278. 
[35]  Muthukkaruppan V.R. , Kubai L. and Auerbach R. , Tumorinducedneovascularization in the mouse eye, J Natl Cancer Inst 69 (1982), 699–705. 
[36]  Sudiono J. and Thalib S. , Angiogenesis intensity within benign andmalignant oral mucosa epithelial tumor, Mol Cell Biomed Sci 2 (2018), 23–27. 
[37]  Ishikura N. , et al., Importance of bevacizumab maintenance followingcombination chemotherapy in human nonsmall cell lung cancerxenograft models, Anticancer Research 37 (2017), 623–630. 
[38]  Grote J. , Susskind R. and Vaupel P. , Oxygen diffusivity in tumortissue (DSCarcinosarcoma) under temperature conditions within therange of 2040C, Pugers Archiv 372 (1977), 37–42. doi: 10.1007/BF00582204. 
[39]  Casciari J.J. , Sotirchos S.V. and Sutherland R.M. , Mathematicalmodelling of microenvironment and growth in EMT6/Ro multicellulartumor spheroids, Cell Proliferation 25 (1992), 1–22. 
[40]  Hotary K. , Allen E.D. , Punturieri A. , Yana I. and Weiss S.J. , Regulation of cell invasion and morphogenesis in a 3dimensionaltype I collagen matrix by membranetype metalloproteinases 1, 2 and 3, Journal of Cell Biology 149 (2000), 1309–1323. 
[41]  Anderson A.R.A. , A hybrid mathematical model of solid tumourinvasion: the importance of cell adhesion, MathematicalMedicine and Biology 22 (2005), 163–186. 
[42]  Anderson A. and Chaplain M. , Continuous and Discrete MathematicalModels of Tumorinduced Angiogenesis, Bulletin of MathematicalBiology 60 (1998), 857–900. 
[43]  Papadopoulos N. , et al., Binding and neutralization of vascularendothelial growth factor (VEGF) and related ligands by VEGF Trap,ranibizumab and bevacizumab, Angiogenesis 15 (2012), 171–185. 
[44]  Welter M. and Rieger H. , Interstitial Fluid Flow and Drug Deliveryin Vascularized Tumors: A Computational Model, PLoS ONE 8(8) (2013), e70395. https://doi.org/10.1371/journal.pone.0070395. 
[45]  Cremasco M. and Wang L. , Estimation of partition, free and specificdiffusion coefficients of paclitaxel and taxanes in a fixed bed bymoment analysis: experimental, modeling and simulation studies, Acta Scientiarum 34 (2012), 33–40. 
[46]  A549 Cell Line: Cell Culture and Transfection Protocol, https://www.a549.com/cellsubcultureprotocol/, Retrieved January 20, 2021. 
Appendices
Appendix
Table 1
Symbol  Parameter  Value  Reference 
D_{u}  Oxygen diffusion coefficient  1.75 × 10^{5} cm^{2}/s  [38] 
γ_{u}  Tumor cell oxygen consumption rate  6.25 × 10^{17} M/cells/s  [39] 
∈  Degree of chemical localization  0.1  Est. in [22] 
δ_{u}  Oxygen decay rate  10^{}3/h  [13] 
p_{u}  Vessel permeability of oxygen  3 × 10^{4} cm/s  [16] 
r  Blood vessel radius  0.001 cm  [16] 
D_{m}  MDE diffusion coefficient  10^{9} cm^{2}/s  [40] 
μ_{m}  Single cell MDE production rate  1*  [41] 
δ_{m}  MDE decay rate  0*  [41] 
δ_{f}  ECM degradation rate  50*  [41] 
μ_{f}  Fibronectin production rate  0.05*  [42] 
D_{v}  VEGF diffusion coefficient  2.9 × 10^{7} cm^{2}/s  [42] 
μ_{v}  VEGF production rate  0.6 nM/h  [16] 
δ_{v}  VEGF decay rate  0.001/h  [14] 
γ_{v}  VEGF uptake rate  0.1 × 10^{4} /s  [42] 
β  Binding rate of VEGF and bevacizumab  35.1 pM  [43] 
D_{B}  Bevacizumab diffusion coefficient  4 × 10^{7} cm^{2}/s  [17] 
δ_{B}  Bevacizumab degradation rate  0.17 day^{1}  [17] 
p  Vessel permeability to drug in tumor tissue  1710^{3} μ/s  [44] 
D_{P}  Paclitaxel diffusion coefficient  4.27 × 10^{6} cm^{2}/s  [45] 
γ_{P}  Paclitaxel uptake rate  27.68 cm^{3}/g· day  [17] 
δ_{P}  Paclitaxel degradation rate  1.23 day^{1}  [17] 
*=nondimensionalized value
Table 2
Symbol  Parameter  Value  Reference 
C_{A}/C_{R}  Attraction to repulsion coefficient ratio  0.3*  [23] 
l_{A}, l_{R}  Attraction and repulsion lengthscales  0.5, 0.1*  [23] 
ξ  Haptotaxis coefficient  2600 cm^{2}/s/M  [41] 
T  A549 lung cancer cell doubling time  22 hours  [46] 
T_{q}  Threshold oxygen level for quiescent state  0.5  Est. in [22, 23] 
T_{n}  Threshold oxygen level for necrotic state  0.4  Est. in [22, 23] 
T_{c}  Threshold value of paclitaxel toxicity  0.39  Est. 
*=nondimensionalized value
Table 3
Symbol  Parameter  Value  Reference 
α  EC chemotaxis coefficient  2600 cm^{2} s^{1} M^{1}  [42] 
λ  EC haptotaxis coefficient  986 cm^{2} s^{1} M^{1}  [42] 
σ  Constant of VEGF chemotaxis sensitivity  1.667 × 10^{10*}  [42] 
E_{max}  Maximum effect of VEGF on EC activation  0.05 h^{1}  [13] 
K_{D}  Dissociation constant of VEGF at VEGF receptor  10^{7} M  [13] 
E_{active}  Threshold for EC activation  0 . 5^{*}  Est. 
*=nondimensionalized value
Table 4
VEGF + Bevacizumab

VEGF + Akt

Akt^{*} + BCl 2

BCl 2^{*} + CXCL 8

ApopStim + Bax

BCl 2^{*} + Bax*

Bax^{*} + Casp 9

Casp 9^{*} + Casp 3

TNF + Casp 8

Casp 8^{*} + Casp 3

Notation: A* = the activated state of protein A; A.B = the compound of proteins A and B; k^{+} = forward rate constant of reaction; k^{} = reverse rate constant of reaction.
Table 5
d [VEGF . Bevacizumab]/dt = β [VEGF] [Bevacizumab] = β [v] [B] 
d [Akt]/dt =  k_{1} [Akt] 






d [CXCL8^{*}]/dt = k_{5} [BCl2^{*} . CXCL8] 
d [Bax]/dt =  k_{6} [Bax] [ApopStim] 


d [ApopStim]/dt =  k_{6} [Bax] [ApopStim] 






d [Casp3^{*}]/dt = k_{11} [Casp9^{*} . Casp3] + k_{14} [Casp8^{*} . Casp3] 
d [TNF]/dt =  k_{12} [TNF] [Casp8] 
d [Casp8]/dt =  k_{12} [TNF] [Casp8] 


Notation: A* = the activated state of protein A; A.B = the compound of proteins A and B; k^{+} = forward rate constant of reaction; k^{} = reverse rate constant of reaction.
Table 6
Akt  [0, 1]  Akt^{*}  0  VEGF.Bevacizumab  0  FasL^{*}.Casp8  0 
BCl2  [0, 1]  BCl2^{*}  0  Akt^{*}.BCl2  0  Casp3  [0, 1] 
CXCL8  [0, 1]  CXCL8^{*}  0  BCl2^{*}.CXCL8  0  Casp3^{*}  0 
Bad  [0, 1]  Bad^{*}  0  BCl2^{*}.Bad^{*}  0  Bad^{*}.Casp9  0 
Casp9  [0, 1]  Casp9^{*}  0  Casp9^{*}.Casp3  0  ApopStim  0.1 
Casp8  [0,1]  Casp8^{*}  0  Casp8^{*}.Casp3  0  TNF  0.1 
All values are in nondimensional form. The value [0, 1] means a uniformly random number between 0 and 1. Notation: A* = the activated state of protein A; A.B = the compound of proteins A and B. Here the initial values of the proteins are set to be a random number between 0 and 1 representing the variability within the cells, while the activated form of these proteins are initially set to be 0.