Towards integration of time-resolved confocal microscopy of a 3D in vitro microfluidic platform with a hybrid multiscale model of tumor angiogenesis

The goal of this study is to calibrate a multiscale model of tumor angiogenesis with time-resolved data to allow for systematic testing of mathematical predictions of vascular sprouting. The multi-scale model consists of an agent-based description of tumor and endothelial cell dynamics coupled to a continuum model of vascular endothelial growth factor concentration. First, we calibrate ordinary differential equation models to time-resolved protein concentration data to estimate the rates of secretion and consumption of vascular endothelial growth factor by endothelial and tumor cells, respectively. These parameters are then input into the multiscale tumor angiogenesis model, and the remaining model parameters are then calibrated to time resolved confocal microscopy images obtained within a 3D vascularized microfluidic platform. The microfluidic platform mimics a functional blood vessel with a surrounding collagen matrix seeded with inflammatory breast cancer cells, which induce tumor angiogenesis. Once the multi-scale model is fully parameterized, we forecast the spatiotemporal distribution of vascular sprouts at future time points and directly compare the predictions to experimentally measured data. We assess the ability of our model to globally recapitulate angiogenic vasculature density, resulting in an average relative calibration error of 17.7% ± 6.3% and an average prediction error of 20.2% ± 4% and 21.7% ± 3.6% using one and four calibrated parameters, respectively. We then assess the model’s ability to predict local vessel morphology (individualized vessel structure as opposed to global vascular density), initialized with the first time point and calibrated with two intermediate time points. In this study, we have rigorously calibrated a mechanism-based, multiscale, mathematical model of angiogenic sprouting to multimodal experimental data to make specific, testable predictions.


Introduction
There is now a mature literature on the mathematical modeling of tumor initiation, expansion, angiogenesis, and invasion [1][2][3][4][5][6][7][8][9][10][11][12]. These models (which may be discrete, continuous, or hybrid) seek to provide new strategies for understanding the underlying biology, and then use this knowledge to make predictions of the spatiotemporal evolution of the disease as well as its response to therapy [13][14][15][16][17][18][19][20]. A fundamental barrier limiting progress in the field, though, is the paucity of studies linking mechanism-based mathematical models with the appropriate experimental data [21][22][23][24][25]. One fundamental reason for this impasse is that most mathematical models require specific and quantitative data sets to calibrate model parameters; data types that are not frequently available. This is especially true in the field of tumor angiogenesis where the mathematical models demand data that are acquired at both high temporal and spatial resolution to resolve vascular dynamics [22,[26][27][28][29][30]. Furthermore, the models themselves may have a myriad of parameters yielding an ill-posed parameter calibration problem [31][32][33][34]. Still, these models, even if uncalibrated against experimental data, can serve as useful tools to guide the experimental study of tumor-induced angiogenesis. For example, Vilanova et al. developed a phase-field model of capillary progression with a discrete tip endothelial cell that guides capillaries based on tumor angiogenic factors and neighboring capillaries [35][36][37]. This model accounts for capillary regression due to anti-angiogenic therapies and qualitatively matches in vivo experimental measurements. Their approach was extended by Xu et al. who incorporated blood flow and initialized the resulting model with photo-acoustic imaging data characterizing the tumor and surrounding vasculature [38]. Travasso et al. developed a phasefield approach describing endothelial, stalk, and tip cells that compares differences in vessel morphology due to tip cell migration and stalk cell proliferation [39]. Their model's predictions were qualitatively compared to experimental data of vascular patterns in response to the production level of angiogenic factors. These efforts are to be commended as some of the first examples of linking mechanism-based, mathematical models of tumor angiogenesis with experimental data. The next step is to calibrate such a model with early time point data, and then make a prediction to directly compare model forecasts to data acquired experimentally at later time points.
Calibrating mathematical models with experimental data allows for the ability to both predict future model outcomes and more rigorously test and understand the biological phenomena in question [23,24,40,41]. Specifically, for tumor angiogenesis, the ability to calibrate parameters that dictate the development of tumor vasculature in a well-controlled, experimental platform would allow the systematic testing of various models of angiogenesis [21,22]. However, obtaining the appropriate data for this task remains a difficult challenge. A substantial focus of experimental angiogenesis research is performed in vitro on 2D monolayers as well as in vivo xenograft models in animals. However, 2D systems do not mimic the complex behaviors observed in vivo, while xenograft systems pose substantial barriers to acquiring the data with adequate spatial and temporal resolution to resolve tumor induced angiogenesis [22,42]. We do note other in vitro works in experimental angiogenesis [43][44][45][46][47]. In the present effort, we address the data demands via a 3D in vitro microfluidic platform [48]. This system has the advantage of emulating the complex architecture of the tumor microenvironment observed in vivo, as well as capturing the spatiotemporal behavior of the cells. The 3D in vitro microfluidic platform we utilized in this study consists of a collagen matrix seeded with inflammatory breast cancer (IBC) cells around a functional blood vessel to represent the IBC tumorvascular interface. Over time, the IBC cells secrete pro-angiogenic factors that diffuse and bind to the functional blood vessel, initiating angiogenic sprouting. Spatiotemporal measurements of cell proliferation and apoptosis, branch lengths and number, anastomosis, and lumen formation in the platform will provide key parameters that drive the mathematical model.
In this contribution, we aim to calibrate an agent-based model of tumor angiogenesis [49] with experimental data obtained from the in vitro vascularized tumor platform [48] described above to predict vasculature sprouting. First, we perform experiments that isolate key components in the process of angiogenesis and use these data to calibrate important parameters in mathematical models of increasing complexity. This sequential approach allows a systematic, stepwise, approach for linking mathematical models with experimental data across spatial scales. The steps (referred to as scenarios in the remainder of this work) we follow are: 1) calibration of tumor and endothelial cell number to time-resolved hemocytometry, 2) calibration of the secretion and consumption of vascular endothelial growth factor (VEGF) by tumor and endothelial cells, respectively, to time-resolved protein concentration measurements, 3) calibration of stalk cell growth rate against time-resolved measurements of angiogenic sprout length, performing a 4) global and 5) local spatiotemporal calibration of the hybrid model of tumor angiogenesis to confocal microscopy data acquired from the in vitro microfluidic platform, and 6) testing the predictive capabilities of the calibrated mathematical model at the global and local scales.  Fig 1A, we show the coupled ordinary differential equation (ODE) models to describe tumor (IBC3) and endothelial (TIME) cell number, as well as the production of VEGF by tumor cells and the consumption of VEGF by endothelial cells. These models are calibrated using cell number over time (Scenario 1) counted via a hemocytometer, and VEGF concentration (Scenario 2) measured via ELISA analysis, shown in Fig 1B. Fig 1C depicts the hybrid multiscale model, initialized from confocal microscopy images observed in our 3D in vitro microfluidic platform, shown in Fig 1D, and calibrated by subsequent time points. In our global analysis, shown in green in Fig 1C and 1D, we calculate global quantities of interest from the confocal microscopy images, after thresholding for area and intensity and projecting the data into a 2D plane. From this 2D image, we calculate the sprout length, vascular density, and vascular volume fraction, which are readily comparable to the angiogenesis model to calibrate and predict the stalk cell divide time (Scenario 3) and the distance between tip cells for activation (Scenario 4). After analyzing the model globally, we assess the ability of the model to recapitulate local vascular features. In our local analysis, shown in peach in Fig 1C and 1D, we select individual sprouts for calibration through a Dice and area threshold over time, and skeletonize the observed vessels (and observe adjacent tumor cells). These skeletonizations, along with the tumor cells, are used to initialize and calibrate the local stalk cell divide time. We then assess the ability of the angiogenesis model to locally predict the final time point (Scenario 5). Finally, we use the calibrated parameter distributions to predict sprout length, vascular density, and vascular volume fraction (globally) and the local structure of the angiogenic sprouts.

