A multi-scale model for determining the effects of pathophysiology and metabolic disorders on tumor growth

The search for efficient chemotherapy drugs and other anti-cancer treatments would benefit from a deeper understanding of the tumor microenvironment (TME) and its role in tumor progression. Because in vivo experimental methods are unable to isolate or control individual factors of the TME and in vitro models often do not include all the contributing factors, some questions are best addressed with systems biology mathematical models. In this work, we present a new fully-coupled, agent-based, multi-scale mathematical model of tumor growth, angiogenesis and metabolism that includes important aspects of the TME spanning subcellular-, cellular- and tissue-level scales. The mathematical model is computationally implemented for a three-dimensional TME, and a double hybrid continuous-discrete (DHCD) method is applied to solve the governing equations. The model recapitulates the distinct morphological and metabolic stages of a solid tumor, starting with an avascular tumor and progressing through angiogenesis and vascularized tumor growth. To examine the robustness of the model, we simulated normal and abnormal blood conditions, including hyperglycemia/hypoglycemia, hyperoxemia/hypoxemia, and hypercarbia/hypocarbia – conditions common in cancer patients. The results demonstrate that tumor progression is accelerated by hyperoxemia, hyperglycemia and hypercarbia but inhibited by hypoxemia and hypoglycemia; hypocarbia had no appreciable effect. Because of the importance of interstitial fluid flow in tumor physiology, we also examined the effects of hypo- or hypertension, and the impact of decreased hydraulic conductivity common in desmoplastic tumors. The simulations show that chemotherapy-increased blood pressure, or reduction of interstitial hydraulic conductivity increase tumor growth rate and contribute to tumor malignancy.

and ECM production/degradation are all important processes within the TME that affect tumor biology 5,[27][28][29] , but are difficult to study using experimental systems. The biochemical signaling pathways are also directly influenced by physical abnormalities in the TME. For example, high interstitial fluid pressure (IFP) within tumors causes increased interstitial fluid flow (IFF) at the tumor boundary, which affects biochemical transport and distribution within the tumor 24 .
Once a primary solid tumor exceeds a threshold volume (~1 mm) 30 , its central cells become starved because the outer surrounding cells deplete the nutrients; in response, the central cells secrete a variety of morphogenic and chemotactic growth factors 31 such as vascular endothelial growth factor (VEGF), angiopoietin-1, and -2 (Ang-1 and Ang-2) 4,32 . VEGF binding to its receptor, VEGFR-2, on ECs induces differentiation of ECs and formation of tip endothelial cells (tECs) that extend from the existing vessel wall to create sprouts in the initial stage of angiogenesis 33,34 . tECs lead the migration of vessel sprouts into the surrounding tissue, biased toward areas of high VEGF concentration (chemotaxis), and traveling through regions where the ECM supports cell-matrix interactions (haptotaxis) 35,36 . During migration, tECs secrete MMPs, which proteolytically degrade structural components of the ECM 5,37 . Following behind the tECs are the stalk endothelial cells (sECs), which are induced by NOTCH signaling to proliferate and form lumens 5,33,34 . As the migrating vessels meet and form connections, this process of VEGF-induced angiogenesis creates a complex network of immature neo-vessels that supports blood flow into the tumor. Another cytokine, Ang-1 causes maturation and stabilization of neo-vessels through binding to Tie-2 of ECs and, in contrast, Ang-2 promotes death of ECs and destabilization of neo-vessels by competing with Ang-1 for Tie-2 binding sites [38][39][40] .
In this work, we developed a comprehensive, fully-coupled multi-scale mathematical model that recapitulates a three-dimensional TME to better understand the transition from avascular growth to angiogenesis and VTG. As shown in Fig. 1, our model includes subcellular-, cellular-, and tissue-level size scales. The subcellular scale is the basic scale of the model, and consists of the biochemical agents: ECM (fibronectin), MMP, VEGF and VEGFR-2, Ang-1 & Ang-2 and their common receptor Tie-2, and CR-agents. At the cellular scale, proliferating, quiescent, and necrotic phenotypes of TCs are determined through a novel CR-based method. Unlike many previous studies that consider oxygen as the only factor determining cell viability, we include multiple CR-agents -oxygen, glucose and carbon dioxide -to determine cellular activities and phenotypes. For the ECs involved in angiogenesis, we use agent-based modeling to determine proliferation, quiescence persistence, death and migration of the cells. At the tissue scale, solid tumor growth and blood vasculature are implemented by a new double hybrid continues-discrete (DHCD) model with the highest degree-of-freedom (DOF) for movement of TCs and tECs in three-dimensional space. At the tissue level, we also model vessel remodeling, which is determined by biochemical and mechanobiological signals from blood shear stress with accurate hemodynamics and hemorheology 57 . We also consider blood pressure explicitly, and calculate interstitial fluid flow, which influences the transport of soluble species.

