Sensitivity analysis of a multi‐physics model for the vascular microenvironment

The vascular microenvironment is the scale at which microvascular transport, interstitial tissue properties and cell metabolism interact. The vascular microenvironment has been widely studied by means of quantitative approaches, including multi‐physics mathematical models as it is a central system for the pathophysiology of many diseases, such as cancer. The microvascular architecture is a key factor for fluid balance and mass transfer in the vascular microenvironment, together with the physical parameters characterizing the vascular wall and the interstitial tissue. The scientific literature of this field has witnessed a long debate about which factor of this multifaceted system is the most relevant. The purpose of this work is to combine the interpretative power of an advanced multi‐physics model of the vascular microenvironment with state of the art and robust sensitivity analysis methods, in order to determine the factors that most significantly impact quantities of interest, related in particular to cancer treatment. We are particularly interested in comparing the factors related to the microvascular architecture with the ones affecting the physics of microvascular transport. Ultimately, this work will provide further insight into how the vascular microenvironment affects cancer therapies, such as chemotherapy, radiotherapy or immunotherapy.


| INTRODUCTION
Oxygen heterogeneity in solid tumors is recognized as a limiting factor for therapeutic efficacy. 1Such heterogeneity arises from the abnormal properties of the vascular microenvironment (from now on VME), but the precise mechanisms of compromised oxygen transport are only partially understood.
Multiscale and multiphysics models 2,3 (i.e., computer models that examine cancer behavior across different spatial, temporal, and/or functional biological scales) are uniquely positioned to capture the space-and time-dependent changes and the heterogeneities that occur in tumor properties, and to provide potentially clinically useful insights.In this context, the multiphysics model adopted here, previously developed by the authors and co-workers, [4][5][6] is a simplified representation of the VME, retaining two fundamental aspects: (i) the physical complexity of flow and oxygen transport; (ii) the complexity and variability of the microvascular network geometry.
8][9] Microvascular topology can differ for several reasons.For example, the dimension of the extravascular volume gaps in the tissue may vary significantly in different microvascular network samples.From the standpoint of vascular architecture, a tumor microenvironment can be related to the presence of a highly heterogeneous spatial distribution of capillaries; in particular, the risk of hypoxia is increased in the low-vascularized regions.In physiological conditions, the network is characterized by a more regular and dense vascularization.Therefore our computational study encompasses capillary structures with high topological variability.
The main purpose of this work is to compare these multi-factorial traits of the model.For rigorously quantifying the influence of the model parameters on selected quantities of interest (QoIs), we set up a sensitivity analysis (SA) approach.Sensitivity analysis 10,11 provides qualitative and quantitative instruments to determine how input variations affect the model responses, and thus rank the relative contribution of each parameter, either singly or in combination with the others.Sensitivities are usually discovered by running the underlying model against different configurations and analyzing the statistical properties of the associated input-output samples.
Sensitivity analysis and uncertainty quantification (UQ) are increasingly becoming required steps in the workflow for the reliable application of complex biological models in bio-engineering design and clinics.As a few examples among many others, we mention UQ studies for cardiovascular applications, 12,13 for brain perfusion models, 14,15 and tumor growth. 16,17n this work, we perform SA studies on the aforementioned physical and morphological inputs of a multiphysics model describing the VME.In particular, we rely on two well-known sensitivity measures, which are Sobol' indices and Morris elementary effects.The former are quantitative measures, which however are usually computationally expensive, whereas the latter provide a qualitative ranking criterion relying only on a small number of model simulations.These methods are applied to the problem at hand to assess the overall influence of the model parameters, as well as their nonlinear effects and interaction with other inputs.Furthermore, we compare these methods with each other, aiming to determine what is a more suitable choice as a SA technique when the computational model costs are significant. 18e are particularly interested in the relation of the VME with tumor treatment, in which oxygen plays an essential role.For this reason, the selected quantities of interest are mainly related to the oxygenation of the tissue.Both physical parameters and morphological ones may have an impact on tissue oxygenation and related cancer therapies.For example, high values of oxygen consumption rate may result in a low oxygen concentration in the tissue; in turn, we expect that the network topology influences the oxygen spatial distribution.Nonetheless, other model parameters may have little or even no impact on all the QoIs, and could be fixed to a reference value without affecting the reliability of predictions.Quantifying these dependencies in a rigorous way will facilitate interpreting the role of VME in cancer therapies, 8,19 such as chemotherapy, radiotherapy, or immunotherapy.

| A GEOMETRICAL MODEL OF THE MICROVASCULATURE AND A PHYSICAL MODEL OF ITS FUNCTION
In this section, we report the models describing the VME.We subdivide the relevant models for the VME into two categories: (i) the model of the vascular geometry; (ii) the physical models of the vascular function.For the former model, starting from the well-accepted approach based on Voronoi diagrams, we propose a new parametrization that is specifically tuned for our sensitivity analysis approach.On the other hand, we adopt the physical models previously addressed in the literature [4][5][6] (also by co-authors of this work and their collaborators), but we believe it is necessary to provide the reader with a basic knowledge of such models in order to understand the results discussed later on.
The physics based models consist of mixed-dimensional 3D-1D model of the VME describing blood flow and oxygen transfer from the microvasculature to the tissue, representing microvessels as one dimensional channels.As discussed in References 20,21 this approximation significantly simplifies the problem at the computational level.

| Geometrical model of the microvasculature
One of the main challenges to model microcirculation is to surrogate the capillary network organization.In particular, the goal is to produce artificial structures reproducing main topological and morphological properties of microvascular networks.Vascular architectures can vary considerably, depending on the spatial scale of observation and on the organ.For the smallest scale, the one of capillaries, Voronoi-based models represent efficient methods to generate 3D synthetic microvascular networks, as stated by Smith et al., 22 due to the space-filling pattern that is generated and the constraint of three connections on each junction node.In order to provide heterogeneity in vessels' spatial distribution, a pruning strategy similar to the one of Grogan et al. 8 is carried out, entailing a more significant vascular variability and reproducing some topological features observed in the tumor microenvironment.
Hinging upon the aforementioned approaches, we introduce specific metrics to measure the variability of a capillary network and to provide structures with peculiar characteristics.Vascular networks can be characterized according to different types of metrics, some are based on counting algorithms, and others derive from more complex analyses.We consider morphometrical metrics (such as the distribution of vessel length, the number of vessels per volume and the total length of vessels per volume, and the interior and boundary vertex density).In particular, we use two metrics for analyzing the vessels structure: • S/V: the lateral surface (S) of a capillary vessel over the tissue volume (V ); • d max : the maximum extravascular distance (EVD) from the tissue domain to the 1D network.
The second is a local indicator that measures the space-filling nature of capillary networks and it has been used to discriminate between healthy and tumor networks. 23