Cell culture and VEGF concentration. Human inflammatory breast cancer (IBC)
MDA-IBC3 cells and telomerase-immortalized human microvascular endothelial (TIME) cells were used in this study. Green fluorescence protein labelled MDA-IBC3 (a HER2+ IBC cell line) were kindly provided by Dr. Wendy Woodward (MD Anderson Cancer Center, Houston, TX) and mKate labeled TIME cells were a generous gift from Dr. Shay Soker (Wake Forest Institute for Regenerative Medicine, Winston-Salem, NC). MDA-IBC3 cells were cultured in Ham's F-12 media supplemented with 10% fetal bovine serum, 1% antibiotic-antimycotic, 1 μg/ml hydrocortisone, and 5 μg/ml insulin. TIME cells were cultured in Endothelial Cell Growth Medium-2 BulletKit (EGM-2, Lonza). All cell cultures utilized in this study were maintained at 5% CO2 atmosphere and 37˚C.
To determine the production of VEGF by MDA-IBC3 cells and consumption of VEGF by TIME cells, cells were seeded with a density of 30,000 cells/cm 2 and 10,000 cells/cm 2 ,  Fig 1A and 1B show the calibration of tumor and endothelial cell number and VEGF concentration to hemocytometer and ELISA measurements over time. Given the calibrated VEGF production and consumption rates, we inform our hybrid multiscale model, shown to Fig 1C, to confocal microscopy images of angiogenic sprouts, depicted in Fig 1D. We analyze the model globally, calibrating model parameters to summary statistics of the data, namely the sprout length and vascular density, and locally, calibrating the local stalk cell divide time to segmented vascular structures. This sequential approach, starting with protein concentration and cell number experiments to inform the hybrid model prior to integrating confocal microscopy images, allows us to utilize experimental data at multiple scales to inform the multiscale nature of the tumor angiogenesis model. respectively, in separate 12 well plates supplemented with EGM-2 media. Media was collected and replaced from the same 4 wells daily for measurement of VEGF concentration and cell growth was imaged with an IncuCyte Zoom (Essen Bioscience, Ann Arbor, MI) every 8 hours for 7 days. VEGF concentration was measured using a Human VEGF Quantikine ELISA Kit (R&D Systems, Minneapolis, MN). The cell number in each well was determined by trypsinization and counting of cells with a hemocytometer on days 1, 2, 3, 5, and 7. For each timepoint, four wells were sacrificed for cell counting and these were separate from the media collection wells for VEGF measurements.
2.2.2 Vascularized 3D in vitro microfluidic platforms. The in vitro 3D tumor microfluidic platforms utilized in this study were composed of extracellular matrix comprised of collagen type I seeded with green fluorescence protein labeled MDA-IBC3 around a hollow vessel lined with mKate labeled TIME cells housed in a polydimethylsiloxane scaffold as described in our published work. 7 mg/ml collagen solution, which has a similar stiffness to that of in vivo breast tumors, was used in fabricating extracellular matrix of the platforms [50][51][52]. MDA-IBC3 cells were seeded at a density of 1 × 10 6 cells/ml in the collagen solution and polymerized around a 22 G needle at 37˚C for 25 minutes. After polymerization, the needle was removed, and the resulting hollow vessel was filled with a solution of 2 × 10 5 TIME cells to form an endothelialized vessel lumen. Flow was introduced using a syringe pump system and a graded flow protocol was used to establish a confluent endothelium as we have previously published [48,51,[53][54][55]. The microfluidic platform is summarized in Fig 2. Fig 2A shows an example microfluidic platform with an axial cross-section depicting the functional parent blood vessel and cancer cells seeded throughout the microenvironment in Fig 2C. The platforms are connected to a syringe pump system, shown in Fig 2B, that allows for the continuous perfusion of media into the collagen microenvironment. Over time, the cancer cells release pro-angiogenic proteins, which diffuse through the microenvironment and cause tumorinduced angiogenesis. An example of these sprouts is shown in Fig 2D. 2.2.3 Image acquisition and processing. Platforms were imaged with a Leica TCS SP8 Confocal Microscope (Leica Microsystems, Germany) on days 3, 5, 7, 9, 11, 15, and 19. The 3D volume of angiogenic sprouts was imaged with z-slices acquired at intervals of 4.28 microns

PLOS COMPUTATIONAL BIOLOGY
Towards integrating multimodal data with a hybrid multiscale model of tumor angiogenesis and an in-plane resolution of 2.254 microns × 2.254 microns in the x-y plane (shown in Fig  2D). The confocal microscopy images were processed at each time point to ensure the same local structures of the microfluidic platform could be compared over the duration of the experiment. This was done via a reference mark. It was assumed that the diameter of the parent vessel did not deviate significantly over the course of the time points used for calibration. The parent vessel periphery was used as a reference to align the segmented region of sprouting in x-y space, depicted in Fig 2D. We assume that the vessels grow primarily in x-y space and that little growth is in the z direction (we comment on this assumption in the Discussion section). Vessel lengths were estimated at multiple time points using the built-in Leica measurement. We average signal intensity across 10 slices in the z direction (approximately 43 microns in total) and thresholded using area and intensity (to filter out migratory endothelial cells that are not a part of maturing vessels) to create a binarized vessel mask. This mask is used to calculate both the global and local vascular quantities for calibration. The global quantities of interest are the vascular density, calculated by summing the vascular mask along the length of the vessel, and the vascular fraction, calculated by dividing the total vascular mask area by the total area in the experimental computational domain. For local calibrations, we calculate the local vessel morphology (i.e., the structure of specific vasculature) over time.
To select specific vessels for local calibration, we developed an algorithm that calculates the Dice and percentage of overlap across time points. This allows for selecting regions that, when aligned, correspond to angiogenic sprouts growing over time and minimizes the effects of random endothelial cell migration. For the selected regions, we utilized a Zhang-Suen thinning algorithm [56] (available at https://github.com/linbojin/Skeletonization-by-Zhang-Suen-Thinning-Algorithm) for skeletonization after dilating the vessels using the 'strel' function in MATLAB (Mathworks, Natick, MA). Given the skeletonization for the model and the data at time t, we calculate the average distance between the model to the data centerline and from the data to the model centerline to use as the objective function for model calibration by using Python's SciPy Euclidean distance transform [57].

Overview of computational methods.
We calibrate three mathematical models for the purpose of predicting vascular response to the secretion of VEGF by tumor cells. (We first describe them qualitatively before introducing their mathematical representations.) The first and second are systems of ordinary differential equations that describe VEGF consumption by TIME cells and production of VEGF by MDA-IBC3 cells. These models describe the concentration of VEGF over time and the evolving cell type associated in production or consumption. The first model describes VEGF consumption of endothelial cells and has three model parameters: the consumption rate of VEGF by endothelial cells, and the growth rate and carrying capacity of endothelial cells. All of these parameters may be determined from the VEGF concentration and cell culture data. The second model describes the growth of MDA-IBC3 cells and their production of VEGF through three parameters: the VEGF production rate and carrying capacity, and the growth rate of the tumor cells. These parameters may also be determined from the VEGF concentration and cell culture data. We utilize the consumption rate of endothelial cells and the production rate of endothelial cells as inputs into the third model, a hybrid multiscale model of tumor angiogenesis. An agent-based model describes the tumor and endothelial cell dynamics and a continuum model governs the dispersion of VEGF. This hybrid model computes the spatiotemporal distribution of vascular sprouts and is readily comparable to the confocal microscopy images (and the calculated quantities of interest at the local and global scales) of the 3D in vitro microfluidic platform. These three models vary in complexity and allow us to isolate model parameters in a series of well-controlled experiments, and ultimately to test the predictive capabilities of the multiscale model at the global and local scales.
2.3.2 Model of VEGF production and consumption. The following two systems of mathematical models describe the consumption and production of VEGF by TIME and MDA-IBC3 cells, respectively. The first system of ODEs is: where l g e is the growth rate of endothelial cells, N e (t) is the number of endothelial cells at time t, θ e is the carrying capacity of endothelial cells, [VEGF] is the concentration of VEGF, and l c e is the consumption rate of endothelial cells. This system assumes that each endothelial cell consumes VEGF at a constant rate, VEGF is not produced by endothelial cells, N e (t) follows logistic growth with carrying capacity θ e , and VEGF decay is negligible. To model the VEGF secretion by tumor cells we write: where l g T is the growth rate of tumor cells, N T (t) is the number of tumor cells, l p T is the production rate of VEGF by tumor cells, and θ V is the carrying capacity of VEGF. Since both N e (t) and N T (t) are measured daily, we calibrate l g e and l g T and utilize the model to predict cell numbers between experimental measurements. To directly compare model predictions to the cell culture and protein concentration data, we take N T (t) and N e (t) to be MDA-IBC3 cells and the generic endothelial cells to be TIME cells, respectively.

Hybrid multiscale model of tumor angiogenesis.
We utilize the hybrid multiscale model of tumor induced angiogenesis developed in [49] to compare the global and local quantities of interest with the confocal microscopy data in [48]. We now summarize the salient features of this model and refer the interested reader to [49]. The multiscale model describes tumor and endothelial cell dynamics using an agent-based model, and the evolution of nutrients and vascular endothelial growth factor (VEGF) using a continuous partial differential equation (PDE) model. Tumor cells act as individual agents and may reside as a proliferative, quiescent, hypoxic, or necrotic cell. They transition from quiescent to proliferative stochastically with a linear dependence on nutrient concentration. Tumor cells uptake nutrients, grow, and proliferate until the supply of nutrients is depleted (i.e., less than a threshold, σ H ), at which point they become hypoxic and secrete VEGF. The VEGF diffuses across the domain and eventually contacts nearby vasculature. The blood vessels, comprised of endothelial cells modeled as agents, uptake VEGF and become activated when the concentration of VEGF is greater than a threshold, α A , and the distance between other tip cells is greater than d tip . The activated endothelial cells then transition to tip cells, which are responsible for the migration of angiogenic sprouts up the concentration gradient of VEGF due to a chemotactic force, where f v is the VEGF force coefficient, d ij the distance between tip cell i and the neighboring endothelial cell j, R a N is the nuclear radius of cell α, R α is the cytoplasmic radius of cell α, and the VEGF potential function rψ is given as The proximal endothelial cells, termed stalk cells, proliferate after time d sc and allow the elongation of the sprouts. The effective sprout elongation rate can be calculated by summing the stalk cell divide time, d sc , and the stalk cell growth time, g sc , divided by the diameter of the stalk cell, 2R sc . Mechanical adhesion and repulsion forces maintain the structural integrity of the vessels and allow the formation of the lumen. Importantly, this formulation of VEGF force properly scales the chemotactic forces acting on the tip cell to maintain vessel structural integrity. Once the sprouts mature, they form complex networks and establish blood flow, enabling the delivery of new nutrients to the tumor. The VEGF and nutrient dynamics are governed by the following reaction-diffusion equations: where D n and D v are the nutrient and VEGF diffusion coefficients respectively, Λ n is the nutrient uptake rate of tumor cells, Λ v is the VEGF consumption rate of endothelial cells, Γ n is the nutrient production rate from endothelial cells that are part of anastomotic loops, Γ v is the VEGF production rate of hypoxic tumor cells. These source and sink terms, which couple the cellular agent-based model to the continuous PDE model at the tissue scale, are defined as follows: where the subscripts pq, h, l, t, s, and e specify proliferative and quiescent tumor cells, hypoxic tumor cells, looped endothelial cells, tip endothelial cells, stalk endothelial cells, and endothelial cells, respectively. Also, ϕ α (x, t) is the volume fraction of cell type α 2 {pq, h, e, t, s, p}, at position x and time t, g c a is the consumption rate of the cell type α, and γ α is the production rate of the cell type α. The volume fraction of each cell type is calculated by computing the total area of each cell type within each element mesh.
In the analysis that follows, we assume that all tumor cells produce VEGF equal to the value calibrated from the cell culture and VEGF concentration computational-experimental scenario (see sections 2.4.1 and 2.4.2) but stop producing VEGF whenever they are within the threshold distance from the sprouting vessel. Thus, we are using this distance as a surrogate for nutrient delivery and can perform our computational-experimental studies independent of nutrient/ hypoxia status. For the calibration of the stalk cell divide time in the global analysis, we use a simplified model to calculate the length of the angiogenic sprout to compare with the experimental data and assume the growth of the sprout is independent of VEGF concentration. This simplified model is identical to the full model but with a significantly reduced domain and only one tip cell that grows perpendicular to the parent vessel, allowing removal of the VEGF field, thereby dramatically reducing computational cost while maintaining the angiogenic growth dynamics.

Bayesian calibration.
Bayesian inference provides a statistical framework for calibrating model parameters by accounting for both uncertainties in the mathematical model and the experimental data. Thus, it is an excellent methodology for determining uncertainties in model predictions of experimental outcomes [58][59][60][61]. We now summarize the salient features of Bayesian parameter calibration.
Beginning with Bayes' Rule, where m denotes the model parameters, y is the experimental data, π posterior (m|y) is the distribution of model parameters given the experimental data, π prior (m) is the prior distribution of the model parameters, π likelihood (y|m) is the conditional probability of the experimental data given the model parameters, and π evidence (y) is a normalizing factor. These probability density functions (PDFs) provide a holistic description of the system (and therefore) the parameter uncertainty. Since π evidence (y) is a normalizing factor, we can rewrite the posterior distribution as, The prior information of the parameters, π prior (m), is chosen on a problem specific basis. If only the bounds a and b of the parameter are known, then the prior is selected as a uniform distribution, π prior (m i ) * U(a, b). Since only the bounds are known in this work, we exclusively use uniform priors. This means that π prior (m i ) is constant and we can further simplify Bayes' Rule to p posterior ðmjyÞ / p likelihood ðyjmÞ: ð11Þ Assuming independent and identically distributed data, where σ is the standard deviation, N is the total (combination of spatial and temporal) number of data points, y i data is the experimental data, and y i model is the model prediction given parameters m. This can also be written as, where ϕ(m) is the misfit to the data. With a given data misfit, Eq (13) is straightforward to calculate. The hyperparameter, σ, characterizes the uncertainty of the experimental data and the inadequacy of the model [62,63]. The posterior distribution for the parameters can be learned by utilizing sampling or quadrature methods to compute π likelihood (y|m) at specific values of the parameters. Therefore, this framework allows us to use experimental data to learn the updated distribution of parameters. These parameter distributions can then be used to propagate uncertainties in the model predictions for a robust description of prediction uncertainty.
We also specifically note that the maximum likelihood parameter values, m � , can be obtained from π likelihood (y|m) by taking the maximum of π likelihood (y|m) with respect to m. This can be particularly useful in multilevel calibration schemes, as some parameters from one scenario are useful in later scenarios. For example, we calibrate the growth rate of endothelial cells, l g e in Eq (1), to use as an input for solving VEGF concentration in Eq (2) where we calibrate l c e without considering the uncertainty in l g e (the consumption of VEGF by endothelial cells), which will be used in subsequent scenarios).

Model calibration and prediction scenarios
We now define the parameter estimation problems in the ODE models (i.e., Eqs (1)-(4)) and the hybrid multiscale model (i.e., Eqs (5)- (8)); the parameters to be calibrated are listed in Table 1. The calibration scenarios results and analyzes are provided at https://github.com/ CalebPhillips5/Vascular_calibration.

Scenario 1: Calibration of tumor and endothelial cell growth.
First, we utilize the hemocytometry data to calibrate the parameters from Eqs (1) and (3). We assume tumor and endothelial cells growth is independent of the VEGF concentration and can therefore separate the coupled systems into two separate calibrations. The parameters calibrated in this scenario are the growth rate of tumor and endothelial cells (l g T and l g e , respectively), and the carrying capacity of endothelial cells (θ e ). We then define the parameter estimation problem for parameters m in both models as: Find m ¼ ½l g T ; l g e ; y e � such that is minimized. N data (t i ) and N model (t i , m), denote the measured and model prediction of tumor or endothelial cells respectively, and N days is the number of measurements. We use Eqs (1) and (3) with calibrated parameters as input functions required to infer the VEGF production and consumption rates for tumor and endothelial cells, respectively, in Eqs (2) and (4). We compare the relative error in tumor and endothelial cell number to show the model goodness of fit.

Scenario 2: Calibration and prediction of VEGF concentration.
For the VEGF calibration problem, we utilize the ELISA data to calibrate the production and consumption rates of VEGF, as well as the carrying capacity of VEGF in the dish (Eqs (2) and (4)). We then define the parameter estimation problem in both models as: is minimized, where V data (t i ) and V model (t i , m), denote the measured and predicted concentration of VEGF, respectively. In this scenario, we employ a naive Bayesian quadrature (i.e., selecting quadrature points uniformly across the domain) approach to minimize ϕ(m) with respect to the parameters m. Since the VEGF production and consumption rates are assumed to transcend the experimental scenario (as characteristics of the cells as opposed to characteristics of the environment), we will sample from these distributions moving forward in both the prediction of the VEGF concentration as well as the prediction of the vascular ABM (section 2.4.4). To easily generate samples, we approximate the PDFs as Gaussian distributions. We calculate the relative error in VEGF concentration for tumor and endothelial cells and show the propagated uncertainty through time.

Scenario 3: Calibration and prediction of global vessel length.
We utilize confocal microscopy images to measure the length of angiogenic sprouts over time. This provides a 1D length measurement that can be used to calibrate parameters that govern the growth and development of the sprouts. We selected vessels that can be observed on day 3 and subsequently tracked during the remainder of the experiment to ensure consistent (i.e., the growth of the same sprouts) sprout growth. This is done via a reference mark on the microfluidic platform to align the slices across time. The key parameters are the VEGF force coefficient, f v , and the stalk cell divide time, d sc . We then define the parameter estimation problem as: is minimized, where L data (t i , s j ) and L model (t i , m) are the observed and predicted length of sprout s j at time point t i , with parameters m at time point t i . The number of sprouts, N sprouts , was chosen such that the same sprouts were tracked throughout the duration of the experiments. The calibration is solved using a Metropolis-Hastings algorithm available in the QUESO (Quantification of Uncertainty for Estimation, Simulation, and Optimization) library [64]. The prediction of L data (t i , s j ) is solved by running the forward model with parameters selected from specific samples used in the Metropolis-Hastings algorithm for calibration and calculating the mean and standard deviation of the resulting quantities of interest, L model (t i , m). We report the relative error between the observed and predicted vessel length over time.

Scenario 4: Calibration and prediction of vascular density.
To conclude our global analysis, we calculate the vascular density moving away from the parent vessel across time. This provides a global description of the total vascular sprouting to infer the probability distribution function of the distance between new tip cells. We then define the parameter estimation problem as: Find m = [d tip ] such that is minimized, where D data (t i , s j ) and D model (t i , m) are the experimentally estimated (from the confocal microscopy data) and model predicted vascular density, respectively, at time t i and voxel s j , with N x being the number of voxels N x . As the parent vessel in the microfluidic platform is not perfectly straight, and some endothelial cells may migrate short distances from the parent vessel without forming mature vessels, we only utilize voxels s j that are beyond 15 voxels (> 35 μm) away from the parent vessel. This mitigates the effects of endothelial cells that do not form mature vessels which is the focus of our model. This calibration scenario is solved by using a naive quadrature approach as the parameter space is only one dimensional and the computational complexity of the model is prohibitive for sampling. A quadrature approach allows us to utilize parallel resources to run thousands of simulations concurrently, while sampling requires several thousand model runs in serial. For the calibration, we use the day 3, 5, and 7 data and the calibrated probability distribution for distance between new tip cells to predict the vascular density at day 9. We note that the computational and experimental domains used in this scenario are identical (676.2 × 3155.6 μm 2 ). We assess the predictive capabilities of the hybrid multiscale model of tumor angiogenesis at the global and local scales by comparing the model prediction, using calibrated parameter distributions, to the observed data at the final time point in each experiment. The uncertainty in the predictions of vascular density and vascular volume fraction is quantified globally in two cases: the 1-parameter case, using only the parameter distribution for distance between new tip cells (d tip ; values sampled from a Gaussian distribution calibrated in Scenario 4), and the 4-parameter case, using the parameter distributions for VEGF production and consumption (l p T and l c e ), stalk cell divide time (d sc ), and the distance between new tip cells (d tip ; values are sampled from Gaussian distributions calibrated from Scenarios 2-4). This is accomplished by sampling from the Gaussian fits of the parameter distributions and running the model with the sampled model parameters to calculate the quantities of interest. The 4-parameter case accounts for uncertainty (and therefore the error) of the calibrations performed in Scenarios 2-4, as well as, the propagation of the uncertainty. In the global analysis, we predict the vascular density and vascular volume fraction, while in the local analysis we predict the centerlines of the angiogenic sprouts.
2.4.5 Scenario 5: Calibration of local stalk cell divide time to predict vessel morphology. We employ the confocal microscopy images (see 2.2.3 Image acquisition and processing) to calibrate the hybrid multiscale model to local vascular structures over time. To capture the spatial differences between the model and the data, the optimization problem involves comparing the average distance between the centerlines of the vessels observed in the data and computed by the model. Since the local analysis focuses on specific vascular sprouts, we initialize the mathematical model from the data at the first time point (day 3 of the experiment), calibrate the growth of the stalk cells from days 5 and 7, and predict the vascular structure at day 9. Since the model and data are not of the same type (i.e., the data is a binary field and the model is a list of cells interacting via forces), initialization is quite complex. We initialize the calibration by setting the direction of VEGF gradient to closely resemble the true vessel skeletonization at day 3. We then define the parameter estimation problem as: Find m ¼ ½d � sc � such that is minimized, where D MD (t i , x j , m) is the average distance between the voxels (with number of voxels N x ) on the centerline of the model to the data and D DM (t i , x j , m) is the average distance between the voxels on the centerline of the data to the model. Across days 5-7 and days 7-9, only two regions had a Dice score above 0.5 and a percentage of overlap above the average of all regions analyzed. To assess the model prediction of local vessel morphology, we implemented the following steps: 1) take 1000 samples from the Gaussian distribution calibrated from Scenario 5 (as well as the VEGF consumption and production rate calibrated in Scenario 2), 2) run the forward model with those 1000 parameter realizations, 3) calculate the centerlines of the forward model solutions, 4) sum over the 1000 centerlines to create a vessel centerline heat map, and 5) threshold the heat map at 10, 50, and 100 to provide the prediction envelope for 1%, 5%, and 10% of the simulations (i.e., voxels covered by 10, 50, and 100 predicted vessel centerlines). The quartiles of the centerline prediction of each region are calculated from MATLAB's 'quantiles'.