Methods
The overall algorithm of the present multi-scale TME model is shown in Fig. 1. Subcellular, cellular and tissue scales and their components are represented by green, yellow, and blue colors, respectively. All model compartments and their communications were designed to facilitate calibration of model parameters based on experimental observations for all stages of the simulations (avascular tumor growth, angiogenesis and vascular tumor growth). Subcellular scale: biochemical agents. According to the CR process, glucose, oxygen and carbon-dioxide influence the production of mitochondrial-synthesized adenosine-triphosphate (ATP), which is the basis of cellular energy 28,58 . Here, we use a set of reaction-convection-diffusion mathematical equations to model absorption-distribution-metabolite-excretion (ADME) of each CR-agent. The ADME equations for concentrations of oxygen, c o 2 , glucose, c g , and carbon-dioxide, c co 2 , are respectively written in Eqs. 1-3:  D is the diffusion coefficient of each respective agent, u ins is the interstitial fluid flow (IFF) velocity, C and R are consumption and production rates, respectively, E is the perfusion rate, and L indicates the location of ECs (L EC ), and TCs (L TC ). IFF and convection of soluble agents is considered throughout the domain. This is especially relevant near the vessels, where the Peclet number can be as high as ~1. Glucose and oxygen diffuse and are advected in the tissue from the source blood vessels. Carbon dioxide is produced in the tissue and taken up by the blood vessels, which are represented as sinks in the model. Assumed to be a function of cellular vitality (ϑ), consumption of glucose and oxygen, and production of carbon-dioxide by TCs are related by the stoichiometry of the CR-reaction according to: Figure 1. The overall algorithm for simulating a three-dimensional TME with subcellular, cellular, and tissue size scales. γ 0 is the maximum consumption or production rate of the CR agent. The transport of oxygen and carbon-dioxide -two hydrophobic molecules -is calculated using a simplified perfusion model based on transvascular pressure, p (the difference between intravascular blood pressure, p lum , and interstitial fluid pressure, p ins ) and vessel diameter, d v 19 . However, transport of glucose, a hydrophilic polar molecule, is more complex 59 . In addition to transvascular pressure and vessel diameter, efflux from blood and influx to the TCs can affect the glucose transfer 60 . Mammalian TCs have transmembrane symporter proteins that automatically intake glucose without the need for ATP consumption 61 . Therefore, we include only a Michaelis-Menten (M-M) model to determine glucose efflux from the blood to TCs. Finally, the convection components of oxygen, glucose and carbon-dioxide transvascular transport are derived through Eq. 5a-c, respectively: is the rate constant for oxygen transport into tissue, f g and km g are the maximum rate and M-M constant of glucose transport from blood, f co 2 is the transport rate of carbon-dioxide, and d c is the characteristic diameter of a neo-vessel. TCs exposed to oxygen concentration below a threshold level, c o ch 2 , secrete VEGF to stimulate ECs of nearby vessels to sprout. Once ECs become associated with the tumor tissue, their VEGFR-2 receptors are expressed 38 . The ADME model of VEGF coupled with VEGFR-2 is given in Eqs. 6-8: R v is the rate of production of VEGF by TCs, f v is the VEGF transport rate between tissue and vessels, k v + is the binding rate of VEGF to VEGFR-2, k v − is the dissociation rate of VEGF from VEGFR-2, v ε is the natural excretion rate of VEGF, r v f is the concentration of free-VEGFR-2, and r v a is the concentration of active-VEGFR-2 bound to VEGF. As such, Ang-1 is secreted by ECs in the vicinity of both healthy and tumor tissue, and Ang-2 is secreted by ECs associated with tumor tissue [38][39][40] . The ADMEs of Ang-1, c a1 , and Ang-2, c a2 , coupled with Tie-2 are mathematically modeled through Eqs. 9-13. R a1 and R a2 are respectively the secretion rates of Ang-1 by ECs and of Ang-2 by ECs associated with tumor tissue. e 0 is the characteristic concentration of ECs in each blood vessel, K a is the carrying capacity coefficient of angiopoietins, + k a1 and − k a1 are respectively Ang-1 binding rate to and unbinding rate from Tie-2, k a2 + and k a2 − are respectively Ang-2 binding rate to and unbinding rate from Tie-2, a1 ε and a2 ε are the natural excretion rates of Ang-1 and Ang-2, respectively. r a f is the concentration of free Tie-2, r a a 1 and r a a 2 are concentrations of active Tie-2 bound to Ang-1 and Ang-2, respectively.
a a a a f a a a a 1 1 MMPs secreted by both ECs and TCs degrade ECM to make space available for TCs and ECs to spread 5,37 . Based on the proteolytic mechanism of MMPs, the coupled ADME equations of MMPs and ECM are presented in Eqs. 14 and 15, respectively, c m and c e are respectively concentrations of MMPs and ECM (fibronectin), R m T , and R m E , are respectively MMP secretion rates by TCs and ECs, ε m is natural excretion rate of MMPs and e ε is excretion rate of ECM by MMPs, respectively. We assume that MMP secretion by TCs (marked by L TC ) dynamically varies based on cellular vitality, ϑ. cellular scale: agent-based cellular phenotypes. TC phenotypes. The concept of cellular vitality, ϑ, is considered to determine TC activities and phenotypes. In the CR-based model presented in Eq. 16, the nutrients in the CR reaction, oxygen and glucose, increase ϑ and the waste product of the CR reaction, carbon-dioxide, decreases ϑ. In this equation, ϕ is a proportionality coefficient, c o ch 2 , c g ch , and c co ch 2 are oxygen, glucose, and carbon dioxide characteristic concentrations, respectively 62 2 is a Heaviside function to ensure that carbon dioxide reduces ϑ when it's concentration, c co 2 , exceeds the characteristic value, c co ch 2 . In this model, the TCs with ϑ below ϑ ch are assumed to be quiescent (unable to migrate and proliferate) and those with ϑ above ch ϑ are active (prone to migrate and proliferate) 19 . As a biological assumption, the quiescent TCs can be converted to both proliferative and necrotic phenotypes based on their ϑ and cellular energy (ψ) values, respectively, and the proliferative TCs can become quiescent based on ϑ; however, necrotic TCs cannot be converted to the other phenotypes.
In the cellular energy model presented in Eq. 17, the quiescent TCs consume ψ at a constant rate, k q c , but the active ones produce ψ at a linear rate related to ϑ with a proportional coefficient, k a p ; they also consume ψ based on a M-M model with maximum rate k a c and M-M's constant 1 19 . Quiescent TCs with negative ψ are converted to necrotic phenotype. The active TCs need a characteristic energy, ψ ch , before they can proliferate into two daughter TCs with half value of ch ψ 63,64 . Here, ψ introduces available units of ATP molecules as a function of mitochondrial activity 19,61 . Based on the CR-reaction, a mathematical model of coupled cellular vitality and cellular energy is presented in Eqs. 16 phenotypes. An agent-based model determines the proliferation, migration, quiescence persistence, and death of the tECs and sECs. The cell phenotypes affect neo-vessel growth dynamics during angiogenesis. The tECs migrate toward positive gradients of VEGF and ECM (fibronectin) 34,36,65 . The sECs migrate into the tECs-generated conduits in ECM and also proliferate and create lumens of the neo-vessels 5 www.nature.com/scientificreports www.nature.com/scientificreports/ sECs can differentiate into tECs in response to high VEGF concentration, and thus generate bifurcating branches from the neo-vessel wall 19,65 . The inactive quiescence and death states are also considered for sECs based on VEGF and Ang-1/Ang-2 concentrations (see Eq. 21). tissue scale: solid tumor growth and angiogenesis. Fibronectin is a provisional ECM component abundant in tumors that serves as an adhesive glycoprotein and modulates mechanical stiffness. It is important for cell adhesion and migration and supports haptotactic migration of adhesive cells, including tECs and active TCs 36 . tECs can also use chemotaxis to migrate in response to VEGF gradients 34,65,66 . Moreover, any live cell can randomly walk in the ECM due to intracellular actomyosin-cytosolic dynamics and extracellular stimuli [67][68][69][70] . Corresponding to these movement mechanisms of tECs and TCs, the tumor growth and angiogenesis are mathematically modeled by Eqs. 18 and 19, respectively. Implementation of the DHCD method to calculate the highest DOF movement probabilities for each cell (tumor and endothelial cells) derived from these equations is presented in Supplementary Material 71 . ρ and ρ TC are respectively EC and TC densities, D EC and D TC are respectively diffusivity of ECs and TCs in the interstitium, α is the saturation coefficient for chemotaxis, and β c and β h are weight coefficients for chemotaxis and haptotaxis, respectively. Accumulation of rapidly dividing TCs increases the mechanical, compressive solid stress, p s [72][73][74] For this additional tumor growth-induced stress, we implement a Gaussian-like function for accumulative systems 19 . Tumor growth-induced solid stresses can also influence tumor growth and angiogenesis. Therefore, we assume that TCs are pushed toward low pressure/stress regions and tumor-induced solid stress acts as a barrier for tEC migration 73,75,76 .
The direct effect of IFF is to drive mass transport of molecular agents. It is generally thought that the pore size of the interstitium is too small for large object such as cells to be carried with the IFF. However, because cell migration is influenced by molecular gradients, IFF indirectly affects cell migration by affecting chemotaxis. Although we and others have shown that IFF can be a signal for cell migration or differentiation 7,77,78 , this mechanism was not a focus of this study. Therefore, we assume that the interstitial flow doesn't directly affect the cellular dynamics and phenotype. tissue scale: vessel growth and remodeling. After new vessels form via angiogenesis, they need to form lumens before flow can proceed. They do this through a process of lumenogenesis, which is controlled by VEGF, Ang-1, and Ang-2. The equation for neo-vessel diameter, d v , is written as Eq. 20, α and θ p are, respectively, the maximum rate and the M-M constant for sEC proliferation; ω m , θ m , and m ρ are the maximum quiescence persistence rate -which is a parameter that controls how long the cell remains quiescent before maturing to form lumen (quiescence-for-maturation) -for sEC, the M-M constants for maturation relative to VEGF and Ang-1/Ang-2, respectively; δ and θ d are the maximum rate and M-M constant for sEC death, respectively. k age is also a positive constant that represents the neo-vessel lumen growth due to aging 19 .
Vessel adaptation. Angiogenic neo-vessels are highly sensitive to biomechanical and biochemical stimuli in the TME. Respectively shown in Eq. 22a-d, the TME stimuli caused by endothelial wall shear stress (WSS), S WSS , Scientific RepoRtS | (2020) 10:3025 | https://doi.org/10.1038/s41598-020-59658-0 www.nature.com/scientificreports www.nature.com/scientificreports/ transvascular pressure, S p , hematocrit-induced metabolite, S HM , and VEGF, S vegf , are the major signals that control neo-vessel dilation/constriction 21,32,80-82 . ( ) S k log log p log (100 86 exp( 5000 ( ( ( ))) )) (22b) p p 10 10 10 In Eq. 22a, τ is the WSS in a neo-vessel segment, which depends on blood dynamic viscosity, blood µ , lumen flow rate, Q lum , and vessel diameter, d v . It is calculated using Hagen-Poiseuille's law: . τ ref is a positive constant to avoid singularity at low WSS. In Eq. 22b, k p is a proportional coefficient and p is transvascular pressure. In Eq. 22c, k m is a proportional constant describing the intensity of metabolic stimuli, Q ref is the maximum flow rate within the neo-vessel network, determined by the inflow supplied by the primary vessels 83 ; H D is the hematocrit of blood, and Q lum is the blood flow rate within the lumen. In Eq. 22d, k e is a proportionality coefficient, k 0 the vessel wall elasticity constant per unit length in the absence of VEGF, c v ch is the characteristic VEGF concentration for which the vessel wall susceptibility to adaptation defined by d c , d 0 stress-free vessel diameter, and τ c is a characteristic elastic stress introduced to ensure calculation stability 21 .
Finally, the TME-induced adaptation of vessel diameter is considered in Eq. 23. In addition to four TME stimuli, a shrinkage term, S Sh , is included to allow the possibility for vessels and decrease their diameter if the other stimuli are not sufficient for growth 21,82 Vessel deformation. We assume that the neo-vessels are mechanically compliant with a constant elasticity, E and compliance power, cp, as well as a collapse pressure, p c 55,84,85 . Therefore, the transvascular pressure, p, and tumor growth-induced solid stresses, p s , through Eq. 24 can deform the neo-vessels and change the vessel diameter from Vessel disruption. In addition to vessel growth and remodeling mechanisms, we consider the possibility of vessel disruption. If a neo-vessel has insufficient blood flow and τ falls below a threshold value, th τ , during a survival time, t s , then the vessel is mathematically eliminated 20,42,86 . Similarly, vessels inside the tumor exposed to VEGF concentrations higher than a threshold value, c v th , during the survival time, t s , are also disrupted 21 .
Tissue scale: Hemodynamics and interstitial fluid flow in the TME. For the fluid dynamics model, we consider flow in three regions spanning from a neo-vessel's lumen to the interstitium: intravascular blood flow, transvascular fluid flow, and interstitial fluid flow (IFF). Because of the very small diameter of newly-generated vessels, the Womersley number of flow is very small, and hence the pulsatile effects of the cardiac cycle can be ignored 87,88 . Therefore, Hagen-Poiseuille's law as an exact solution of the Navier-Stokes fluid dynamics equation can be applied for intravascular blood flow. IFF is calculated using Darcy's law assuming low porosity of healthy and tumor tissues 18,22,52,87 . Finally, transvascular fluid flow is calculated by Starling's law, which couples IFF and intravascular blood flow 18,22,52 . The continuity equation for intravascular blood flow is written as Eq. 25. In this equation, Q lum is the blood flow rate in the lumen calculated as the difference between intravascular blood flow rate based on Hagen-Poiseuille's law, Q IBF , and transvascular fluid flow rate based on Starling's law, Q TFF , respectively, shown in Eqs. 26 and 27.