| Microvascular flow and hematocrit transport
We model the fluid flow within the vasculature and the surrounding tissue through the 3D-1D problem (ℱ).In general, such model consists of two main ingredients, the governing principles (essentially mass conservation) and the constitutive laws for flow in the vasculature and in the interstitial space.Concerning the latter, the interstitium is modeled as a saturated porous media, employing Darcy's equation.We use Poiseuille equation suitably modified to account for nonlinear blood rheology for describing the blood flow within the vasculature.The effect of the red blood cells at small vascular scales is a complex problem that requires complementing the balance laws of fluid mechanics with several closure models. 24We account for red blood cells through the hematocrit model (ℋ) that describes their influence on viscosity, by taking into account the Fahraeus-Lindqvist effect, 25 and their nonuniform split when traveling through a bifurcation, namely phase separation, the plasma skimming or Zweifach-Fung effect.For the latter effect, many closure models are available 25,26 and can be used in the forthcoming computations.For modeling consistency among Fahraeus-Lindqvist and Zweifach-Fung effects, with use here the formulas of Pries and Secomb. 25Fluid filtration through the microvascular wall is described by the classical Starling equation, accounting for hydraulic and osmotic pressures.Moreover, lymphatic drainage is included as a nonlinear function of the interstitial pressure. 4The problem is complemented by the mass conservation equations in the two domains (the tissue Ω and the vasculature Λ), and it reads as follow: Þin Λ, Starling 0 s eq: where subscripts t and v stand for the tissue and the vasculature respectively, precisely u t , u v are the fluid velocities, p t ,p v are the fluid pressures, μ is the viscosity, K is the porous media permeability, R is the radius of the vessel, ϕ V is the fluid extravasation, ϕ L is the lymphatic drainage, H is the discharge hematocrit, L p is the hydraulic conductivity of the microvascular wall and π v À π t is transmural the osmotic pressure jump.In the context of the 3D-1D formulation of blood flow, δ Λ is the Dirac delta function located on the centerline of the vascular channels, named Λ.
Equation (1) 6 is not sufficient to uniquely determine the variation of hematocrit in the network.It must be combined with suitable conditions for the conservation of hematocrit at the junctions and at the boundary of the network.As (1) 6 is a pure transport equation, it is well known that we have to prescribe a constraint on hematocrit at each inflow point of the network branches.For the internal junctions, we exploit the mass conservation of hematocrit.Let us consider a generic junction with multiple branches joining at a single node.Given the orientation of the flow in each branch, we subdivide the branches into K out ¼ card K out j outflow ones and K in ¼ card K in j inflow branches.We prescribe as many constraints as the number of inflow branches, namely K in .Mass conservation always provides one constraint that is, The previous equation is not sufficient to close the problem in the case K in > 1.The simple case K in ¼ 1 identifies anastomoses, where one, two, or multiple outflow branches merge into a single inflow one.In this case, since all the terms on the left-hand side are known, the hematocrit value on the right is uniquely determined.In the case of bifurcations, namely K in ¼ 2, the problem can be solved using a flow split model.In particular, we adopt here the one proposed in Reference 25.Since we exclude the presence of trifurcations or more complex configurations, this approach will be entirely sufficient to determine the distribution of hematocrit in the network.Without loss of generality, let us consider the classic Y-shaped configuration, where one channel divides into two branches.We denote by the subscript f the quantities related to the parent channel and with α, β the daughter branches.Given the blood flow rates with Ã ¼ f , α, β and the outflow hematocrit H f , we aim to determine H α and H β , which provide hematocrit values at the inflow of the bifurcation branches.In practice, we define, and we calculate these fractions by means of the following model where, A, B are fixed parameters determined in Reference 27, logit and X 0 is the fractional blood flow rate under which any RBC will flow into the daughter branch α.Finally, the desired hematocrit levels are determined as We note that the microvascular flow model presented here is affected by a limitation about the description of Starling forces.These forces affect the fluid exchange between the vascular network and the tissue, are determined by the osmotic pressure jump, π v À π t in Equation (1) 5 , which are assumed to be fixed.Because of the interaction with the lymphatic drainage, ϕ L in (1) 2 , Starling forces may be locally perturbed from their equilibrium state.From the modeling standpoint, this effect must be, in future, added to the model, for example, as in Linninger et al. 28,29 From the computational standpoint, this improvement entails to complement the model with transport equations for large proteins, for example, albumin, similar to the hematocrit model.Further details regarding the modeling equations and the numerical solution of the problem can be found in Reference 4. The problem is then complemented by defining values for the vascular pressure jump between the endpoints of the vascular network.This determines the main orientation of the vascular flow and the separation of the vascular endpoints into inlets and outlets.The value of hematocrit is then prescribed at the inlet points of the network.For tissue, mixed conditions representing resistance to flow are enforced to account for the environment surrounding the domain of interest.

| Oxygen transport
Oxygen is delivered to tissue, where is depleted by the cellular uptake.The transport equations for the two domains comprise diffusion, consumption, and advection, leveraging data from the models (ℱ) and (ℋ).Besides the above-mentioned component, we introduce three particular phenomena: (i) the oxygen binding to hemoglobin, (ii) the oxygen uptake dynamics, and (iii) the oxygen exchanges through the microvascular walls.
For the first item, we recall that oxygen is present in the system as solved in plasma/interstitial fluid (C v or C t ) or bound to hemoglobin (C HbO 2 ).The total content of oxygen in the vasculature is the sum of the two.We assume the binding dynamics as an instantaneous process, described by the Hill's equation: where, k 1 is a constant value equal to N MCHC (being N the Hüfner factor, which represents the oxygen binding capacity of human hemoglobin and MCHC is the Mean Corpuscular Hematocrit Concentration), α pl is the oxygen solubility in plasma, p s 50 is the oxygen partial pressure at hemoglobin half-saturation, and γ is the Hill exponent.
The second item addresses the cell uptake in the tissue by the well-known Michaelis-Menten equation 30 : where, V max is the maximum oxygen consumption rate, α t is the oxygen solubility in the tissue, and p m 50 is the Michaelis-Menten constant, that is, the oxygen partial pressure at half consumption rate.The third item regards oxygen exchanges with the tissue, governed by diffusion and advection through the microvascular wall, modeled as a semi-permeable membrane 31 : where, Υ O 2 is the vascular wall permeability.We point out that for small molecules the diffusion across the endothelial layer largely dominates over transport.For example, using the parameters reported in Table 2 for oxygen, we see that the coefficient 2πR Υ O 2 is four orders of magnitude larger than ϕ V .The second term on the right-hand side of (4) is useful in case the semi-permeable model is used for the transfer of large molecules, such as proteins.
Comprehensively, the problem reads as follow: The model describes the oxygen transport through the microvasculature and its delivery to the tissue (ϕ O 2 ), where it reaches the cells and cell uptake occurs.The problem is complemented with boundary conditions defining oxygen concentrations.The network inlets and outlets are determined according to the orientation of blood flow, in turn, determined by the enforcement of the vascular pressure jump between the endpoints of the network.At the network inlets Dirichlet type boundary conditions are prescribed.At network outlets, we set null axial derivative of the oxygen concentration (i.e., null diffusive flux).Finally, we apply null concentration gradients at the artificial interfaces with the tissue.We note that these boundaries are anyway open to oxygen convection by means of the transport field u t .