Calibration Scenarios 1 and 2: Cell number and VEGF concentration
Fig 3 provides an overview of how cell number and ELISA analysis are integrated with the ODE models for the calibration of VEGF production and consumption rates. Cell number was calculated via hemocytometer at 24, 48, 72, 120, and 168 hours post seeding. The first time point initialized the TIME (Fig 3B) and IBC3 (Fig 3C) growth models and the subsequent four time points were used to calibrate the models for cell number over time (i.e., Eqs (1) and (3)). The number of cells are independent of VEGF concentration allowing the tumor and endothelial cell number calibration to be carried out first with N e (t) and N T (t) used as input functions into the production and consumption models of VEGF, respectively. The fit of the growth models is shown in Fig 3B and 3C with growth rates of 8.9 × 10 −3 h −1 and 3.7 × 10 −2 h −1 for the IBC3 and TIME cells, respectively, and a carrying capacity of 2.67 × 10 5 for the TIME cells. In the ELISA experiments, every 24 hours for 7 days, VEGF concentration is measured, and media is replaced; note that this action causes the total concentration of VEGF in the plate to return to baseline (i.e., 1100 pg/mL) after each measurement. The VEGF experiments for TIME and IBC3 cells are carried out in different plates, so the calibration of each cell line is not coupled. We utilize naive Gaussian quadrature to calculate the probability distribution of VEGF production by TIME cells, as well as the VEGF consumption and carrying capacity of the IBC3 cells. The maximum likelihood values were obtained using MATLAB's 'lsqnonlin' function, and we then utilized a parameter sweep using 100 values for each of the three parameters (total of a million points) taken about the maximum likelihood values. The resulting posterior distributions are shown in Fig 3D. The fit Gaussian distributions are given by l c e � 10 À 7 � Nð9:1; 2:3Þ and y e � Nð4035; 164:6Þ pg; where N(a, b) denotes a Gaussian distribution with mean a and standard deviation b. In Fig  3E, the IBC3 and TIME VEGF concentrations (with standard deviations) are shown as the red and blue dots, respectively. The maximum likelihood value for each parameter (see Section 2.3.4 is l c e ¼ 8:9 � 10 À 7 and y e ¼ 4035 pg: The means (and associated σ prediction envelope) of the best fit of TIME and IBC3 VEGF models are shown in blue and red, respectively. The calculated distributions and maximum likelihood values of VEGF production and consumption rates will be used in subsequent calibration and prediction scenarios described below in Sections 3.2-3.4. The relative errors between the models and experimental measurements are summarized in Table 2.   Table 2. Error between the ODE model output and cell culture and VEGF concentration data. All errors are relative errors between the model output and experimental data. The dashes (-) denote where the data is not available to calculate the quantity, and the asterisk ( � ) denotes when the data records a value of less than 2% of the baseline concentration of 1100 pg/mL, thereby precluding computing the relative error (i.e., it would be near infinite).

