Hybrid model of tumor growth, angiogenesis and immune response yields strategies to improve antiangiogenic therapy
Introduction
Blood vessels that supply tissues with nutrients and oxygen cannot meet the high demands of rapidly proliferating cancer cells during solid tumor progression, resulting in hypoxia and increased production of angiogenic factors. This increased demand for oxygen is met by the sprouting of new blood vessels from existing ones, a process known as tumor-induced angiogenesis1. Cancer cells under hypoxic conditions, overexpress tumor angiogenic factors, with the most common being the vascular endothelial growth factor (VEGF) that induces the formation of new vessels. These neo-vessels, however, have abnormal structures and are often over-leaky due to their rapid formation and immaturity, which renders them dysfunctional and not well-perfused1. This abnormality of the tumor vasculature contributes to hypo-perfusion, which in turn results in immunosuppression and limited transport of drugs and immune cells into the tumor1,2. Therefore, the tumor vasculature plays a crucial role in solid tumor initiation, growth, immune evasion, and response to therapies3. Normalization of the tumor vessels with anti-angiogenic drugs (e.g., anti-VEGF antibodies) is a clinically applied treatment to restore the structural abnormalities of the vessels and thus, improve tumor perfusion, oxygenation, and vascular function4,5,6. However, it is still unclear which patients will benefit from such treatments. To provide a mechanistic insight into the effects of anti-angiogenic treatment on tumor progression and immune cell behavior, we developed a mathematical model to investigate the interplay of tumor growth, angiogenesis, anti-angiogenic drugs, and the immune system.
Numerous models have contributed to our understanding of tumor microenvironment dynamic, particularly emphasizing the intricate interplay between vascular networks, interstitial fluid flow, and drug transport. Early work explored the impact of vascular compliance and leakiness on tumor blood flow patterns, revealing hindered drug delivery in central tumor regions due to non-uniform blood supply and interstitial fluid pressure7. Anderson and Chaplain8 employed a discrete model to simulate angiogenesis, considering endothelial cell motion influenced by angiogenic and matrix macromolecule factors and generating vascular networks resembling solid tumor vasculature. Subsequent studies coupled intravascular and interstitial flows, highlighting the influence of vascular leakiness on systemic flow patterns and tumor microenvironments, potentially contributing to cancer cell metastasis9. Additional models coupled tumor growth with angiogenesis, offering insights into the temporal and spatial dynamic of avascular and vascular growth stages10. Also, investigations combining in vivo experimental data with discrete vasculature simulations explored the impact of vascular normalization on nanomedicine delivery, revealing size-dependent effects of nanoparticles11,12,13,14. Furthermore, models employing discrete vascular remodeling over time emphasized the significance of mechanical forces in tumor-induced angiogenesis, highlighting the impact of the treatment timing and vessel wall permeability on the effective delivery of drugs15,16,17.
To this end, we built on previous research to develop a hybrid model that accounts for a continuous 3-dimensional tumor growth model, a discrete model of tumor-induced angiogenesis vessel remodeling, and a physiological-based kinetics model that simulates the immune cells transport and anti-tumor immune response. The current framework extends previous research by incorporating a detailed model of immune response and anti-tumor immunity that involves also the transport of immune cells and cytokines among the upper and lower trunk, the lymph and blood circulation, as well as the growth of a solid tumor (Fig. 1 and Supplementary information). The model also calculates the perfusion within the vascular network, which allows the quantification of functional vessels that directly affect the infiltration of the immune cells into the tumor. A comparison of the components incorporated in our modeling framework with previously published pertinent models7,8,9,10,11,13,15,16,17,18,19,20,21 is shown in Supplementary Table 1. To our knowledge, the current model is the first to include a non-regular vascular network, vascular flow through the network, interstitial flow, dynamic angiogenesis, tumor growth, and the immune system.
Schematic representation of the model compartments (a) and the main interactions of immune cells, cancer cells, cytokines, angiogenesis, anti-angiogenic treatment, and oxygen (b). a The compartments of the model and the immune cells and cytokines that are circulated through the entire mouse body. The mouse was divided into the following compartments: the tumor, upper trunk, lower trunk, lymph circulation, and blood circulation. The tumor is modeled as a 3D compartment with a spherical geometry, and the rest of the mouse (upper trunk, lower trunk, blood, and lymph circulation) are modeled as homogenous (0D) compartments. All the compartments are connected with blood and lymph circulation and the balance of all variables was kept through all compartments. Dendritic cells, antigen-presenting cells, CD8+ effector T-cells, and cytokines (immunosuppressive and pro-inflammatory) can enter the tumor via the tumor vessels and return back to circulation via the tumor-draining lymphatic vessels. The dendritic cells that are activated to present the tumor antigen can be transported to the lymph node and activate naïve T-cells against the tumor antigen. b The oxygen (O2) as a nutrient is supplied by the blood vessels and fuels the cancer cells proliferation. Hypoxia induces the VEGF production which attracts the nearby vessels to move through the VEGF gradient. Cancer cells induce immunosuppressive cytokines which can induce immune cell exhaustion which can be inhibited by pro-inflammatory cytokines. CD8+ effector T-cells can induce cytolysis in the cancer cells, and dendritic cells (DC) and antigen-presenting cells (APC) can phagocytosis the cancer cells. The antigen-presenting cells activate the naïve T-cells like CD8+ T-cells against tumor antigen (CD8+ effector T-cells). CD8+ effector T-cells and APC produce pro-inflammatory cytokines. Further details are provided in the Supplementary Information. Created with BioRender.com.
We first compare model predictions and provide mechanistic insights to published experimental data that involved anti-VEGF treatment in a murine glioma model. Then we present model predictions for the effect of the anti-VEGF treatment on vessels’ blood velocity, oxygen concentration, and tumor growth, and provide insights into how anti-VEGF treatment can improve blood vessel functionality. Finally, we present the spatiotemporal distribution of antigen-presenting cells, CD8+ effector T-cells, and pro-inflammatory and immunosuppressive cytokines, following anti-VEGF treatment.
Results
Model validation with experimental data of anti-VEGF therapy
First, we set out to compare model predictions with data from the response of mouse GSC005 glioblastoma to anti-VEGF treatment4. Figure 2a presents the experimental and simulated protocols. The simulation starts on day −10 of the experiment, a day after tumor implantation, as the model needs a small finite tumor volume to start with. After day −10 the simulation and the experimental protocol overlap, as shown in Fig. 2a. Anti-VEGF treatment started on day 0 (day 10 of simulation), and the last dose was given on day 9 (day 19 of simulation).
a The timeline of the experiment and the simulation (created with BioRender.com). b Sensitivity of the model by increasing the immune response (killing potential of immune cells) and the transport (extravasation) of immune cells. c Validation of the model, comparison of the simulated (dashed curves) with experimental tumor (continuous curves). The average values of the experimental tumor growth data are presented along with their standard errors. The simulated curves are presented with a range (shaded areas) by varying +−5% and +−2.5% (light-shaded and darker-shaded areas, respectively) the tumor growth rate constant (parameter lc) and cancer cells’ random mobility (i.e., their invasiveness to the surrounding tissue, parameter Dc.) from their baseline values.
Figure 2b shows the comparison of the experimental tumor growth curves with the model predictions for a baseline immune response and for increasing the model parameters related to the extravasation of immune cells into the tumor and their killing potential. Specifically, the blue dashed line represents the case of a baseline immune response without any vascular normalization. By further increasing the extravasation of the immune cells into the tumor compartment and the killing potential of immune cells (i.e., increased immune response) the tumor can be eliminated (orange dashed line). Given that the experimental tumor volume of the anti-VEGF group is lower than the control group and at the final stages of tumor growth both curves have a similar slope (Fig. 2b), we surmised that the lower tumor volume in the case of anti-VEGF was due to a reversible effect during treatment. The reversible effect of anti-angiogenic treatment could be a transient delay of angiogenesis which causes a reduction of oxygen and nutrient supply and not due to a significant anti-tumor immune response. A robust anti-tumor immune response would persist even after the last dose of anti-VEGF and tumor volume would follow a similar behavior to the orange curve in Fig. 2b, which is the effect of a strong anti-tumor immune response. In the pertinent literature, it is also well established that proper doses of anti-angiogenic treatment can induce vessel normalization, which is defined by the pruning of some immature/dysfunctional vessels and reduction in vessel permeability18,19. These changes in tumor vasculature have the potential to improve the efficacy of treatments, such as immunotherapy and/or chemotherapy. However, vascular normalization treatment alone often is not effective20, which explains the increase in tumor growth after the termination of anti-angiogenic treatment observed in the experimental data in Fig. 2b.
The analysis of the experimental data along with the predictions of the model and the established mechanisms of anti-angiogenic treatments18,19,21,22 leads us to the following hypotheses: (i) anti-VEGF treatment could mitigate and temporally delay the angiogenesis and (ii) normalize the vasculature by reducing the permeability of the vessels and/or pruning immature vessels. This in turn would affect overall tumor perfusion and oxygenation. Furthermore, in the experimental data of Fig. 2, vascular normalization with anti-VEGF treatment could affect immune cell infiltration but it cannot lead to a sufficient increase in the immune response in order to cause robust anti-tumor effects (i.e., the immune stimulation by anti-VEGF was not sufficient by itself to significantly reduce tumor volume). To prove our hypotheses and the underlying mechanisms, we set to validate the experimental data with our model by incorporating the necessary mechanisms and the appropriate parameters for these mechanisms. The simulated and experimental curves for each group are presented in Fig. 2c where all the parameters of the immune system, cell extravasation, and tumor growth are the same for both groups, and only the anti-angiogenic treatment changed according to the experiment (for the values of the model parameters see Supplementary Table 2). For the simulated growth curves a range of tumor volumes was calculated by varying by +−5% and +−2.5% the tumor growth rate constant and cancer cells’ random mobility (i.e., their invasiveness to the surrounding tissue) from their baseline values. These parameters were selected due to their clinical implication in glioblastoma and other cancers. Specifically, more aggressive tumors, such as high-grade glioblastoma, proliferate faster and are more invasive23,24,25, and thus, they exhibit a higher tumor growth rate constant and cancer cells’ random mobility. These types of cancer cells are also able to shift from a proliferative to an invasive phenotype based on the availability of nutrients/oxygen or space26,27. This shifting of phenotypes arises naturally in our model because the tumor growth rate constant and cancer cells’ random mobility are parameters that define the maximum proliferation and invasion rates, not the actual rates; the actual rates depend also on the availability of oxygen and space which determine if the invasion or proliferation will prevail.
An improvement in tumor response can be achieved by anti-VEGF monotherapy, but the effect is reversible according to both simulated and experimental tumor curves and it is observed only within the time window of anti-VEGF treatment. The increase in the proliferation and invasion parameters by +5% or +2.5% leads to greater tumor volumes of 50 mm3 and 40 mm3, respectively (see the upper range of blue-shaded areas in Fig. 2c). The use of anti-VEGF treatment for such cases results in a higher reduction of tumor volume by 30 mm3 to 40 mm3 compared to the less proliferative and invasive tumors (-5% or -2.5%), where the reduction in tumor volume ranges from 10 mm3 to 2 mm3. Overall, the trend in experimental and simulated curves is in good agreement, which is considered acceptable given the complexity of the biological system. Further improvement of the validation can be achieved by varying extensively the parameters that affect the angiogenesis, the tumor growth, and the immune system which is computationally demanding and will not change the main conclusions for the underlying mechanisms. Because of the large variation in the experimental data of tumor growth, even within the same treatment group, we chose to follow the trend of the data rather than overfit our model by extensively varying all model parameters. Overfitting the experimental data will not improve the predictability of our model, on the other hand, it will capture the random fluctuations in the data. To justify our approach, we also present the individual experimental data points along with the model predictions (Supplementary Fig. 1).
Effect of anti-VEGF treatment on blood vessels’ functionality
Next, we investigated the effect of anti-VEGF treatment on blood velocity and distribution of vascular density. The results presented in this section are model predictions from the simulations for control and anti-VEGF treatment (Fig. 2c). We assumed that anti-VEGF not only affects the tumor vasculature by blocking VEGF-induced angiogenesis but also reduces the permeability of the vessels (parameter Lp).
Figure 3a, b presents the distributions of the average cross-sectional velocity in the vessels with and without anti-VEGF treatment, and for high and low vessel permeability (Lp), respectively, whereas Fig. 4 depicts the spatial representation of the vascular density. The distributions of the average velocity are generated by exploiting the discrete vascular network and intravascular pressures at two different time points (at day 16 and the end of the simulation) and at different regions of the tumor compartment (i.e., peritumoral vs. intratumoral region). Interestingly, all the distributions are shifted to higher blood velocities when the permeability of the vessels is reduced due to anti-VEGF treatment. By comparing the velocity distribution in the intratumoral and peritumoral regions, we infer that it is mostly the intratumoral blood velocity that is affected by anti-VEGF treatment, and the effect is greater at the early stages of the angiogenesis (Fig. 3a, b).
The distribution of average blood velocity in the vessels with and without anti-VEGF (Lp = 6 10−9 and Lp = 6 10−10 [m/(mmHg s)], respectively). a on Day 16 and b at the end of the simulation. The distributions are shown for different regions of the tumor compartment: entire tumor, peritumoral and intratumoral regions.
The spatial distribution of the vascular density and functional vascular density of neo-vessels after treatment for control and anti-VEGF on Day 16 (a) and at the end of the simulation (b).
Given that the blood velocity in the vessels determines their functionality, we calculated the percentage of the functional vessels that have a sufficient average cross-sectional velocity greater than 0.1 mm/s28 for the distributions of Fig. 3a, b. The results for the percent of functional vessels are presented in Table 1 for the entire tumor as well as for the intratumoral and peritumoral regions, with and without anti-VEGF treatment at two time points: on Day 16 and at the end of the simulation. The results show that anti-VEGF increases the percentage of the functional vessels by about 10%, 3%, and 23% for the entire tumor, peritumoral, and intratumoral regions, respectively, on Day 16. These values, however, are reduced significantly until the end of the simulation to 4%, −3 and 10%. This temporal improvement in vessel functionality with anti-VEGF treatment can be explained by the vascular normalization window14,21,29.
In Fig. 4 the spatial distributions of the entire vessels and the functional vascular density are presented (a) on Day 16 and (b) at the end of the simulation with and without anti-VEGF treatment. Without anti-VEGF (Control), the vascular density always covers a greater volume compared to anti-VEGF treated tumors where the angiogenesis is delayed due to the anti-VEGF treatment. We note also that although on Day 16 the coverage of functional vascular density looks similar, a great proportion of vessels are dysfunctional in control as it is shown in Fig. 3 and Table 1. At the end of the simulation, the difference in the functional vascular density is more pronounced: for control, the functional vessels are localized and randomly spread, whereas for anti-VEGF treatment, the functional vessels are spread more homogeneously inside the tumor.
The effect of anti-VEGF on oxygen concentration and the amount of viable cancer cells in comparison to the control group on day 16 is shown in Fig. 5. The anti-VEGF treatment delays the angiogenesis by reducing the VEGF and this delay affects the oxygen supply which affects the oxygen concentration on the tumor periphery and reduces tumor growth rate. In control, the oxygen concentration is almost uniform, which makes the tumor growth rate independent from oxygen concentration. Finally, anti-VEGF-treated tumors have lower tumor volumes and cancer cells proliferate mostly at the center of the tumor until the termination of anti-VEGF treatment. However, cancer cells migrate in the radial direction inside healthy tissues during anti-VEGF treatment and if vascular normalization is not combined with proper anti-cancer therapy or a sufficient immune response, the cancer cells start to proliferate rapidly after the termination of the anti-angiogenic treatment. The uninhibited collective migration and/or regrowth of tumor cells on the periphery of the tumor can be explained by the “go or grow” mechanism26,27,30 and the shift of cancer cells from proliferative to invasive and again to proliferative phenotype. The “go or grow” mechanism arises naturally in the current model due to the dependence of cancer cells’ growth rate on oxygen concentration and available space (i.e., term (it {l}_{c}frac{{c}_{{ox}}}{{K}_{{c}_{{ox}}}+{c}_{{ox}}}{varPhi }_{c}{varPhi }_{f}) which depends on oxygen concentration ({c}_{{ox}}) and fluid volume fraction (it {varPhi }_{f}) that becomes zero when there is no available space for growth, for more details see the first reaction term of viable cancer cells in Supplementary Table 3).
The spatial distribution of oxygen concentration [mole/m3] and viable cancer cells (volume fraction) for control and anti-VEGF on day 16 of the simulation.
Spatiotemporal distribution of immune cells and cytokines
Subsequently, we set out to investigate the effect of improved perfusion caused by anti-VEGF treatment on the immune system. An increase in blood velocity (i.e., perfusion) can improve the transport of immune cells to the tumor site, and as more immune cells arrive at the tumor, their infiltration inside the tumor region might be increased. Figure 6 presents the spatial distribution of CD8+ effector T-cells, viable and dead cancer cells for control and anti-VEGF treatment at the end of the simulation. The CD8+ effector T-cells are more homogeneously distributed for the anti-VEGF treated tumors compared to the control group, presumably due to improved perfusion (Fig. 4). The anti-VEGF treatment has the highest volume fraction of CD8+ effector T-cells at the center of the tumor. On the contrary, untreated tumors have a high-volume fraction of CD8+ effector T-cells located randomly in the periphery similar to the distribution of the functional vascular density (Fig. 4). Also, the CD8+ effector T-cells induce cancer cell cytolysis, which results in a high accumulation of dead cancer cells and a reduction of viable cancer cells as shown in Fig. 6. For the control group, the dead cancer cells are mostly in small areas where the oxygen levels are too low (i.e., necrosis due to hypoxia). Overall, anti-VEGF treatment not only reduces tumor growth during treatment but also normalizes the vascular network, which results in a more intense and uniform extravasation of immune cells, an increased amount of dead cancer cells, and a reduction in viable cancer cells. Despite the increased anti-tumor immunity that normalization treatment can induce, in the experimental data considered in our study, the immune response was not sufficient to restrict tumor growth. It seems that the immune system needs even more time to reach sufficient anti-tumor effect and thus, the anti-angiogenic treatment should be combined with targeted immunotherapies31 and other non-targeted anti-cancer treatments to further boost anti-tumor immunity.
The spatial distribution of CD8+ effector T-cells, viable and dead cancer cells volume fraction for control and anti-VEGF at the end of the simulation.
The activation of the adaptive immune system is also shown by the time evolution of antigen-presenting cells (APCs) and CD8+ effector T-cells at the peripheral compartments as well as from their average value at the tumor compartment (Fig. 7a, b). Although initially, the APCs and CD8+ effector T-cells in anti-VEGF treated tumors increase at a lower rate than in the control group, there is a point where both cell types exceed the corresponding values of the control group (Fig. 7a, b). The initial lower rate of activation of APCs and CD8+ effector T-cells can be explained by the lower tumor volume and lower tumor antigen, which is caused by the delay of the angiogenesis after the anti-VEGF treatment. At longer times, the activation of APCs and CD8+ effector T-cells increases due to improved tumor perfusion caused by the normalization of the vasculature (anti-VEGF treatment), which in turn causes a higher and more uniform extravasation of immune cells inside the tumor compartment. Also, the amount of pro-inflammatory cytokines is increased with anti-VEGF treatment, whereas the amount of anti-inflammatory cytokines is decreased which could be explained by the immunosupportive effect of vascular normalization31.
The time evolution of antigen-presenting cells, CD8+ effector T-cells, pro- and anti-inflammatory cytokines: a in the peripheral compartments (upper and lower trunk, lymph, and blood circulations) and b the tumor compartment for the control versus anti-VEGF treatment.
Discussion
In this study, we formulated a comprehensive mathematical model capturing the dynamic of solid tumor growth, VEGF-driven angiogenesis, anti-angiogenic drugs, and pivotal cancer-immune-cell interactions. Our model was compared with previously published experimental data4 and served as a tool to examine critical mechanisms influencing the effectiveness of anti-angiogenic treatment in solid tumors and their effect on the local (tumor microenvironment) and systemic immune response. The following mechanisms arise from the mathematical model and the pertinent literature18,19,21,22 (i) anti-VEGF treatment could mitigate and temporally delay the angiogenesis, resulting in reduced tumor growth, (ii) normalization of the vasculature by reducing the permeability of the vessels and/or pruning immature vessels can increase the functional vascular network with a more uniform distribution inside the tumor, (iii) vascular normalization can improve immune cells infiltration and their uniform distribution within the tumor. However, prolonged anti-VEGF treatment can cause negative outcomes because during anti-angiogenic therapy, cancer cells transition towards an invasive phenotype and migrate inside the healthy tissue, excessive vessel pruning increases the levels of hypoxia with immunosuppressive effects, and at the end of the anti-angiogenic treatment, cancer cells can rapidly proliferate inside the healthy tissue by transitioning back to a proliferative phenotype. Overall, the positive contribution of the anti-angiogenic treatment must balance the possible negative effects by choosing a proper treatment protocol.
Interestingly, we found that anti-angiogenic treatment increased the number of functional intratumoral vessels by 23% and 10% on day 16 and at the end of the simulation, respectively. The computational results indicate that the percentage of functional vessels with anti-angiogenic treatment is increased within a period that defines the normalization window29. The normalization window can be defined as the time interval that optimizes tumor perfusion with the use of anti-angiogenic drugs. Our mathematical model can thus be employed to provide insights into the mechanisms of vascular normalization strategies combined with targeted immunotherapies and/or other non-targeted anti-cancer treatments given that vascular normalization alone cannot boost anti-tumor immunity, whereas targeted immunotherapy alone cannot be effective due to the heterogeneity of cancer cells’ antigens23,32.
While our developed hybrid model significantly advances the understanding of complex interactions within solid tumors and the basic immune responses, it is important to acknowledge certain limitations inherent to its design and assumptions. Firstly, the model’s complexity, incorporating a combination of continuous 3D tumor growth, discrete tumor-induced angiogenesis, and immune response compartments, introduces computational challenges that may limit the feasibility of extensive parametric analysis. The stochastic calculation of vessel formation in the discrete domain, while capturing the inherent randomness of angiogenesis, poses challenges in terms of computational efficiency. Furthermore, our model assumes a simplified representation of the tumor microenvironment and host tissue, overlooking potential heterogeneities that exist within a real tumor and its environment. Another possible effect of angiogenesis and anti-angiogenic treatment is the potential dynamic of acid production and variations in H+ concentrations in the tumor. The acidic environment in the tumor can induce immunosuppression33,34 and it should be investigated both with experiments and mathematical modeling. Mathematical modeling requires the incorporation of acid production and aerobic and anaerobic glycolysis as well as their effects on immune response, which would further increase the complexity of the current model. Additionally, the model assumes a uniform distribution of immune cells and cytokines throughout the four compartments that comprise the entire body, neglecting potential regional variations that may impact the local immune response of specific organs. These assumptions were made in order to reduce the unknown parameters for each organ/compartment. It must be noted that our model considers only one tumor antigen but in reality, the continuous mutations of tumor cells could result in more than one antigen which can affect the immune response. Despite these limitations, our model serves as a valuable tool for exploring key mechanisms influencing the effects of vascular normalization and immune response in solid tumors, offering insights that can guide further refinement and development of more sophisticated models and optimized treatment in the future.
Overall, we present here a methodology to include immune-system response to hybrid-type models composed of a discrete tumor vasculature and a continuous solid tumor growth. Our model aligns well with clinical observations, particularly in the context of anti-VEGF treatment for high-grade glioblastoma (relapsed or progressing glioblastoma). As highlighted in our study and supported by the literature35,36,37, anti-VEGF therapy has demonstrated promise in prolonging survival and improving the quality of life for patients with relapsed or progressing glioblastoma. Notably, the parametric analysis in Fig. 2c, focusing on crucial mechanisms, such as the tumor growth rate and cancer cell random mobility (i.e., their invasiveness to the surrounding tissue), indicates increased efficacy of anti-VEGF treatment, especially in high-grade gliomas, with a greater volume reduction compared to less proliferate and less invasive gliomas. Furthermore, the reduction in the interstitial fluid pressure (Supplementary Fig. 2) post-anti-VEGF treatment, as indicated by our simulation results, correlates with the observed clinical outcome of reduced edema in glioblastoma patients with anti-VEGF treatment36. While acknowledging the challenges of directly comparing simulated interstitial fluid pressure values with clinical data due to the absence of such data in the literature38; the observations of the clinical trial and the model are in agreement. By considering that anti-VEGF monotherapy did not induce sufficient immune response as found in clinical trials and predicted by our model, we propose exploring combination therapies of anti-VEGF treatment with targeted immunotherapy and non-targeted anti-cancer therapies, to directly target cancer cells with known and unknown tumor antigens. Our model could be helpful in studying the underlying mechanisms of such combination therapies and provide insights for better treatment outcomes.
Methods
The model was built in COMSOL Multiphysics v5.6 (Burlington, MA) and interacts with MATLAB codes for the calculation of angiogenesis and intravascular pressure. Box 1 presents an overview of the algorithm, the interactions between MATLAB and COMSOL Multiphysics as well as the core equations. The mathematical framework, the complete set of equations, and the solution strategy are provided in detail in the Supplementary Information.
Continuous tumor domain model
In the 3D continuous tumor compartment, the solid ((it {varPhi }_{s})) and fluid ((it {varPhi }_{f})) phases conservation is calculated39 as well as the volume fraction of viable cancer cells ((it {varPhi }_{c})), dead cancer cells ((it {varPhi }_{d})), CD8+ effector T-cells ((it {varPhi }_{{T}^{E}})), antigen-presenting cells ((it {varPhi }_{{APC}})) and dendritic cells ((it {varPhi }_{{DC}})). The conservation of all cells in the tumor compartment was calculated as solid phases by considering their accumulation, convection with solid velocity (({{bf{v}}}_{s})) in the conservative form to satisfy the volume change, diffusion, and reaction (the reaction terms are given in Supplementary Table 3); see Box 1 conservation of solid phases. The mechanical stress is also calculated by decomposition of the deformation gradient tensor (({bf{F}})) in elastic deformation (({{bf{F}}}_{e})) and inelastic deformation due to tumor growth (growth deformation gradient tensor, ({{bf{F}}}_{g}))12,40,41,42. The inelastic growth deformation gradient tensor was calculated by the growth stretch ratio (({lambda }_{g})). The growth stretch ratio depends on the total accumulation of the solid phase (i.e., viable minus dead cancer cells). The 2nd Piola-Kirchhoff stress tensor was calculated using the constitutive equation of the Neo-Hookean material and solving for the momentum balance; see equations in Box 1 and section “Solid mechanics—mechanical stress in tumor compartment” in Supplementary information for more details. The interstitial fluid pressure ((p)) is also calculated by Darcy’s equation and the conservation of medium (solid and fluid phases), which considers the fluid supply by blood vessels (pre-existing and neo-vessels) and the fluid absorption by the lymphatic vessels; see conservation of medium in Box 1 and section “Conservation of solid and fluid-phase in the tumor compartment” in Supplementary information for more details. The mass balance of oxygen (({c}_{{ox}})), cytokines pro-inflammatory (({c}_{p})) and immunosuppressive (({c}_{a})), and angiogenic factors VEGF (({c}_{{VEGF}})) are solved as diluted species; see conservation of diluted species in Box 1. The conservation of diluted species was calculated by considering their accumulation, convection with average bulk velocity (({{bf{v}}}_{{ab}})) in conservative form to satisfy the volume change, diffusion, and reaction; the reaction terms of these species are given in Supplementary Table 4. Matrix metalloproteinases (({c}_{{MMPs}})) and fibronectin a part of the extracellular matrix (({c}_{{EXM}})) are also calculated (their equations are given in the section “Mass balance of diluted species in the tumor compartment” of Supplementary information).
Discrete model of tumor-induced angiogenesis
The process of angiogenesis is initiated to satisfy the increased demand for oxygen by the cancer cells. The motion of discrete vessels is calculated stochastically in the 3D discrete model in MATLAB with probabilities P0 − P6. The tips of the vessels have random mobility, a motion due to the gradient of soluble angiogenic factors, such as VEGF (chemotaxis), and a motion in response to the gradient of bound molecules within the extracellular matrix (haptotaxis). All these mechanisms define the values of the probabilities P0 – P68,9,43. The probabilities and the motion of discrete vessels are calculated in the 3D discrete model whereas the VEGF and extracellular matrix values are calculated in the 3D continuous model (COMSOL Multiphysics) and transferred to the 3D discrete model (Box 1). For more information see the section “Discrete mathematical model of endothelial cells motion (angiogenesis) in tumor compartment” in Supplementary information. The vascular density (({Sv})) and functional vascular density ((S{v}_{f})) are transferred in the 3D continuous model. The vascular densities (({Sv}) and (S{v}_{f})) are calculated by adding the pre-existing vascular density and the total surface of discrete vessels inside the volume of a tissue element divided by this volume (see Box 1). Only discrete vessels with an average velocity higher than 0.1 mm/s contribute to the functional vascular density28,44. For the discrete neo-vessels, the volumetric fluid flow rate in each vessel and transvascular fluid flow rate into interstitial space are also calculated by considering the mass balance of fluid. All these quantities depend on both the vascular fluid pressure (({p}_{v})) and interstitial fluid pressure ((p)). The vascular pressure (({{p}_{v}}_{j})) was calculated for each j-th discrete vessel node as it is presented in Box 1. The fluid flow through vessels’ walls is calculated by Starling’s approximation by using the ({{p}_{v}}_{j}) and (p) which are calculated in the 3D discrete model and 3D continuous model, respectively12. The coupling of all these quantities was satisfied by transferring the total flow across the vessels’ wall in the 3D continuous model, which contributes to the calculation of the interstitial fluid pressure and then the fluid pressure is transferred back to the 3D discrete model, see the schematic representation in Box 1 and for more information see section “Calculation of volumetric fluid flow rate in vessels and transvascular fluid flow rate into interstitial space of the tumor compartment” in Supplementary information. The calculation of flow through vessels allows the calculation of vessels’ average cross-sectional velocity and subsequently the functional vascular density which is used to calculate the infiltration of the immune cells into the tumor compartment. It must be noted that the diluted species like oxygen supply depends on the vascular density (not the functional vascular density) because the oxygen is a small molecule and can be transported by diffusion without the need for sufficient perfusion like immune cells.
Modeling of transport of cells and cytokines between compartments
The body of the mouse is divided into compartments, and they are connected as presented in Fig. 1 and Supplementary Fig. 3. All the immune cells (CD8+ effector T-cells, antigen-presenting cells, and dendritic cells) and the pro-inflammatory and immunosuppressive cytokines are transported through the different compartments. The mass balance for all these variables in the 0D compartments is given by the following equation:
Where the first term is the accumulation of each quantity in each 0D compartment, the second term is the incoming quantity from other 0D compartments, the third term is the outcoming quantity to other 0D compartments, the fourth term describes the reactions that take place inside the 0D compartments and the final term connects the blood and lymph circulations (0D compartments) with the 3D tumor compartment by integrating the exchange of mass from 0D compartments to the 3D tumor compartment44. The incoming and outcoming volumetric flows and the volume of each compartment are given in Supplementary Table 5. For more details see the section “Immune cells and cytokines transport through pharmacokinetic/pharmacodynamic compartments” of the Supplementary information. Their reaction terms (({R}_{{y}_{c}})) as well as the mass transported in and out of the tumor compartment (({T}_{{y}_{c}})) are defined in Supplementary Table 6.
Responses