| Quantities of interest
When dealing with models based on parametric PDEs, it is crucial to understand how variations or uncertainties in the inputs affect selected outputs of interest.These are quantities associated with the solution of the problem that can be useful from the application viewpoint.For the problem under consideration, the focus is on the oxygen transport model, described by (O), that provides the oxygen concentration C t : Ω !ℝ þ as one of the state variables.Then, the spatial average C t , yielding a global description of the tissue oxygenation, is evaluated as a QoI and computed as follows: Furthermore, the spatial variability of the oxygen concentration is considered as QoI, and is identified by the difference between the maximum and the minimum concentration ΔC t : with the goal of quantifying possible local hypoxic effects in the tissue.The other state variable describing the oxygen transport is the concentration of dissolved (free) oxygen in the vessels C v : Λ !ℝ þ , defined for each node in the 1D mesh of the vascular network Λ. Proceeding analogously to the oxygen in the tissue, we define as QoIs the spatial average C v : and the spatial variability of the oxygen concentration in vessels measured by

| Comparison of the model with real data
Validation (in a broad sense) is a fundamental step in the application of any theoretical model.Although a thorough validation of the proposed model is not in the scope of this work, it is mandatory to provide evidence that the predicted quantities fall in a physiological range.To this purpose, we guide the reader through previous works where the results of the computational model proposed here (combining the mathematical equations and the numerical solver) have been compared with real data.In Reference 40, the authors and co-workers compare the outcome of model ( 1) (obtained using a tissue slab of 540 Â 740 Â 400 μm in which is embedded a network of 28 branches of diameters ranging from 56.4 to 18.8 μm, similar but not equivalent to the one analyzed here), with the data of microvascular flow for the category of arterioles provided in Table 1. 41Concerning oxygen transport, from a general standpoint, our results in T A B L E 1 Mesh dependence study for the interstitial domain Ω, evaluating the QoIs C t and ΔC t .

Tetrahedra
Total nodes CPU time (min) terms of tissue oxygen partial pressure also agree with values reported in the literature for the brain tissue.In our tests, the average oxygen partial pressure in the tissue is about 60 mmHg, with a spatial variation of approximately 35 mmHg (see Figure 12).These data are in agreement with values reported in the literature for brain tissue. 42We also remark in the discussion that the range of variation of hematocrit (H in , see Table 2) ensures in any case good blood oxygenation.As a result, our simulations fit well into a microvascular environment with good oxygenation.

| SENSITIVITY ANALYSIS OF VASCULAR NETWORKS
In this section, we provide an overview of the techniques used in this work to perform global SA, namely variance based methods (see Section 3.1) and screening methods (see Section 3.2).These are two families of approaches for SA with different characteristics.The former enjoys a more rigorous theoretical framework and the computed indicators convey more information about the parameters, with the price of a large amount of samples required to form appropriate statistics.The latter has more heuristic foundations and correspondingly provides weaker results, leading to a significantly lower computational cost.From the methodological standpoint, the purpose of this work is to compare these strategies and establish what is a good tradeoff between reliability of the SA results and the computational cost to obtain them, in the particular case of multi-physics models of the microvascular environment.
T A B L E 2 Input parameters used for the numerical tests and ranges of variation of the physical and geometrical input factors for the SA.
Oxygen partial pressure at half consumption rate mmHg --0:5 À 1 Percentage of negative seeds over positive Voronoi seeds

| Variance-based methods: Sobol's indices
Variance-based methods represent a powerful tool to quantify the relative importance of individual factors or groups and are derived from a suitable high-dimensional model representation. 43,44The theoretical concepts introduced are mainly based on. 11,45Let us consider the probability space Ξ,A, P ð Þand the random variables , where from the modeling perspective each random variable corresponds to one of the input parameters of the oxygen transport model.We assume all the parameters to be independent and uniformly distributed, and normalized in the range 0, 1 where U denotes a uniform probability distribution.From a more abstract perspective, a generic model can thus be seen as a scalar-valued, nonlinear function of the random input vector , where the model realization y x ð Þ corresponds to a selected output of interest associated with the PDE solution computed for x ¼ x 1 , …,x k ð Þas an input parameter.A primary variance-based measure of sensitivity is given by the first Sobol' index which measures the effect of varying X i alone averaged over variations in all input parameters.The first Sobol' index S i represents the main contribution of the input factor X i to the variance of the output: the higher is S i , the greater is the influence of X i on the output Y .In fact, according to the law of total variance, we can write where the left-hand side corresponds to the model variance once a source of variation (i.e., X i ) has been fixed.Thus, the more influential is X i , the smaller is , hence the greater must be S i .Note that P i S i ≤ 1.Nonetheless, for non-additive models, it is not possible to separate the effects of its inputs on the output and one should look for higher-order interactions.Under few assumptions, the variance of the model output can be decomposed as where, V ij , V ijl are the contributions to the variance of the interactions between X i , X j , X i ,X j , X l , etc. up to the term V 123…k considering the interactions between all variables, showing how the variance of Y can be decomposed into terms related to each input and to the interactions between them.To avoid the calculation of all the interaction terms, that would require the evaluation of 2 k À 1 indices, thus becoming computationally demanding (e.g., in what follows we will use k ¼ 7 ending up with 127 possible interaction terms), we can rely on the total effect indices Þdenotes the random vector of all input factors but X i , and corresponds to the (average) variance of Y left if each factor in X $i could be fixed to its true value, so that the smaller is S T i , the less influential is X i and can be arbitrarily fixed within its range of uncertainty without appreciably affecting the output of interest.Unlike S i , it holds P i S T i ≥ 1 since S T i 1 and S T i 2 , for i 1 < i 2 , take both into account the interactions between X i 1 and X i 2 .
In conclusion, the first Sobol' index and the total effect indices are rigorous mathematical indicators that quantify the variance of outputs under the assumption of a given variability of inputs.As addressed in the next section, the computation of these indicators entails a significant computational cost required by the approximation by means of quadrature formulas of multi-dimensional integrals, appearing in the definition of expected value and variance, where each function evaluation correspond to a solution of the physical model For most models, the sensitivity indices cannot be calculated analytically and must be estimated by suitable techniques.One of the most used and efficient procedure is the Saltelli method 10 (based on the approaches originally proposed in Reference 43), which avoids the cumbersome computation of multidimensional integrals by brute-force.
Let N S ℕ be a prescribed integer, known as base sample.The Saltelli method consists of the following steps: 1. generate a N S Â 2k matrix of random parameter realizations obtained, for example, from a Sobol' quasi-random sequence, 2. define two matrices of data A,B ℝ N S Âk each containing half of the samples, that is, where, x i is the realization of the ith random variable X i .