Days
Error IBC Error TIME Error IBC3 VEGF Error TIME VEGF   N(147.8, 18.6). The Gaussian fit of d sc is utilized moving forward (specifically, in Scenario 3.4) to assess the uncertainty in the model prediction of vascular density and volume fraction. As the increase in vessel length is only due to the division of stalk cells, we can ignore the uncertainty contribution of the VEGF force. This constant acts to maintain physiological sprouting that does not break the vessel but balances the forces to allow vessel elongation; i.e., if this force is higher than the adhesion between endothelial cells, the tips cells would disconnect from the sprout. Therefore this uncertainty contribution does not affect the resulting vascular network.   Fig 5B. Fig 5C depicts the calibrated PDF and Gaussian fit of the distance between new tip cells. Fig 5D-F show the density calculated from the left and right side of the microfluidic platform shown in blue and magenta, respectively, and the best fit of the model shown in red for days 3, 5, and 7, respectively. Specifically, the model best fit is calibrated using days 3 and 5 and then temporally evolved to compare day 7 with the observed experimental data. The vertical dashed lines represent the location away from the parent vessel that we begin to use for calibration (i.e., we utilize the voxels to the right of the dashed vertical line, see Scenario 4, Section 2.4.4). The relative error of the best fit in the calibration for days 3, 5, and 7 are 23.5%, 11.1%, and 18.5% respectively.