TFF v p l um ins
In Eq. 25, N = 26 is the number of peripheral vessel lattice nodes adjacent to the central vessel node, and β describes the direction of lumen blood flow; it is +1 for outlet flow from a peripheral node and −1 for inlet flow to a peripheral node. In Eq. 26, p lum is intravascular blood pressure, d v is the neo-vessel diameter, L is the length of a neo-vessel segment which can be equal to 1, 2 , or 3 times the lattice length up to the neo-vessel path line in three-dimensional space, and µ blood is the dynamic viscosity of non-Newtonian blood as a function of neo-vessel diameter, d v , and blood hematocrit, H D . In Eq. 27, L p is the hydraulic conductivity of the neo-vessel wall, which is a constant for healthy tissue and increases in tumor tissue based on vessel diameter 20,55 ; p is transvascular pressure, σ is the osmotic reflection coefficient, lum π and ins π are colloid osmotic (oncotic) pressures of the intravascular plasma and interstitial fluid, respectively. Oncotic pressure is a type of osmotic pressure that depends on a concentration difference across a semipermeable membrane structure. Oncotic pressure can only be transmitted by the fluid phase of the tissue, and not through solid elements such as matrix, cells, or the vessel wall 89 .
The continuity equation for IFF given in Eq. 28 shows both incompressibility of plasma fluid by zero divergence of interstitial fluid velocity, u ins , and leakiness of the neo-vessel wall due to a source term in vascular tissue. Darcy's law determining IFF is also shown in Eq. 29.