construct k matrices
formed by all columns of B except the ith column taken from A, that is, 4. compute the output of the model for all the parameter vectors given by the rows of A,B and C i , for i ¼ 1, …, k, thus obtaining the N S -dimensional vectors, . for i ¼ 1, …, k, estimate the ith first-order sensitivity index as follows: where, y A ¼ 1 6. finally, for i ¼ 1, …,k, estimate the ith total order index as follows: As a result, computing a full set of S i 's and S T i 's indices for k input factors requires S runs that would have been necessary for the brute-force approach.Nonetheless, the accuracy of the estimates heavily depends on N S , which can vary from a few hundred to a few thousand, so that, the main drawback of this method is still its computational cost.

| Screening methods: elementary effects
For computationally demanding models, a more efficient alternative to variance-based methods is given by screening methods, such as the Elementary Effects, proposed by Morris in Reference 46.This method allows us to identify which input factors are (i) negligible, (ii) linear and additive, or (iii) nonlinear or involved in interactions with other factors, using a relatively small number of model evaluations.By considering wide ranges of variations for the inputs and averaging over a number of local derivative approximations, this method overcomes the main limitations of one-at-a-time designs and derivative-based methods.
Given a partition Ω ℓ of the parameter domain 0, 1 ½ k into ℓ ℕ steps, the elementary effect associated with the ith input factor and computed at x Ω ℓ , for i ¼ 1, …,k, is defined as where, x þ e i Δ Ω ℓ , being e i the unit vector of all zeros but the ith component, which is equal to one, and g a user-defined grid step size in the parameter space, which is kept constant for all parameters.Since EE i x ð Þ, for i ¼ 1, …, k, quantifies a local behavior, as it depends on x, to obtain a more global sensitivity measure one needs to estimate statistics of their distribution, that is, F i $ EE i , constructed by randomly sampling different x from Ω ℓ .Moreover, to avoid cancelation effects that can occur when F i contains both positive and negative elements, the mean of the distribution of the absolute values of the elementary effects, that is, G i $j EE i j, is also computed.Thus, for r ℕ sample points x 1 ð Þ , …, x r ð Þ , we obtain the following sensitivity measures for X i , Note that, in order to guarantee an equal probability sampling for each F i and G i , i ¼ 1, …, k, it is convenient to choose the step value Δ in the middle of the admissible interval, that is Δ ¼ ℓ=2 ℓ À 1 ð Þ, being ℓ an even number.
The mean μ i (or μ Ã i ) represents the magnitude of the ith effect and can be used to quantify the individual influence of X i on the output.We remark that μ Ã i represents a valid substitute for the total index S T i previously defined.We also observe that the comparison between μ i and μ Ã i provides information on the signs of the effects of X i on the output.The variance σ i measures the spread of EE i and thus estimates the effect of the ith input due to nonlinearities and interactions with the other factors.In particular, low values for σ i can be interpreted as the fact that the effect of X i on the output is independent of x, that is, of the values taken by the other inputs, whereas a large variance is a symptom of strong interactions between inputs or nonlinear effects.
A qualitative classification of Morris indices can be done using the σ versus μ Ã plane, see Figure 1.Here, the measure μ Ã i provides information about the importance of the ith input parameter and the ratio σ i =μ Ã i is an indicator of linear dependence.Using these indicators, we obtain a classification of the inputs 47 from the scatterplot of the sensitivity indices EE i , 46 for i ¼ 1,…, k.More precisely, relying on statistical properties, three thresholds for σ i =μ Ã i have been identified.A smaller value of the σ i =μ Ã i distinguishes factors with a more linear and monotonic behavior.In particular, the case σ i =μ Ã i ¼ 0 is associated to a true linear response.Conversely, a large ratio detects inputs with nonlinear effects and interactions with the other parameters.

| Computation of the sensitivity indices: Morris sampling
Since the computation of each elementary effect requires two sample points, a naive design would require 2rk model evaluations to construct μ i , μ Ã i and σ i , for i ¼ 1, …, k.However, Morris 46 has introduced a special sampling strategy that allows us to rely on only r k þ 1 ð Þ samples by employing r trajectories of k þ 1 points, each differing from the neighbor in only one component, thus providing k elementary effects per trajectory.
Each trajectory can be seen as a k þ 1 ð ÞÂk sampling matrix B Ã built such as, for every j ¼ 1, …, k, there are two rows that differ only in the jth entry.To construct B Ã , the most straightforward strategy would be to choose a random starting sample x Ã Ω ℓ , and define where, J mn is a m Â n matrix whose elements are all equal to 1, for some m,n ℕ, Δ is the step value and B is a k þ 1 ð ÞÂ k strictly lower triangular matrix of ones.In this way, the orientation matrix B Ã has the form where, However, the elementary effects associated with this matrix would not be random.To overcome this drawback a k-dimensional diagonal matrix D Ã of randomly chosen entries in À1, 1 f gand a k Â k permutation matrix P Ã are introduced, so that the sampling matrix is finally defined as F I G U R E 1 Typical parameter classification derived from visual inspection of the elementary effect statistics.The ratio σ i =μ Ã i provides information regarding linearity and monotonicity for the ith input factor.
Here P Ã is used to permute the order of the directions in the trajectory, whilst D Ã takes into account the fact that B iþ1,j , that is the value of the jth entry of the i þ 1 ð Þ-th point in B changing with respect to the previous point, for some i ¼ 1, …, k and j ¼ 1,…, k, can be either increased by Δ, that is, B iþ1,j ¼ B i,j þ Δ, or decreased, that is, B iþ1,j ¼ B i,j À Δ (or, more specifically, B i,j ¼ B iþ1,j þ Δ).Finally, one elementary effect per input can be computed by subtracting the model evaluation at two consecutive rows, divided by Δ.For further details, we refer to Reference 11.

| Morris approach for group of parameters
When interested in ranking groups of inputs the Morris method can be suitably extended, based on the idea to move all factors of the same group simultaneously.In fact, using the sensitivity measure μ Ã alone allows us to take into account the case in which two or more factors have been changed in opposite directions.With this aim, we define the absolute elementary effect computed in x as where e x denotes the point obtained from x by increasing or decreasing each entry of the ith group by Δ.
To compute the estimated mean the sampling strategy described in Section 3.2.1 has to be slightly modified.First of all, one needs to define the k Â g matrix G, being g the number of different groups, whose element G ij is equal to 1 if the ith factor belongs to the jth group and 0 otherwise.Using a g-dimensional permutation matrix P Ã , the matrix B Ã of trajectories has now dimension g þ 1 ð ÞÂk and is defined as Note that in this case, since the order of the directions in the trajectory has to be rearranged in a way that is coherent with the group subdivision, only G is multiplied by the permutation matrix P Ã .