Scenario 4: Calibration and prediction of vascular density
https://doi.org/10.1371/journal.pcbi.1009499.g005 the vessel mask along the length of the vessel. Since the left-and right-hand side are not able to be compared directly with the computational domain, we align the vascular density calculated from the image and the model initial conditions (see the blue and magenta lines in Fig 5D-5F for an example alignment). In Fig 5C, we show the calibrated PDF of the distance between new tip cells in blue and the Gaussian fit of the PDF, which may be readily sampled for analyzing the uncertainty in the prediction of the model. In Fig 5D-5F, we show days 3, 5, and 7 of the calculated density of the left and right sides of the microfluidic platform (blue and magenta, respectively), and the model best fit obtained during calibration over days 3 and 5 (red). The day 7 model best fit is the temporally evolved solution of the model best fit calibrated from days 3 and 5. The vertical dashed line represents the location from which we begin to make use of the vascular density observed in the microfluidic platform for the calibration. That location is chosen to minimize the effects of migratory endothelial cells close to the channel wall that do not form mature angiogenic sprouts (see Section 2.4.4. The Gaussian fit of the PDF is d tip * N(243.9, 18.1) microns. The error in the calibration and prediction is shown in Table 4. We utilize 1000 samples from the probability distribution of the distance between new tip cells, as well as the probability distributions obtained in scenarios 2 and 3, to rigorously quantify the uncertainty in the predictability of the model. The prediction results are shown in Fig  6. In Fig 6A-6C, we show the 95% prediction confidence intervals propagating uncertainty using one parameter distribution for distance between new tip cells shaded in light red and using four parameter distributions (distance between new tip cells, VEGF production and consumption rates, and stalk cell growth rate) shown in light blue. The vascular density calculated from the left and the right sides of the microfluidic platform are shown in blue and magenta, respectively, and the model best fit is shown in red. The vertical black dotted line denotes the distance away from the vessel that we begin to use for calibration. The total uncertainty (summed over space) was 97.4%, 66.2% and 45.2% higher in the four-parameter case than the 1-parameter case on days 3, 5, and 7, respectively. Fig 6D shows the calculated average vascular fraction from the data in black, and the predicted average vascular fraction for the one parameter (blue) and four parameter (red) uncertainty analyses. In Fig 6E, the vascular fraction of the data, 1-parameter, and 4-parameter scenarios are shown in blue, orange, and yellow, respectively.