u L S V p At the vessel wall
In the extravascular space ins ins ins u ins is the velocity of IFF, K ins is the interstitial hydraulic conductivity (IHC) of the TME, S V / is the surface area of the neo-vessel per unit volume for mass transport in the interstitium 18,50 . By replacing Eq. 28 in Eq. 29, we can derive the Poisson-Laplace's equation for IFP shown in Eq. 30.

the vessel wall
In the extravascular space

(30)
ins p ins lum i ns 2 tissue scale: hemorheology. Because of the non-Newtonian behavior of blood, its dynamic viscosity as a function of vessel diameter, d v , and blood hematocrit, H D , is calculated using Eq. 31. Indeed, we assume that blood is a non-Newtonian suspension red blood cells within a Newtonian plasma fluid 57,82,83,90 . In Eq. 31, µ blood and µ plasma are the dynamic viscosity of blood and plasma, respectively; n µ is the relative dynamic viscosity of healthy blood with a normal hematocrit, H D n , defined in Eq. 32 as a function of vessel diameter, d v , and f H ( ) D is a function relating viscosity to hematocrit, see Eq. 33.
is an empirical formula first presented by Pries, et al. 80 based on the behavior of blood in very small vessels.  The hematocrit distribution can change at bifurcations. Generally, the faster daughter branch (branch 1) has more hematocrit 18,90 . If the velocity ratio of two branches exceeds a threshold value, r th , all the blood cells enter the faster daughter branch 18  tional domain of the three-dimensional TME. A three-dimensional double hybrid continuous-discrete (DHCD) method was applied to solve the mathematical equations of the TME model (see Supplementary Material). In Fig. 2, a schematic of the multi-scale TME including subcellular, cellular and tissue scales is shown. Two separate cellular lattices with 3 × 3 × 3 nodes are constructed to determine the location of each tumor and endothelial cell, and a finite difference mesh with the same size as the lattices is selected to solve continuous parts of the DHCD (see Supplementary Material). The distance between adjacent nodes is 50 µm, which means we can assume each cell is like a sphere with diameter 50 µm or a cube with side 50 µm. The equations for the continuous part of the DHCD (all model equations except Eqs. (18) and (19) for tEC and TC, respectively) are normalized, discretized, and solved by a finite difference method (FDM) to determine the spatiotemporal distributions of TME biochemical and biomechanical factors. Based on these distributions, the motion probabilities of each TC and tEC are calculated to determine tumor progression in the TME (see Supplementary Material online for detailed computational approach and model setup). The parameters of the model used for computational results have been listed in Table S of Supplementary Material. The mathematical model is completely robust and independent of the selected parameters and thus can be used to simulate various tumor types with specific conditions. We demonstrate the robustness of model by changing some parameters, which can represent some of the most likely clinical problems of cancer patients. For the lattice resolution, we checked the mesh independency by performing simulations on lattices 20% smaller and 20% larger than the 50 µm used in the present study. Reducing the lattice size by 20% changed results www.nature.com/scientificreports www.nature.com/scientificreports/ by less than 10%; therefore, to reduce computational cost we chose the 50 µm lattice to present the results. A 50 µm spacing is also consistent with the average size of cancer and endothelial cell.
initial and boundary conditions. The initial concentrations of glucose, oxygen, carbon-dioxide, and ECM were assumed to be homogeneous and equal to c g ch , c o ch 2 , c co ch 2 , and c e ch , respectively. The initial concentrations of VEGF, Ang-1, Ang-2, and MMP are assumed to be zero. The free and active VEGFR-2 and Tie-2 are also homogeneously initialized to zero. At the boundaries of the computational domain of the TME, a Dirichlet boundary condition was used for each agent with value equal to its initial concentration. The TME was seeded with five tumor cells located at the center of the computational cube and a hypothetical circular primary vascular network with a radius approximately 5 mm. The locations of initial sprouts on the circle of primary vessels are determined based on VEGF concentration but spaced randomly according to NOTCH induction more than 50 µm apart 21 . The biomechanical factors IFF velocity and IFP, intravascular blood flow velocity and pressure, and WSS were initially set zero in the entire computational domain and boundaries. For these parameters, a Dirichlet boundary condition with zero value was set on all boundaries of the TME domain. At the inlet of the neo-vessels connected  www.nature.com/scientificreports www.nature.com/scientificreports/ to the primary vessels, we chose p lum = 30 mmHg, consistent with the range of primary vessel pressures reported in the literature 19,21,50,56 . To cover all three biological stages of tumor progression including ATG, angiogenesis (ATG-to-VTG), and VTG, we simulated 35 days of tumor growth and angiogenesis.