| GENERATION OF IN-SILICO MICROVASCULAR NETWORKS
The generation of synthetic vascular networks that preserve the morphometric properties of real ones is a challenging task.Some approaches among the most advanced 42,48 separate the problem into two steps: first arterial and venous network trees are defined and then suitable microvascular closure techniques are applied to seal the previously defined parts.Here we pursue an objective different from the one of reproducing realistic anatomical models.To perform sensitivity analysis we need to generate networks that exactly match some desired metrics.It is reasonable that, in this context, some additional simplification can be accepted.To tackle the task of reproducing the patterns of microvasculature at the smallest scale, essentially the one that is previously addressed by the microvascular closure, we chose one of the simplest and also widely used approaches for microvascular network generation (or closure), based on Voronoi tessellation. 49,50This method provides planar networks that are subsequently perturbed in three-dimensions.Furthermore, the chosen approach only provides straight connections between nodes (junction points) of the vascular graph.Curved small vessels or more importantly vessels with high tortuosity are not encompassed in our vascular models.
We generated vascular networks with prescribed values of the chosen metrics, that are, the capillary surface over the tissue volume S=V and the maximal extravascular distance d max .The generation of artificial networks with desired d max could be achieved by means of an inefficient trial-error strategy.Nonetheless, a proper sampling of each indicator is not granted, since the metrics cannot be sampled independently one from the other.To overcome these limitations we propose a novel indicator of the vascular architecture that is independent from S=V and is highly correlated to d max .Therefore, SA can be performed by sampling two independent metrics, one represented by S=V and the other one strictly depending on d max .

| Generation of (quasi) 3D vascular networks satisfying independent metrics
To fulfill the task of generating reasonable capillary networks for a sensitivity analysis study hinging upon constrained Voronoi-mesh based diagrams we start from Grogan et al., 8 where a similar objective has been previously addressed.In Reference 8 a Voronoi tessellation of a 2D domain is generated starting from a uniform spatial distribution of point seeds.The resulting tessellation emulates a common physiological capillary network, characterized by vessels with uniform spatial distribution.This procedure prevents degenerated links, only anastomoses or bifurcations are admissible.To perturb homogeneous distribution of vessels, Grogan et al. 8 use a pruning strategy where points are randomly selected from a normal distribution with the mean centered at the domain midpoint.Then the closest vessels to these points are progressively removed until a target density is reached.We pursue here the purpose of generating vascular networks with nonhomogeneous spatial distribution of vessels that satisfy prescribed constraints on the metrics S=V and d max .For this reason, we slightly modify the procedure presented by Grogan et al.The main difference of our approach to Grogan et al. is that we first generate a nonuniform distribution of seed points with the desired properties to match the metrics, then we build the Voronoi tessellation on it, as illustrated in Figure 2.
The starting point is to define an initial set of seed points p i f g (positive seeds) in the domain S ¼ 0, 1 ½ Â 0, 1 ½ , whose coordinates are obtained from a uniform distribution of pseudo-random numbers, that is, Then, to increase the disorder level we introduce a set of SEEDS À ð Þ points q i (negative seeds), whose polar coordinates r i ,

and center
x q 0 , y q 0 ¼ 0:5,0:5 ð Þusing the following probability distributions (see Figure 2A): For the application to vascular networks, we fix λ ¼ 16.Polar coordinates are then converted to Cartesian ones and the associated points subsequently undergo a translation by a vector of coordinates x q 0 , y q 0 .For each of these points, the nearest neighbor from the group of positive seeds is found and removed.This process is repeated sequentially, discarding as many positive seeds as negative ones.Afterward, the center of the disk q 0 is added to the remaining set of points, from which the constrained-Voronoi mesh is constructed (Figure 2B).
The distribution of the negative seeds based on the Poisson distribution is motivated by the need of generating networks with variable maximum extravascular distance and constant vascular density.The negative seeds distributed in this way weaken the vasculature density in a circular region at the center of the domain.By controlling the radius and the intensity of negative seeds in this region, it is possible to control the maximum extravascular distance and the vascular density independently.
The result of this procedure is a Voronoi tessellation with higher vascularization at the boundaries of S, consisting of a vascular void in the middle area (Figure 2C).Given a fixed set of positive seeds, the extravascular gap distance increases along with the ratio between the number of negative and positive seeds.
As a next step, a suitable perturbation of the Voronoi tessellation in the out-of-plane direction provides a 3D structure.We assume that the network is immersed in a 3D slab with a fixed height h ¼ 0:15, in dimensionless coordinates (150 μm).Therefore, an additional constraint is enforced to the z-coordinates of each Voronoi-mesh point: z i 0, h ½ , 8i ¼ 1, …, N, with N representing the number of nodes.The considered process consists of three phases: 1.For each mesh node, a z-coordinate is prescribed through a uniform probability distribution in the support 0:05h,0:95h ½ , in order to fulfill the spatial constraint in the vertical direction.
2. Each perturbed internal junction node is then shifted again vertically by a distance determined from a uniform distribution with support in À1:3h,1:3h ½ .3. The overflowing branches are cut off at the slab boundaries in order to enforce the boundedness by the domain depth on each mesh point.
By applying the described pipeline, we construct an in-silico microvascular network with desired vascular density and controlled level of pruning.An example of a 3D regular Voronoi network and an irregular one, derived from the proposed algorithm, are displayed in Figure 3.
In particular, fixing the ratio % , we finally obtain a synthetic indicator that controls accurately the dimension of the largest vascular void and, consequently, can surrogate the maximum extravascular distance for a sensitivity analysis.To support this claim, we have performed a correlation analysis between % We generated a large number of networks, which is associated to a different value of % Then, we evaluated the d max for each network to estimate the Pearson correlation coefficient (PCC) between the two indicators.The scatterplot reported in Figure 4 and the value PCC ¼ 0:9731 ≈ 1 show high positive linear correlation between d max and the ratio of negative and positive seeds.This entails that it is possible to perform SA using % with a proper sampling of the geometrical input parameters.

| Criteria for admissible vascular networks applicable to the model
In terms of spatial distribution, the network generation algorithm allows to determine networks with a desired variability.To model irregular and disordered networks, large vascular voids are produced in the middle area of the domain of interest.The volume gap diameter can be up to $ 500 μm if at most the 75% of positive seeds is removed.
The vascular generation algorithm requires to be complemented by imposing additional constraints on specific features of the vascular network.In particular, it is crucial to prescribe suitable vessels radii that have a relevant impact on the fluid and hematocrit mass balance at the junctions.We also impose a control over the number of inlets and outlets for each network.

| Prescription of vessels radii
Vessel radii are assigned in a physiological range of 2 À 6 μm , 51 with an average radius of the capillary vessels equal to 4 μm.To avoid vascular configurations that may become hardly solvable, we limit the variability to the range 2:5 À 5:5 μm.
In the first step of the pipeline for the prescription of vascular radii, the inflow and outflow nodes are assigned to the vascular graph, influencing the average blood flow direction in the 3D-1D multiscale model.Precisely, we assume Subsequently, the graph nodes are visited, starting from the inlets and moving in the direction of the flow.We assign the radii for the inspected branches depending on the network connectivity.In particular, if a parent vessel bifurcates, the daughter branches radii are a fraction determined from a uniform distribution in the range 95% À 99% of the inlet branch radius.Otherwise, if two parent vessels connect, the outgoing vessel radius is selected from the uniform distribution in the range 101% À 105% of the largest inlet vessel radius.
If a computed radius exceeds the aforementioned interval limits, the respective minimum or maximum of the admissible range is enforced, avoiding vascular networks with too large or too narrow channels.