Scenario 5: Calibration of local stalk cell divide time to predict vessel morphology
While sprout length, vascular density, and volume fraction over time are reasonable metrics for analyzing global features of vascular structure, additional local analysis is needed to assess the model's ability to recapitulate local features in the vascular network. Fig 7 shows the local region selected from the microfluidic platform for the calibration of stalk cell growth rate. Table 4. Error in calibration and prediction of vascular density. We compare the relative error in the calibration best fit, the 1-parameter, and 4-parameter prediction cases.

Days
Calibration best fit % error   Fig 6A-C show the uncertainty in the model prediction assuming uncertainty in one parameter calibrated from scenario 4 (shaded in light red, outlined in black) and four parameters (shaded in light blue, outlined in black) calibrated from calibration scenarios 2, 3, and 4, along with the mean of the data. In the 1-parameter case, we consider only the distance between new tip cells with a PDF shown in Fig 5C. In the 4-parameter case, we consider distance between new tip cells, VEGF production and consumption rates (shown in Fig 3D), and stalk cell growth rate (shown in Fig 4B). Fig 6D presents the volume fraction calculated from the data in black, and the predicted volume fraction from the 1-parameter and 4-parameter case in light red and blue, respectively. Fig 6E depicts the prediction of volume fraction compared to the data (blue), and the 1-parameter and 4-parameter cases in orange and yellow, respectively, with the corresponding standard deviations shown in black. The average error in vascular volume fraction is 20.2% and 21.7% of the 1-parameter and 4-parameter case, respectively. We also note the drastic increase in uncertainty moving away from the parent vessel in the 4-parameter case (Fig 6A 120-180  Confocal images are taken on days 3, 5, 7, and 9, and processed in MATLAB using a Dice score and area threshold to select regions for local calibration (see 2.2.3 Image acquisition and processing). Once selected by our sorting algorithm, we calculate the skeletonization of the local vascular structures using a Zhang-Suen parallel thinning algorithm [56]. We use the day 3 skeletonization, shown in Fig 8A to initialize the model (Fig 8D) by segmenting the tumor cells from the surrounding region and calculating the trajectories of the tip cells shown in the data to formulate the initial conditions of the ABM that match the skeletonization of the data. The stalk cell growth rate is calibrated by running the model forward and comparing the centerline distance between the data at days 5 and 7 (Fig 8B and 8C, respectively) and the vasculature predicted by the model with the maximum likelihood stalk cell growth rate, shown in Fig  8C and 8F. The resulting probability distribution function of stalk cell growth rate (blue) and the Gaussian fit (black) are shown in Fig 8G. The maximum likelihood value for stalk cell growth rate in this scenario is 10.86 hours and the Gaussian fit is d sc * N(10.85, 2.4) hours. We then take 1000 samples from this Gaussian distribution and calculate the resulting vasculature for days 5, 7, and 9 as shown in Fig 8. In Fig 9, the columns denote day 5, 7, and 9, respectively. The rows depict the 1%, 5%, and 10% prediction envelope (top to bottom, described in 2.2.3 Image acquisition and processing). In Fig 9B, the data and the model calculated centerlines form an anastomosis, predicted in 3.3% of simulations. Anastomosis is also predicted on day 9 in Fig 9C, 9F and 9I in 16.7% of simulations, recapitulating the vessel structure observed in the data at day 7. In Fig 9J, the prediction of the average centerline distance from model to data is calculated from the 1000 samples of the distribution in Fig 8G. Fig 9K presents the normalized average centerline distance from the data to model (normalized to the length of the sprout in day 3). The prediction quartiles of this region and another (shown in S1 Appendix) is shown in Table 5. Finally, Table 6 shows the calibrated model parameters and their Gaussian fit. In particular, we note that s e is the global sprout elongation rate which is independent of cell diameter and is calculated from the stalk cell divide time and the growth time of stalk cells, g sc .  Fig 8A-C show the centerlines calculated from the model best fit (red), the data (green), and the overlap (yellow) on day 3 (used to inform the initial conditions), day 5, and day 7 (used to calibrate the stalk cell growth rate). The ABM is shown in Fig 8D-F with tip cells in green, stalk cells in cyan, endothelial cells in red, tumor cells releasing VEGF in orange, and tumor cells not releasing VEGF in blue. In Fig 8G, we show the calibrated PDF and the Gaussian fit of the stalk cell growth rate. The best fit model recapitulates the general structural features of the data without allowing additional sprouts to form. We also note that from day 7 (depicted in Fig 8F) Fig 9B) and Day 9 (Fig 9C and 9F and 9I). Fig 9J and 9K shows the prediction of the average centerline distance from model to data and from data to model, respectively. The average centerline distance is normalized by the length of the longest sprout in this region at day 3. This results in a normalized length of less than 10% from model to data and a normalized length of less than 20% from data to model. However, at day 9 the complexity in the vascular network, specifically the vascular remodeling that eliminates the anastomosis, is beyond the capabilities of the model to predict.
https://doi.org/10.1371/journal.pcbi.1009499.g009 Table 5. Error in the prediction of vessel centerlines. We compare the relative error in the prediction of vessel centerlines calculated from simulations of local vasculature (shown in Fig 9). C dist (d, m) is the average centerline distance between the skeletonized vessel in the data, d, and the skeletonized vessel predicted by the model, m.