Results and Discussion
Morphology of the tMe. To analyze the various stages of tumor growth and angiogenesis, we selected days 10, 15, 20, 25, 30, and 35, which span ATG, angiogenesis (ATG-to-VTG), and VTG (Fig. 3) (see Movie S1 for a video of three-dimensional tumor growth and angiogenesis). Before day 10, tumor cells are supplied by nutrient diffusion symmetrically, and they expand as a symmetric semi-spherical structure. Between days 10-15, VEGF produced by the tumor creates a gradient that reaches the primary vessel and causes sprouting angiogenesis to generate neo-vessels. The neo-vessels elongate, branch, deform and extend toward the tumor during days 15-25. After day 25, the ATG-to-VTG transition is completed as neo-vessels penetrate into the tumor. During VTG (days [25][26][27][28][29][30][31][32][33][34][35], heterogeneities emerge because of the non-uniform neo-vasculature and the resulting asymmetric distributions of nutrients and growth factors. This causes spatial variations in cancer cell proliferation, resulting in distortion of the semi-spherical geometry. This transition from symmetric to asymmetric tumor shape mimics the evolution of an actual tumor. A qualitative comparison between our results and two biological observations of solid-tumor morphology 91,92 is illustrated in Fig. 4. Panels a-c show the initial stage of growth, with ATG transitioning to VTG, with very low density of angiogenic neo-vessels. Subsequently, VTG is characterized by a dense network of tumor-penetrating neo-vessels ( Fig. 4d-f). In agreement with in vivo observations, our simulations demonstrate that VTG can be characterized by an asymmetric and more invasive morphology compared to the symmetric, spheroidal morphology of ATG. In addition to recapitulation of three major stages of tumor progression (ATG, angiogenesis and VTG), deconstructing a reliable morphology of a solid tumor is the most important goal of a three-dimensional computational model because asymmetries and heterogeneities create complex chemical and mechanical gradients that greatly influence tumor biology and physiology.
We next quantified the number of angiogenic neo-vessels and their diameters at different times during tumor growth (Fig. 5a). Nascent neo-vessels have diameters less than 5 µm, while subsequent growth and maturation increase vessel diameters and cause lumenogenesis, all regulated by environmental stimuli. Vessel diameters can also decrease if there is a lack of positive angiogenic stimuli (see Eqs. 22a-d and 23). Deformation of the compliant neo-vessel can be caused by tumor growth-induced solid stresses, which are generally opposed by transvascular fluid pressure gradients (see Eq. 24). Vessel integrity can also be compromised by low levels of WSS or high levels of VEGF, resulting in local changes in vessel diameter. During ATG, some neo-vessels grow to as large as 50 µm in diameter with the majority falling in the diameter range of 40-45 µm (Fig. 5a). During VTG, branching of the neo-vessels increases near the tumor boundary and inside the tumor where VEGF concentration is high (Fig. 5a). This is because the VEGF induces sECs of neo-vessels to differentiate into tECs and generate new branches. Behind the leading edge of angiogenic sprouts, vessel maturation causes an overall increase in vessel diameters, with some reaching 65 µm in diameter. In contrast, further flow-mediated vessel remodeling from Day 30 to Day 35 causes a reduction in the largest diameters from 65 µm to 60 µm (see Eq. 23).
In addition to vessel morphology, we can examine the evolution of cancer cell phenotypes in the simulations (Fig. 5b). During ATG, active cells proliferate, but nutrient limitations in the tumor center cause some cells to become quiescent starting around Day 10. As angiogenic neo-vessels approach the tumor, but before tumor vascularization, the active cells have completely converted to quiescent cells because of insufficient nutrients, and the total number of cells plateaus; later in the transition from ATG to VTG, some quiescent cells become necrotic due to chronic nutrient deprivation (see Eqs. 16 and 17, and section 2.2). The number of necrotic cells increases until nutrient concentrations start to rise again due to tumor vascularization. Interestingly, the necrotic www.nature.com/scientificreports www.nature.com/scientificreports/ core persists due to rapid growth at the tumor periphery, which depletes nutrients before they can reach the deeper cells. During VTG, vessels bring more nutrients, and active cells begin to proliferate. After a time lag of a few days, the quiescent cell population increases in a similar fashion, and finally the necrotic cell population starts to increase again. The total number of cells linearly increases during ATG, plateaus during angiogenesis (ATG-to-VTG), and then exponentially increases during VTG. All these four tumor population dynamics agree with experimental observations and computational analyses 19,30,55,93 and are typical of solid tumors initiated far from primary-vessels 43,56 . Distribution of angiogenesis-associated tc vitality. We next examined the distributions of oxygen, glucose, CO 2 , and cellular vitality during tumor growth and angiogenesis (Fig. 6). The transition of the tumor geometry from a symmetric, semi-spheroid to an asymmetric, amorphous mass is apparent in the results. The active and quiescent TCs located at the tumor periphery regions around the necrotic core -central spherical region consistent during days 25-35 -(see Fig. 6, "Vitality"), consume more oxygen and glucose and, consequently, secrete more CO 2 . The limited supply of oxygen combined with high oxygen consumption by the peripheral TCs results in insufficient oxygen inside the tumor. However, glucose can be carried by neo-vessels to the boundary, but not to the central regions. Glucose delivered by the neo-vessels is sufficient for the peripheral cells, but not the central core of the tumor (Fig. 6, "Glucose"). The production rate of CO 2 by the peripheral TCs is higher than its removal by intratumoral neo-vessels, so CO 2 is concentrated inside the tumor, as shown in Fig. 6, "Carbon Dioxide". Diffusion into the center with no vessels for convection/removal combined with inactivity of necrotic cells for production/consumption causes some accumulation of CO 2 , oxygen and glucose in the tumor necrotic core. Distributions of other biochemical agents and biomechanical factors within the TME are presented in Supplementary Material.

Effect of biochemical blood abnormalities on tumor progression. Because cancer patients often
have abnormalities in their blood chemistry due to tumor growth or treatments, we next used the model to predict how these abnormalities might affect evolution of the tumor and its vasculature. We simulated hyperglycemia by doubling the glucose perfusion rate and hypoglycemia by halving the glucose perfusion rate. Hyperoxemia and hypoxemia were simulated by doubling and halving the oxygen perfusion rate, respectively, and hypercarbia and hypocarbia were produced by halving and doubling the carbon-dioxide perfusion rate, respectively. These conditions were compared with each other, and with normal blood (Fig. 7). Hyperoxemia and Hyperglycemia resulted in similar neo-vessel distributions, with the fewest vessels of all the conditions (Fig. 7a). This is because tissues with high levels of oxygen and glucose do not need to recruit as many new vessels to supply the tumor. With hypoxemia and hypercarbia vessel densities are similar, but hypercarbia resulted in more neo-vessel maturation. Hypoglycemia and hypocarbia both resulted in maximum neo-vessel number, with similar distributions. The conditions that produced angiogenesis that most resembled normal blood in terms of sprouting and branching were hypoxemia and hypercarbia. Considering neo-vessel maturation, hyperoxemia produced results similar to www.nature.com/scientificreports www.nature.com/scientificreports/ normal blood. Examining the populations of vessels that completely matured, all cases except hypercarbia have approximately equal vessel distribution. It is noteworthy that these differences in angiogenic neo-vessel distributions emerge because each blood abnormality affects a specific aspect of cellular respiration, which affects oxygen availability to the TCs and, consequently, tumor-induced VEGF.
Before VTG commences, blood abnormalities have no significant effect on the numbers of active, quiescent, necrotic and, consequently, total tumor cell (Fig. 7b-e). However, at the end of the transition from ATG to VTG, small effects of the abnormalities -especially on necrotic cells in the hypoglycemia blood condition -can be seen. Hypoglycemia and hypoxemia reduce the active and quiescent cells in the same manner, but hypoglycemia increases necrotic cells more than hypoxemia. These results suggest that hypoglycemia has a more potent positive effect on tumor growth than hypoxemia. Overall, these abnormalities reduce the total tumor cells (Fig. 7e). Hyperoxemia, and to a lesser extent hyperglycemia, increase active, quiescent, necrotic, and thus total cells. Interestingly, the necrotic core increases in size in response to both low and high oxygen/glucose conditions. In the case of low oxygen or glucose, this is probably due to insufficient penetration of nutrient-carrying vessels into the core. In the case of high oxygen or glucose, it is because of nutrient depletion by the outer layers of tumor cells, which starves the deeper cells. As another factor for promoting tumor growth, hypercarbia slightly increases active and quiescent cells and more strongly increases necrosis due to accumulation of more carbon-dioxide inside the tumor, leading to cell toxicity and death. In our simulations, hypocarbia produces tumor growth similar to the normal condition.

Effect of blood pressure and interstitial hydraulic conductivity (IHC) on tumor progression. to
further test the robustness of our model and its ability to predict a range of TME parameters, we next investigated the effect of blood pressure on tumor growth and angiogenesis by changing the inlet intravascular blood pressure from 30 mmHg to 20 and 40 mmHg (Fig. 8a,b). Increasing blood pressure skews the angiogenic neo-vessel distribution towards more large neo-vessels with diameters 55-65 µm (Fig. 8c) and significantly increases solid tumor size (Fig. 8b). Increased blood pressure is commonly encountered in cancer patients during chemotherapy treatment especially anti-angiogenesis drug treatment 94 . Therefore, there is a competition between drugs and drug-increased blood pressure to control tumor size.
In another set of simulations, we varied IHC (K ins , see Eqs. 29 and 30) to see how tissue hydraulic conductivity might affect the tumor. We found that decreasing IHC results in fewer angiogenic neo-vessels and increased tumor size (total number of TCs; Fig. 8c,d). Consequently, this condition was most efficient in terms of tumor growth-performance (TGP), defined as the number of total TCs divided by the number of total angiogenic neo-vessels (see Fig. 9). Our results show that tumor growth is much more sensitive to IHC reduction than angiogenesis (Fig. 8c,d): The baseline IHC level and the abnormal cases with twice IHC, 10% increase and 10% decrease in IHC all have approximately the same vessel patterns (except a small difference in nascent and very small neo-vessels with diameters 5-15 µm). However, these different IHC levels produced completely different tumor sizes (Fig. 8d). Comparing tumor size at day 35, it is clear that decreasing the IHC of the TME can result in www.nature.com/scientificreports www.nature.com/scientificreports/ a bigger tumor, and there is a lag in growth if we increase IHC (Fig. 8d). These results highlight the importance of only one aspect of desmoplasia -IHC reduction 95 . However, desmoplasia involves other modifications to tissue structure, cell populations and biochemistry, so more detailed information and additional modeling is warranted.
Mechanistically, decreasing IHC increases interstitial fluid flow, in agreement with experimental measurements [95][96][97][98] . Increasing IFF velocity increases transport of oxygen, glucose and CO 2 , helping tumor cells to survive, proliferate and generate a larger mass. On the other hand, this would also help drug delivery, and a richer tumor with more oxygen concentration generates less VEGF, inducing fewer angiogenic neo-vessels (Fig. 8a).
We next compared the various blood chemistry, blood pressure and IHC conditions in terms of tumor growth-performance (TGP), defined as the number of TCs per angiogenic neo-vessel (Fig. 9). The results of various blood chemistry, blood pressure and IHC conditions have been normalized by the values of normal TME, BP = 30 mmHg, and K ins , respectively. For the biochemical blood abnormalities, the highest TGP was seen for hyperoxemia, while similar low values were produced by hypoxemia and hypoglycemia (Fig. 9a). Hyperglycemia also had a relatively high TGP, while hypercarbia and hypocarbia produced TGP values slightly higher and lower, respectively, than the normal condition. Increasing blood pressure increased TGP (Fig. 9b); therefore, high blood pressure can be a risk factor for cancer patients 94 . Finally, IHC had the most dramatic effect on TGP: the lowest and the highest TGP resulted from the highest IHC and lowest IHC, respectively (Fig. 9c). Another parameter study to show the effect of subcellular level components including ECM/MMPs, VEGF and angiopoietins on tumor growth performance has been presented in the Supplementary Material file.

Validation of tumor Growth and Angiogenesis Metrics and trends
Some model predictions were compared and validated against the experimental studies, summarized in Table 2 as well as Fig. 10. The size of the avascular solid tumor before the onset of VTG and the three different phases of tumor growth are two common experimental observations evident in most tumor subtypes. According to Table 2, rows 1-3, and Fig. 10, the calculated results of our model are in good agreement with experimental measurements in terms of avascular solid tumor size and the dynamics of vascularization and vascularized growth.
Based on clinical studies, hyperglycemia is one of the important risk factors for cancer patients; it has been shown to increase tumor growth, causing larger tumors compared with control (5-110%) (Vasconcelos-Dos-Santos, et al. 99 ). Our model reproduces this result, predicting a 53% increment in tumor size ( Table 2, Row 4). Similarly, as mentioned above, desmoplasia has been shown to affect tumor growth and   www.nature.com/scientificreports www.nature.com/scientificreports/ transport of nutrients and drugs. By decreasing interstitial hydraulic conductivity, desmoplasia increases interstitial fluid pressure (Mpekris,et al. 95 ), and this increases tumor size by ~28%. We could reproduce the desmoplasia effects; theses simulations resulted in higher IFP and faster tumor growth with desmoplasia (size increased by 25-57%; Table 2, Row 5, Row 6).
The model also correctly predicts important vascular parameters such as vascular segment length, bifurcation density and diameter of angiogenic neo-vessels (Table 2, Rows 7-9).

conclusion
We presented a systems biology mathematical model of three-dimensional, multi-scale tumor microenvironment. We demonstrated the capabilities of the model by investigating how different TME conditions influence tumor progression during the various stages of tumor growth, including avascular tumor growth (ATG), angiogenesis, and vascular tumor growth (VTG). The model accounts for major biochemical and biomechanical factors in the TME, and these are all coupled across the subcellular, cellular, and tissue size scales. The model uses a comprehensive agent-based approach for the cell-signaling pathways in order to couple the various biological scales, creating a robust model able to recapitulate tumor progression under a variety of environmental conditions. We demonstrated the robustness and versatility of the model by simulating several common patient-specific TME disorders including biochemical blood abnormalities (hyperglycemia/hypoglycemia, hyperoxemia/hypoxemia, and hypercarbia/hypocarbia), increased blood pressure (an adverse event seen in chemotherapy patients) and altered interstitial hydraulic conductivity. The simulation results show that hyperoxemia/hyperglycemia blood abnormalities, chemotherapy-increased blood pressure, and reduction of interstitial hydraulic conductivity produce the most aggressive tumors with high tumor growth performance (TGP) values; these would be predicted to produce the worst patient outcome. All these blood abnormalities are possible in cancer patients, especially those receiving chemotherapy or radiotherapy treatments. We also elucidated the effects of subcellular level components including ECM, MMPs, VEGF and angiopoietins on tumor progression. The simulations show that the maximum TGP is produced when MMP level is low and VEGF secretion is high, and the minimum TGP results from high MMP secretion rates. One of the goals of developing this three-dimensional model was to accurately estimate the time delay before the onset of VTG (and thus the time of increased malignancy). Once the initial solid tumor is detected clinically, this can be determined under different environmental conditions using the model. This type of model is ideal for investigating the delivery and pharmacodynamics of anti-angiogenic and anti-cancer drugs, and their combinations to determine more effective treatment strategies to delay or prevent the switch to VTG, or to more effectively target solid tumors already in VTG. The model agrees well with experimental observations of tumor growth and angiogenesis, reproducing key features of angiogenesis, cancer cell metabolism and tumor asymmetric growth.

Data availability
The supporting data are availabe in the online supplemental files entitled 'Supplementary Material' , and 'Movie S1' . Source files of the TME model are accessible in Zenodo (DOI: 10.5281/zenodo.3609333).