| Inlets and outlets of the vascular network
A relevant factor that can significantly influence the computational model under consideration is the number of vascular inlets ∂Λ IN and outlets ∂Λ OUT .Indeed, by fixing the geometry vascularization and the dimension of the largest showing high positive linear correlation between the two markers, as also stated by the value PCC ¼ 0:9731.vascular gap, the network generation algorithm could provide structures with high clustering of vessels in the inflow regions, leading to increased oxygenation in the peripheral area.Analogously, a richer set of outflow regions facilitates a decrease of oxygen concentration for both the capillaries and the tissue.Since the aim of this work is to assess the effect of the network structure itself on the oxygen transfer process, we modify the vascular network generation algorithm by imposing a control over the sets ∂Λ IN and ∂Λ OUT , depending on the vascularization and the dimension of the largest void area.In particular, defined we have performed two linear regressions for the responses N inlets and N outlets , assuming S=V and % as predictors, based on 150 supervised vascular networks where the parameters N inlets and N outlets fell in the admissible range.The results of this regression are reported in Figure 5 and in the corresponding formulas.Provided the linear regression models, for a given desired sample of S=V and of % we estimate the expected values of N inlets and N outlets , generate a new vascular network, and accept it or reject according to the correspondence of the effective number of inlets and outlets with the predicted value, AE20 additional inlets/outlets.

| NUMERICAL RESULTS
After a brief discussion of the model setup, in this section, we present the results of the sensitivity analysis performed using the methodologies of Sections 3.1 and 3.2.Once the most influential inputs on the quantities of interest have been determined, we briefly show their effect on the oxygen distribution in the VME.

| Numerical simulations setup
We employ the finite element method to discretize the microvascular model (ℱ þ ℋ þ O), and use a fixed point method to linearize the equations containing nonlinear factors.We initialize the microvascular environment with the parameters presented in Table 2. Full details related to model derivation, discretization, solution and validation are reported in Reference 6.The model is applied to describe a 3D domain Ω, represented by a slab of edge 1 mm and thickness 0:15 mm and discretized through a structured mesh of tetrahedrons.The uniform grid consists of 7200 tetrahedrons, for a total of 1764 number of nodes.This setup is supported by a mesh dependence study.More precisely, we have evaluated the relative error between the spatial average of the oxygen concentration in the tissue for the current mesh and the numerical solution C t Ã computed after mesh refinement.A similar analysis has been done for the total spatial variation ΔC t with respect to ΔC Ã t .We report in Table 1 the values obtained considering a vascular network with vascular density given by S=V ¼ 6000 and % The difference between the coarser and the finer mesh is almost negligible in terms of the QoI C t , and reasonably low in terms of ΔC t .Conversely, the computational time required to perform a simulation triples when considering the finer mesh, instead of the coarser one.In order to further highlight F I G U R E 5 Linear regression models fitting the number of inflow (left) and outflow vessel nodes (right) with respect to the vascularization (S=V ) and the dimension of the largest void area % the negligible differences in accuracy between the meshes, Figure 6 displays the distribution of C t along a line ϕ running along the diagonal of the tissue slice (associated to z ¼ 0).The interstitial oxygen concentration is plotted for each 3D mesh refinement with respect to a normalized coordinate s that runs along the line.The profiles turn out to be almost overlapping.Therefore, since there is not a significant increase in accuracy, we decided to rely on the coarser 3D mesh for our study, which has a cell width comparable to one used in previous studies (Δx of 50 μm 8 and 40 μm 52 ).
Regarding the spatial discretization of the 1D domain Λ associated to the vascular network, we fix the number of subintervals for each branch to N subintervals ¼ 5 equispaced segments.Indeed, the considered geometries are highly vascularized, entailing a large number of elements for the 1D mesh increasing N subintervals .The discrete problem is assembled and solved by an in-house code 32 developed using the GetFem++ library. 53The numerical tests were performed on clustered computational resources, consisting of a processor AMD EPYC 7301 16-Core Processor with 2 sockets and 16 cores.

| Sensitivity analysis setup
For the SA studies, we consider k ¼ 7 input parameters, 5 of which are physics-related and 2 linked to the metrics of the vascular architecture.The input ranges of variation are reported in Table 2, spanning from physiological to pathological values.In these ranges, we assume uniform probability distributions for each parameter.For what concerns the vascular networks, once provided the values of synthetic indicators S=V and % , each artificial network is generated through an iterative process in which the number of seeds varies within a fixed range, in order to obtain a network whose vascular density is such that the difference from the given target S=V is under a fixed uncertainty threshold (1 m À1 , corresponding to the 0:02% with respect to the reference value of S=V 6000 m À1 ).
For the computation of the elementary effects, we point out that they exploit normalized QoIs derived from the oxygen transfer model.Once each snapshot of the solution of (O) has been computed using dimensional physical quantities, we rescale the computed oxygen maps on the characteristic value of oxygen concentration C, specified in Table 2, in order to retrieve the model realizations Y ¼ y X ð Þ, so that the sensitivity indices provide a relative measure on the impact of the input factor on the QoIs.In particular, the nonlinear function of the input sample vector X encloses a normalization with respect to C before performing the evaluation of each QoI.
The sensitivity indices have been computed in MATLAB.The source code for the computation of Morris elementary effects is based on the MATLAB toolbox SAFE, 54 further extended to handle groups of parameters, whereas that for Sobol' main and total indices is based on the Dakota toolbox, 55 originally written in C++.
First of all, we perform a convergence analysis using Morris method to retrieve an optimal number of trajectories r Ã , in the sense that the mean μ Ã i and the standard deviation σ i of the ith elementary effect, for all i ¼ 1, …, k, are unaffected when considering a higher number r > r Ã of trajectories.
The results reported in Figure 7 emphasize a different overall behavior of the indices μ Ã i and σ i when calculated for the physical parameters or the geometrical inputs: while for the former a plateau is reached just with a small number of trajectories, regarding the geometrical inputs, a slower convergence of both μ Ã i and σ i is observed.This is reflected, in particular, by higher values of the standard deviation σ i .Overall, the fluctuations of μ Ã i and σ i detected for small values of r are mitigated by increasing the number of trajectories.Nonetheless, the overall ranking of the inputs is preserved.Hence, from now on, we assume r Ã ¼ 60 to be a sufficient number of trajectories to compute accurate Morris indices.