Discussion
We have rigorously calibrated a hybrid, multiscale model of tumor angiogenesis, both with global and local vascular structure, with experimental data obtained from an in vitro vascularized tumor platform and tested the ability of the model to make accurate longitudinal predictions. This was done by first utilizing hemocytometry and ELISA to calibrate mathematical models of cell growth and VEGF concentration, respectively, yielding parameter distributions for VEGF production and consumption rates (Fig 3). The maximum likelihood values are used in subsequent calibration and prediction scenarios. Secondly, we analyzed global features such as sprout length, vascular density, and vascular volume fraction over time, to calibrate the global growth rate of angiogenic sprouts and the distance between new tip cells, shown in Figs 4 and 5. We then reported the errors and rigorously accounted for the uncertainty in model predictions, in Table 4 and Fig 6, respectively. The mean error in the prediction of the 1-and 4-parameter cases was 21% and did not increase with time. This is of note because the model predicts from day 0, not from previous time points (i.e., from day 5 to day 7), making prediction significantly more challenging (see Section 2.4.5). Finally, we analyzed the ability of the model to recapitulate local vascular structure by segmenting and skeletonizing specific vessels over time, initializing the model with confocal images at day 3, calibrating with images at day 5 and 7, and testing the predictability of the model against day 9 observations, as summarized by Fig 7. Our calibrated hybrid multiscale model, informed by hemocytometry and VEGF concentration data, can recapitulate local vascular structure longitudinally in our in vitro microfluidic platform. The calibration and prediction error, shown in Table 5, was on the order of 10% normalized centerline distance. The normalized length of the centerline distance from data to model increases across time and indicates that the hybrid ABM using the initial conditions of two tip cells (see Fig 9D), while quantitatively and qualitatively matching the description of vasculature at day 7, fails to capture the complexity observed at day 9. This is partially due to the restriction of new tip cells in the local calibrations, as we do not allow for new sprouts. However, the complexity across time points of the local sprouts dramatically increases and requires further study from both a modeling and experimental standpoint.
To the best of our knowledge, this work represents the first study to rigorously calibrate an agent-based model of tumor angiogenesis with multimodal experimental data to forecast vascular growth that is then directly compared to the corresponding data. Importantly, other efforts have pioneered the integration of data in other ways [65,66]. For example, Perfahl et al. employed multiphoton microscopy to image angiogenic vasculature in an in vivo dorsal skin fold chamber [7]. These images furnished the initial conditions of a cellular automaton model of vascular tumor growth, with blood flow, and vascular remodeling. However, no model calibration or comparison to experiments was performed. Similarly, in Xu et al., photoacoustic imaging of a murine xenograft model initialized a phase-field model of tumor angiogenesis [38]. While the authors did not utilize time-resolved vascular data, this experimental-computational approach showcases a first step toward longitudinally integrating in vivo data into angiogenesis models, through initialization, calibration, and ultimately, prediction. The present study incorporates initialization, calibration, and prediction through utilizing in vitro data; however, the experimental and computational approaches in this contribution can be made model-agnostic (and, even, data-agnostic). In particular, all the calibration methods were based on gradient free methods, and all the solvers are independent of the model, data, and the system under investigation.
Many of the calibrated parameter values are readily comparable to previous experimental studies and can be used in future modeling efforts. Specifically, discrete or hybrid mathematical models that include tip cell velocity, which is comparable to stalk cell divide time in our work, can be informed from the values reported in Table 6 [37][38][39]67]. However, in continuous models, specifically in the landmark Anderson 1998 model [5], chemotaxis during angiogenesis depends on the magnitude of the gradient of VEGF. As the magnitude of this value is unknown in our experiments, direct comparisons to the Anderson model (and any model with a purely continuous description of endothelial cells) is not feasible. While our calibrated cell doubling time (stalk cell divide time plus stalk cell growth time) is 2-3× less than what is observed in other experimental protocols, we use a simplified cell geometry with a cell diameter of 10 microns, 2-3× less than what has been observed in endothelial cells. This puts our doubling time within the range measured by others [68,69]. With regards to endothelial cell migration rates, endothelial cell velocity (in the context of random walks) has been measured to be 10-50 microns per hour in various conditions, while our model calculates an average sprout elongation rate (velocity) of 2.0 microns per hour. This difference is almost certainly because of the very different experimental scenarios considered in the literature which use a 2D framework and neglects the necessity of cells needing to degrade the extracellular matrix for migration. Our system is a 3D framework and explicitly accounts for degradation of the extracellular matrix [48], which necessarily results in a lower average velocity.
A key experimental limitation of the present study is the lack of longitudinal measures of various proteins (e.g., VEGF) and cell markers (e.g., to assign endothelial cells as tip, stalk, phalanx phenotype) in the in vitro microfluidic platform. We have addressed this by devising smaller calibration and prediction scenarios (with corresponding experiments) that can isolate the salient phenomena (e.g., VEGF production from tumor cells and growth of stalk cells) that we assume are fundamental properties of the cell and therefore independent of the experimental scenario. Another area for advancement is to determine the generalizability of both the local and global calibration results by applying the approach to multiple microfluidic platforms. We utilized one microfluidic platform over time, yielding two regions for global analysis and two local regions for calibration and prediction. Another fundamental limitation is the approximation that parameter distributions calibrated in simpler experimental/computational scenarios are valid for more complex scenarios. However, due to the complexity in both the 3D microfluidic platform and the tumor angiogenesis model, we are not able to isolate quantities of interest to inform individual model parameters in these complex scenarios. For example, we do not have additional temporal measurements of VEGF concentration to inform model parameters. We believe that by approximating these parameters using readily-available quantitative data from simpler scenarios (i.e., the 2D experiments) we can more accurately approximate the real parameter values in the 3D microfluidic scenario than by attempting to calibrate all the parameters using the limited 3D data. We note that in our global analysis in Scenario 4, as we utilize an ABM for angiogenesis, the model cannot be uniquely initialized from day 3 or day 5 and is thus restarted from day 0 with no angiogenic vasculature for each simulation. This is because the model is not continuous and cannot be populated from the vascular density at days 3 or 5, but the endothelial cells in the ABM must start from a parent vessel.
A key limitation in the modeling is that our calibration was only performed on 2D data. As out of plane effects could play a significant role in vessel sprouting, future efforts will be focused on extending the formalism to 3D. The main limitation to extending to 3D is the characterization of the angiogenic sprouts by endothelial cells that make up both sides (in 2D) of the sprout. This can be overcome by describing the vessels as nodes (which describe the radius of the vessel) connected by edges (which describe the directionality of the vessels), an approach taken by [70,71]. Another key modeling limitation is the use the maximum likelihood values from Scenarios 1 and 2, as opposed to Bayesian inference in Scenarios 3-5. If Bayesian inference is used in Scenarios 1 and 2, calculating the likelihood becomes a matter of sampling from the parameter distributions which increases the computational cost of the subsequent calibrations by the number of samples needed to characterize those distributions. In future work, we aim to utilize more experimental replicates to further calibrate and validate our model and to devise experimental protocols with anti-angiogenic and radiation therapies to model the effects of these treatments on angiogenic sprouts in vitro [72].

Conclusion
We have calibrated and determined the ability of an agent-based model to make accurate predictions of tumor angiogenesis by systematically incorporating data of different scales to inform model parameters. To the best of our knowledge, this represents the first effort to calibrate a mechanism-based mathematical model to spatially and temporally-resolved experimental data of angiogenesis, thereby enabling predictions of future vessel development that could then be directly tested against observation.
Supporting information S1 Appendix. Model invertibility, Gaussian fit of parameters, and time evolved solution of the ABM. In this appendix, we show the invertibility of the models by calibrating with in silico data. We then explore the accuracy of the Gaussian fits of the parameter distributions calibrated in Scenarios 1-5. The longitudinal prediction of the ABM is presented after calibration during Scenarios 1-4 and another region of local vasculature is presented. (DOCX) S1 Fig. Model invertibility of ODE models. Panels (A)-(C) depict the parameter distributions from calibrating VEGF production, consumption, and carry capacity, respectively, to 5% Gaussian noised VEGF concentration data. Panel (D) shows the goodness of fit of the VEGF concentration in the tumor (IBC3) and endothelial (TIME) experiments, shown in red and blue, respectively. The relative error between the parameters used to generate the data and the calibrated parameters are 9.3%, 2.0%, and 5.8%. Each column depicts the data (green) and model prediction (red), with overlap in yellow, of days 5, 7, and 9, respectively. Each row shows the voxels predicted by the centerlines using different thresholds, with the voxels in the 99% prediction (top row), 95% prediction (middle row), and 90% prediction (bottom row) of the simulations. In Panel (C), the direction of vessel growth in the model shows the VEGF gradient going away from the direction of data and toward the hypoxic cells shown in Fig 9F. This leads to the overall directionality of the vessel to be misaligned with the data. Panels (J) and (K) shows the prediction of average centerline value from model to data and from data to model with an averaged normalized length difference of less than 2% and 12%, respectively. The average centerline distance is normalized by the length of the longest sprout in this region at day 3. (TIF) S1