| Sensitivity analysis varying all parameters independently
The analysis of the elementary effects of all the parameters computed using the Morris sampling (without making any distinction between the geometric parameters and the physical ones) is reported in Figure for the QoIs related to the tissue and the vascular oxygen concentrations.The plots show the influence of each parameter in terms of the mean μ Ã and the standard deviation σ of the elementary effects following the template described in Figure 1.
According to these results, which will be discussed in detail later on, we highlight some general trends.In all the considered scenarios the elementary effects associated with H in and Pm 50 have low μ Ã and low standard deviation σ, thus implying that the associated inputs have negligible effects on the chosen outputs.For this reason, we deflate the parameter space from k ¼ 7 to k ¼ 5 variables, in all the forthcoming analyses.Conversely, we observe that the geomet- and S=V are not negligible for all the QoIs.
In Table 3, we plot the inputs ranked with respect to the mean μ Ã i , for i ¼ 1, …, k, of the elementary effects.For all QoIs, we subdivide the range of μ Ã i , with i ¼ 1, …,k, into four equally spaced sub-intervals (improperly named quartiles) and rank each input variable according to the quartile where it falls.We notice some differences with respect to the ranking between QoIs based on the mean value (C t , C v ) and the total variation (ΔC t ,ΔC v ).The physical parameters have a greater influence on the mean QoIs, while the geometrical inputs have a greater impact on the total variation of vascular and tissue oxygenation.Furthermore, exploiting the classification of the inputs through the ratio σ=μ Ã F I G U R E 7 Statistics of the elementary effects, precisely the mean μ Ã (left) and the standard deviation σ (right), are reported for all the QoIs, that is, C t (top-left) and C v (top-right), ΔC t (bottom-left) and ΔC v (bottom-right).
presented in Figure 1, we notice higher linearity and monotonicity for the physical parameters rather than for the geometrical ones.
Using the deflated parameter space we perform a Sobol' analysis, fixing in this case the sample size equal to N ¼ 150 in order to obtain reliable results and imposing H in ¼ 0:45 6 and Pm 50 ¼ 1 mmHg 6 without affecting the output variance.The results of Sobol's analysis are shown in Figure 9.We show one plot for each of the QoIs, where we report the calculated main effect and total effect related to each input variable.The main trends already highlighted by the analysis of elementary effects are confirmed.The physical parameters mostly influence the mean-QoIs, whereas the geometrical ones affect the total variations.Moreover, for the input variables affecting the vascular geometry, we remark a significant difference between the total and the main effects, suggesting that the total influence of these variables is more effective when combined with the variation of other parameters.
Finally, the elementary effects and the Sobol's analysis are directly compared in Table 3. There, we compare the ranking of the input variables computed by the two methods and we also report the computational cost of the two approaches, measured by the number of input-output evaluations required to perform the analysis.We note that each evaluation involves the numerical solution of the models (ℱ þ ℋ þ O) described in (1) and (5).Table 3 illustrates one of the main methodological results of this work.It shows that, for the problem at hand, the elementary effects and the Sobol's analysis provide information about the sensitivity of the model that is qualitatively similar.However, as expected, the former method requires a computational cost that is about four times lower than the latter (3.9 reduction).More precisely, the elementary effect method uses r Ã M þ 1 ð Þ¼360 model evaluations, while Sobol's required N M þ 2 ð Þ¼1050 simulations.

| Sensitivity analysis varying physical and geometric parameters by groups
Finally, the Morris approach for groups of factors is employed to rank, as a whole, physical parameters and geometric ones.This new analysis relies on a number of simulations of the oxygen transfer model comparable to the previous case, that is r ¼ 120 trajectories associated to g ¼ 2 groups, hence entailing a total number of input realizations equal to r gþ 1 ð Þ¼360 simulations.From the results reported in Figure 10, we observe that physical inputs have a larger impact on all the considered QoIs with respect to the geometrical ones, although the influence of the latter is not negligible.The interpretation of these results hinges on two observations: on one hand the number of physical parameters is more than twice that of the geometrical ones, and on the other hand, at least one among V max and C v,in is a highly influential factor on the selected output of interest.Another interesting fact is that the geometric inputs are more relevant for the total variation of the oxygen concentration.

| Impact of the parameters on the tissue oxygenation
With the previous results we have established a ranking between the input parameters through a rigorous SA, detecting the most influential input factors with respect to specific QoIs of the oxygen transport model.These results, however, do not directly inform us about how variations of the parameters perturb the actual oxygen maps locally.We provide a visualization of the general interplay between significant parameters and oxygen maps in Figure 11, where we visualize the relation between a QoI, namely ΔC t , the most relevant inputs according to Figure 8 and the variations of the oxygen maps.Precisely, along the rows of Figure 11 we vary the vascular metrics, and T A B L E 3 Ranking of the input factors with respect to Morris and Sobol' indices with a comparison of their computational cost in terms of input/output model evaluations.

Elementary effects method
Sobol' indices along the columns we vary the oxygen consumption V max .Variations are designed to mimic the transition from baseline normal conditions of the VME towards a plausible tumor scenario, where the oxygen consumption is reduced and the vasculature becomes less dense and more irregular.In each cell we report the value of ΔC t and, at the bottom of the table, we calculate the ratio of the QoI of each cell with respect to the baseline value, denoted by A 1,1 ¼ A 0 1,1 ¼ 1:14510 À3 , referring to the table reported in Figure 11, where the first four boxes are named as matrix A and the last four correspond to the matrix A 0 .
According to these results, we point out some important facts that will be thoroughly discussed later on.Figure 11 shows that, while the influence of V max on ΔC t is mostly linear, the response of the same QoI to the variations of vascular metrics is significant but it does not depend linearly on the inputs.This effect is clearly highlighted by the QoI ratios reported at the bottom of the table and also confirmed by the plots of Figure 12.Moreover, it appears that the two vascular metrics influence the oxygen distribution very differently.S=V affects the oxygenation more uniformly, while has a more localized effect.Interestingly, we see that the cells with indices A 1,1 and A 2,2 of the first half of the table are characterized by a similar ΔC t (QoI), but the oxygen distribution is very different due to the combined variation of V max and % In conclusion, this result illustrates the complex relation between the vascular network topology and the oxygen maps.

| DISCUSSION
This work provides relevant results on two different levels, the methodology and the modeling of the VME.For the former, we propose reliable and efficient SA approaches for parametrized problems with microstructure.For the latter, we quantify the relevance of physical and geometric factors on the oxygenation of the VME.
For the methodological results, we show that, in our case, the computation of the sensitivity indices using Morris approach provides a reliable estimate of the relative ranking of all parameters.When compared to more reliable variance-based approaches, of which we consider here the well-known Sobol's indices method, the former does not underestimate the sensitivity of any input: in fact, according to Table 3, it only slightly overestimates the sensitivity of some F I G U R E 1 1 ΔC t estimates obtained varying each of the metrics and the physical parameter V max with respect to a fixed factor scale; the contours highlight the variation of the output with respect to the reference configuration provided imposing V max ¼ 1:8 Á 10 À4 mL O2 =mL B and considering a vascular network associated with % 75% and S=V ¼ 6600.The first four boxes of this table are named as matrix A and the last four are named as matrix A 0 .Correspondingly, at the bottom A ij =A 11 denotes the ratio between the values of the table, using the usual matrix notation.
of the input factors.Nonetheless, the Morris method provides a significant advantage in terms of the computational cost.We believe that this conclusion is general, at least in the framework of the multi-physics models of the VME, such as the one proposed here and denoted before by (ℱ þ ℋ þ O).On the basis of these considerations, we have shown that the Morris method can be used to perform sensitivity analysis of many QoIs related to the tissue microenvironment, for example in the context of radiotherapy, 6 chemotherapy, 56 immunotherapy, 57 or even hyperthermic treatment of tumors. 58,59rom the modeling standpoint, the previous SA 5 on the microvascular fluid balance combined with the results presented here, confirms the importance of interactions between parameters.This further motivates our efforts to move in the direction of multiphysics models of the microvascular environment.Simplified approaches mainly based on oxygen diffusion may not capture the importance of these interactions.This observation becomes even more relevant for those applications where the modeling complexity increases, such as the analysis of local drug delivery to the VME. 56rom the standpoint of applications, the main result of this study consists in showing that the mean value of the interstitial and vascular oxygen concentration is mostly affected by the physical parameters of the model.We also provide a clear ranking of the most influential parameters, reported in Table 3.This conclusion is also clearly confirmed by the application of the Morris method by groups, where we quantify the effect of two groups of inputs, the physical parameters and the geometric parameters.Figure 10 shows that the influence of the physical parameter cluster dominates in all the QoIs on that of the geometric parameters.This analysis confirms the trend already reported by the SA of individual inputs.For C t and C v the superiority of the physical parameters with respect to the geometric ones is evident.However, this result must be interpreted keeping in mind that the dimension of the two clusters that are compared is not well balanced, consisting of five items in the group of physical parameters and only two items in the group of geometric parameters.Figures 8 and 11 show that the influence of the physical parameters on the QoIs is mostly linear, while the relation between the geometric inputs and the oxygenation outputs is nonlinear, overall not easily predictable.Figure 11, in particular, shows that the parameters % and S=V strongly affect the spatial distribution of the oxygen delivered to the microvascular environment.
This study also suggests that the variations of hematocrit have generally a negligible influence on the oxygenation of the tissue, and this statement is in fact true for all the QoIs considered.This result can be explained by observing that the range of variation of hematocrit (H in , see Table 2) ensures in any case a good blood oxygenation.Abrupt variations of dissolved oxygen concentration are prevented and damped by the abundant oxyhemoglobin available.We also observe that the model is stationary, as a result, transient and local variations of hematocrit are not accounted for.
A possible limitation of these results consists in not including the vascular radius in the input parameters of the SA.We recall that the model used here accounts for variable vascular radii, according to the algorithm presented in Section 4.2.However, the variable radius distribution was the same for all the simulations.This study could be extended by introducing a suitable parametrization of the vascular radius distribution, by one or more parameters that could then become inputs of the sensitivity analysis.Nevertheless, we envision some possible drawbacks of this generalization.Primarily, adding the radius to the set of geometric parameters would make it significantly more difficult to sample independently the vascular density, the vascular voids and the radius as these three quantities are naturally interacting.As a result, there is the risk that adding the vascular radius and its distribution to this analysis would play a confounding factor.Moreover, this extension may also increase the computational cost of the whole workflow, requiring eventually the introduction of reduced order models to guarantee its feasibility.
The results of this SA have also an impact in parameter identification, by providing guidelines on what parameters can be reliably determined on the basis of some observations.In this context, we point out that the combination of C t and C v is related to the oxygen measurements coming from imaging data, primarily functional MRI.Then, our study suggests that imaging data and the corresponding average oxygen level at the voxel scale are primarily related to the physical state of the microenvironment, and only mildly affected by the microvascular architecture.

| CONCLUSIONS AND PERSPECTIVES
The main result of this work is the quantification of the relative role between geometric and physical factors affecting the VME.First, we have assessed independently all the input factors of a multi-physics model, ranking them with respect to different quantities of interest generally related to tissue oxygenation.This analysis showed that the interactions among factors are important and that the system response to geometrical inputs is typically nonlinear.Second, we have quantified the relative sensitivity of two groups of input parameters, the ones affecting the transport and diffusion phenomena at the basis of oxygen transfer, and the ones affecting the microvascular architecture.As groups, these two main categories of input factors play significant roles, although for oxygenation the physical parameters dominate over the geometrical ones.However, we have the feeling that the scenario may be substantially different for quantities of interest related to other phenomena such as drug release.
In conclusion, we observe that modeling the VME for cancer therapies is reaching new perspectives, see for example, 16,17 where the direct/forward modeling approach is part of an outer loop involving sensitivity analysis, uncertainty quantification, optimization, or inverse modeling.This study is our first step in this direction, using the multiphysics model involving microvascular flow (ℱ), hematocrit transport (ℋ) and oxygen transfer (O).In the near future, we foresee an expansion of it leveraging on the interaction of in-vivo data, in-vitro experiments, and advanced numerical techniques, such as reduced order modeling, efficient sampling methods and/or multi-fidelity approaches that complement the high performance methods applied here for the discretization of the governing equations at the basis of the multi-physics model of the VME.

F
I G U R E 2 (A) Distribution of positive and negative seeds in the 2D domain S; (B) distribution of positive seeds after removing sequentially the closest one to each negative seeds, with respect to the Euclidean distance; (C) constrained Voronoi-mesh built from the remaining positive seeds.that the boundary endpoints with x ¼ 0, y ¼ 0 or z ¼ 0 represent the set of the vascular inlets ∂Λ IN , where higher values of the vascular pressure are set, while the remaining boundary vertices are the vascular endpoints ∂Λ OUT , where lower values of the vascular pressure are prescribed.The inlet branches radii are set to the mean value, equal to 4 μm.

F I G U R E 3
Example of the progression from uniformly distributed vascular networks to more disordered vascular networks, with an increased maximum extravascular distance (d max ), from left to right.

F
I G U R E 6 Distribution of interstitial oxygen concentration C t x ð Þ (right) along a diagonal line ϕ on the tissue slice associated to z ¼ 0, varying the refinement of the 3D mesh.The plot on the right is determined with respect to the normalized distance from the origin of the line.

F I G U R E 8
Statistics of the elementary effects associated to the interstitial oxygen concentration C t (top-left) and ΔC t (top-right) and to the vascular oxygen concentration C v (bottom-left) and ΔC v (bottom-right).
Sobol' indices associated to the interstitial oxygen concentration C t (top-left) and ΔC t (top-right), and to the vascular oxygen concentration C v (bottom-left) and ΔC v (bottom-right).F I G U R E 1 0 Statistics of the elementary effects for the groups of physical and geometrical parameters, associated to the interstitial oxygen concentration C t and ΔC t , and to the vascular oxygen concentration C v and ΔC v .

F
I G U R E 1 2 (Top row) Boxplot of C t with respect to SEEDS À ð Þ SEEDS þ ð Þ (left)and the vascular density S=V (right).(Bottom row) Boxplot of ΔC t with respect to SEEDS À ð Þ SEEDS þ ð Þ (left) and the vascular density S=V (right).