Uncertainty propagation and sensitivity analysis: results from the Ocular Mathematical Virtual Simulator

. We propose an uncertainty propagation study and a sensitivity analysis with the Ocular Mathematical Virtual Simulator, a computational and mathematical model that predicts the hemodynamics and biomechanics within the human eye. In this contribution, we focus on the e ﬀ ect of intraocular pressure, retrolaminar tissue pressure and systemic blood pressure on the ocular posterior tissue vasculature. The combination of a physically-based model with experiments-based stochastic input allows us to gain a better understanding of the physiological system, accounting both for the driving mechanisms and the data variability.


Introduction
The interest in patient-specific mathematical models applied to biomedical problems has greatly increased in the last years.In particular, the need for a better understanding and knowledge of quantities in the medical context has raised tremendously the complexity of the mathematical models employed to describe such physical systems.However, a crucial aspect to guarantee that the model and its numerical solutions are meaningful from the biomedical viewpoint is how inherent uncertainties are incorporated, as recently discussed for instance in [34].In this direction, several works that studied the impact of uncertainties in the domain of cardiovascular disease modelling showed particular promise for elucidating the complex interplay between hemodynamics, biomechanics, and electrophysiology.Examples include arterial hemodynamics [6,8,24,2], cardiovascular simulations [48,35,32], electrophysiology [28], possibly coupled with electromechanical simulations [20] and/or hemodynamics [3].
To the best of our knowledge, the eye's mathematical and computational modelling is still at its early stages, as recently reviewed in [14].Biomechanical and fluid-dynamical aspects are of particular relevance for several clinical conditions [18], but numerous factors influence their complex coupling, and the underlying mechanisms are still elusive.Also, there is an intrinsic difficulty of isolating these factors in a clinical setting and measuring their contribution [55].The present work focuses on the interaction between the main ocular vessels' hemodynamics, intraocular pressure and the retrolaminar tissue pressure, which is directly related to the cerebrospinal fluid pressure.Among several interesting contributions in this area, we mention those closely related to our work.The first mathematical model that simultaneously accounts for blood flow in the central retinal vessels, blood flow in the retinal microvasculature, retinal blood flow autoregulation, biomechanical action of intraocular pressure on the retinal vasculature, and time-dependent arterial blood pressure was introduced in [16].A theoretical model to study the effects of intraocular pressure elevation on the central retinal artery hemodynamics was proposed in [15] and extended to account for the central retinal venous hemodynamics and the retinal microcirculation in [4].However, none of them explicitly accounted for uncertainties and variabilities in the model parameters.Only a few modelling works include a stochastic analysis framework and focused on the production and drainage of aqueous humour flow [53] and its coupling with ocular hemodynamics in a simplified manner [41].
With these premises, we present in this contribution an uncertainty quantification and a global sensitivity analysis for the main parameters involved in the mathematical and computational framework called the Ocular Mathematical Virtual Simulator (OMVS) that we have developed [43,42].The clinical relevance of results provided by the OMVS is described in [45] or, more extensively, in [42,Chapter 13], and it has been confirmed by an independent population-based study including nearly 10000 individuals [54].The reduced version of the OMVS model employed in the present study originates from [16] for the retinal circulation.It has been extended to include a reduced model for blood flow perfusion in the lamina cribrosa.Preliminary findings from a simplified uncertain quantification study were published as a peer-reviewed conference abstract in [22].As a significant step forward, we present hereafter a detailed global sensitivity analysis, using the Sobol' sensitivity indices.
The paper is organized as follows.The mathematical and computational model is described in Section 2.1, the uncertainty quantification and sensitivity analysis approach is presented in Section 2.2 and the input data in Section 2.3.The results of our study are depicted in Section 3 and discussed in Section 4. Finally, conclusions and future perspectives are outlined in Section 5.

Methodology
In the next two sections we describe the deterministic mathematical and computational foundations of our study, as well as the uncertainty quantification (UQ) and sensitivity analysis (SA) methods we incorporated in the OMVS framework to account for the stochastic features of the system.

Mathematical and computational model
The OMVS is a complex modelling framework that couples hemodynamics, biomechanics, and fluid dynamics in the eye to visualize and estimate in a non-invasive way ocular biofluids and tissues characteristics that are difficult or not accessible with standard investigation methods.The contributions developed within this framework can be subsequently utilized to isolate single risk factors and quantify their influence on the multi-factorial disease process.
To achieve this goal, the full OMVS is designed with a multiscale architecture, that aims at preserving the natural systemic features of blood circulation, while providing detailed views on sites of particular interest from the clinical viewpoint, such as the lamina cribrosa.The lamina cribrosa is a sponge tissue in the back of the eye that has a crucial role from the hemodynamical and neurological viewpoints.Specifically, this membrane is thought to help maintain the balance between the pressure inside the eye (intraocular pressure, hereafter denoted IOP) and behind the eye (the retrolaminar tissue pressure directly influenced by the intracranial pressure, denoted RLTp), which may influence the ocular blood flow.In addition, the lamina cribrosa acts as a scaffold for the retinal ganglion cell axons and the central retinal vessels and feeds RGC axons through its vascular network.The IOP is easily measurable with a Goldmann applanation tonometer.This instrument is based on the Imbert-Fick principle, which affirms that the pressure inside a dry thin-walled sphere corresponds to the force required to flatten the sphere surface divided by the flattening area.
In practice, we developed three model formulations of increasing complexity: we started from a 0-dimensional reduced-order description of the system, progressively adding the coupling with a porous media model for the lamina cribrosa and finally incorporating the effects of the deformation of the ocular tissues.More precisely, the OMVS combines (see Fig. 1): 1. System I (Fig. 1a): a circuit-based (0D) model for blood flow in the retinal vasculature, central retinal artery (CRA), and central retinal vein (CRV); 2. System II (Fig. 1b): a three-dimensional (3D) porous media model for the perfusion of the lamina cribrosa; 3. System III (Fig. 1c): a 3D isotropic elastic model for the biomechanics of the lamina cribrosa, retina, choroid, sclera, and cornea.
In the present contribution, UQ and SA analyses require intensive evaluations of the physical-based model.Therefore, we have employed a reduced version of the OMVS, accounting for the hemodynamical description provided by System I (Fig. 1a), coupled to a 0D adaptation of System II (Fig. 1b).For the proposed study, System III (Fig. 1c) has not been considered.In this manner, the model (i) provides a multiscale hemodynamics overview of the overall system, while maintaining a relatively accessible mathematical complexity and low computational costs; and (ii) combines information on ocular sites for which quantitative data are availablee.g.blood flow in the central retinal artery -and crucial ocular areas that are not accessible with clinical imagese.g.lamina cribrosa perfusion.
The 0D reduced version of the OMVS, see Fig. 2, exploits the electric analogy to fluid flow in complex vascular network [10].In this context, electric potentials correspond to fluid pressure, electric charges correspond to fluid volumes, and electric currents correspond to volumetric flow rates; the resistors and capacitors represent hydraulic resistance and wall compliance, respectively.Writing the constitutive equations characterizing the circuit elements and the Kirchhoff laws of currents and voltages leads to a system of ordinary differential equations whose solution provides the time-dependent profiles of pressures at the circuit nodes and flow rates through the circuit branches.We emphasize that IOP plays a crucial role in the description of the vein collapsibility, which is modelled in the OMVS by Starling resistors [51].Namely, when the external pressure is higher than the internal blood pressure, veins collapse, therefore dropping down the blood flow.
The network is constructed as an extension of a previous model for the retinal circulation, proposed and validated in [16] and [5].The vasculature is divided into six main compartments: central retinal artery (cra), arterioles (r,a), capillaries (r,c), venules (r,v), central retinal vein (crv), and the lamina cribrosa (lc).Each compartment includes resistances (R) and capacitances (C).The intraocular segments (R cra,3 , R cra,4 , R r,v1 , R r,v2 , R crv,1 , R crv,2 ) are exposed to the IOP and the retrobulbar segments are exposed to the RLTp (R cra,1 , R cra,2 , R crv,3 , R crv,4 ).The explicit conservation and constitutive laws, as well as parameters involved in the description of the retinal circuit, follow directly from work by [16] and [5].To the initial model, we have added a simplified description of the hemodynamics in the lamina cribrosa, involving two resistors lcRin and lcR, and one capacitor C 5 , with the following values: lcRin = 78181.9mmHg s cm −3 , lcR = 23988.25mmHg s cm −3 , and C 5 = 0.000000753 cm 3 mmHg −1 .Also, in the original circuit, the external pressure on the resistances R cra,3 and R crv,2 is the effective stress exerted by the lamina on these vessels, which has been computed via a simplified fluid-structure interaction model [15] describing the CRA/CRV interaction with the lamina cribrosa.In the current version of the  model, we do not account for this contribution in a similar manner but rather adopt a simplified approach, in which the external pressure corresponds to IOP.Further extensions could incorporate this dependence, but our choice was dictated by the possibility of computing stresses directly from the System III component of the OMVS.
The reduced model thus obtained can predict the hemodynamics within the lamina cribrosa, the retinal vasculature, and the central retinal vessels based on the key inputs described in Tab. 1 and displayed as coloured dots in Fig. 2. For the proposed study we fixed the value of the pressure at cavernous sinus (blue node), while we choose as input variables for UQ and SA the systolic and diastolic pressure at the ophthalmic artery (SP and DP, red node), the intraocular pressure (IOP, yellow node) and the retrolaminar tissue pressure (RLTp, green node).
Regarding the outputs, and in light of the clinical application in view, we will focus on the quantities of interest in the UQ and SA analysis listed in Tab. 2. Note that these quantities of interest will reflect the behaviour of the system at different time instants through the cardiac cycle.
The mathematical model previously described has been implemented in OpenModelica [11], an open-source Modelica-based modelling and simulation environment intended for industrial and academic studies of complex dynamic systems.Model results have been obtained using DASSL [26] with a tolerance of 10 −6 , a time step of 10 −3 and total simulation time of 8 s.DASSL (Differential/Algebraic System Solver) is an implicit, high order, multi-step solver with a step-size control based on backward differentiation formula (BDF).These features allow it to be stable and fit to be used for a wide range of models.Its first development can be in found in [26].
The system reaches a periodic state after the first cardiac cycle.However, we consider our output the last simulated cardiac cycle.Then, we retrieve the CRA, CRV, and the lamina cribrosa blood flow at three specific instant during the last cardiac cycle, namely the peak systolic time, the end of the systole and the end of the diastole (see Fig. 3).The choice of these three particular time instant in the cardiac cycle is driven by their interest from a clinical perspective.Moreover, in view of the model validation, the measurements of the blood flow is very often taken at peak systole and end diastole [19].

Uncertainty quantification and sensitivity analysis approach
Uncertainty propagation.The construction of a reliable model for ocular biofluid dynamics involves several steps with inherent uncertainties, among which: parameter inference from uncertain experimental data, model personalization to the same subject data at different time instants or to different individuals, etc.Therefore, a major challenge is to assess how these sources of uncertainty impact the clinically relevant outputs of the simulation and ultimately affect the confidence in the model predictions.
In the present contribution, we adopted the following approach: for the set of key inputs of the reduced version of the OMVS model, the uncertainty is represented by a probability density function (pdf), which quantifies the probability of a given parameter to reproduce a specific observation.We next develop a forward uncertainty quantification, also known as uncertainty propagation, to investigate how these input uncertainties are propagated to the outputs via the computational model.Thus, the combined information between previous modelling knowledge and assumptions on the prior pdfs allows us to obtain the posterior pdf for the quantities of interest selected from the clinical perspective.Remark 1.A major challenge that needs to be addressed is the choice and representation of the prior pdf in light of the quantitative and qualitative information available from generally noisy data, see for instance [40] and further discussion in Sec.2.3.
Sensitivity analysis.This part aims to determine synthetic measurements of which key inputs are the most influential on the quantities of interest selected among the outputs of the computational model without making assumption on the model and taking into account the continuous nature of the input parameters.To this end, we adopt the stochastic framework of global sensitivity analysis, which considers the input parameters {X j } j∈{1,...d} to be random independent variables with uncertainty modelled by a probability distribution, and employed to compute the random output Y.We have not considered the sensitivity analysis of the dynamic process Y(t).The methods described hereafter can be adapted to time t, and we could study the sensitivity indices with respect to time.However, it is more useful -easier to interpret -to perform sensitivity analysis at important characteristics of the response time series, in our case, peak systole, end systole, and end diastole, as discussed in Sec.2.1.
To quantify the influence of the variations of X j on the variations of Y, we compute the so-called Sobol' sensitivity indices originally proposed in the seminal paper [50], see also [27].More precisely, we define the first-order indices as where (i) var denotes the variance and E the expected value; (ii) var[Y] corresponds to the variability of Y with respect to the overall uncertainty including non-linear effects; (iii) var[E(Y|X j )], the variance of the conditional expectation E(Y|X j ), corresponds to the main or first order effect of X j ; it means that if Y is sensitive to X j , E(Y|X j ) is likely to vary a lot and hence var[E(Y|X j )] as well.
Another useful index is the total Sobol' index, defined as where X (− j) = (X 1 , . . ., X j−1 , X j+1 , . . ., X d ) and S − j denotes the sum of the indices where X j is not involved.Additionally, the effect due to specific interactions between the j th and the k th factors (k j) can be measured by second-order Sobol' indices, and so on for high-order interactions, see for more details [27].
Several approaches have been proposed to numerically compute these sensitivity indices, as reviewed for instance in [27].In the present work, we adopted the following two strategies: (i) a Monte Carlo-type approach and an estimator proposed in [46] on the basis of a combinatoric argument, and (ii) a Fourier amplitude sensitivity test (FAST) [47], which is a spectral method based upon the Fourier decomposition of the model response.The computational cost of the first approach for first order and total order Sobol' indices is of (d + 2)n model evaluations [46], where d is the input space dimension and n is the sample size; as for the FAST method the computation cost for first and total order indices is of d n model evaluations [47].
This variance-based approach implies intensive sampling, but it allows us to explore the input factors' full uncertainty ranges.We provide the two strategies as a way to ensure the reliability of our estimates.Indeed, given a sampling size, they may vary when repeating the estimations, and, in the Monte-Carlo approach, they can even be negative, although we used an implementation that mitigates this issue.
All the results on the UQ and SA analysis presented hereafter are carried out exploiting the Python statistical library Open-TURNS [1].

Input data
Both mathematical methods described before need as prior knowledge the statistical distribution of the input.We detail and critically discuss in the sequel the choices we propose, based on assumptions deduced from the experimental and clinical literature.Recall that the input of our model are systolic blood pressure (SP), diastolic blood pressure (DP), IOP and RLTp.
Blood pressures.For SP and DP data we refer to the paper of Sesso and co-authors [49] where these two quantities showed a normal distribution.In particular SP has a mean of 124.1 mmHg and a standard deviation of 11.1 mmHg (Fig. 4a), whereas DP has a mean of 77.5 mmHg and a standard deviation of 7.1 mmHg (Fig. 4b).Using SP and DP as input parameters raises an issue for the sensitivity analysis, notably the Sobol' index study, namely that the inputs have to be assumed independent, see Sec. 2.2.For IOP, RLTp, we can make this assumption reasonably, however, this is not valid for SP and DP [13].To overcome this problem we used Mean Arterial Pressure (MAP) as input variable in the Sobol study that has a normal distribution with mean 93 mmHg and standard deviation of 7.6 mmHg.In this case we also reconstructed the MAP distribution starting from the normal distribution of SP and DP and the following relationship: to check its normal distribution assumption (Fig. 4c).From Eq. (2.3) and from the correlation assumption between SP and DP made in [13], we can then reconstruct SP and DP starting from the MAP.Finally, the MAP is independent of the other two inputs (IOP, RLTp) and can be used as a sensitivity analysis parameter.
Intraocular pressure.There is a considerable discussion about normal and lognormal probability density functions in literature, and which of the two can better represent biological phenomena [25].
To explain our modelling choices we consider different IOP distributions using the data recovered from the clinical work of Suh and collaborators [52] (mean µ normal = 14.7 mmHg, standard deviation σ normal = 2.8 mmHg).Starting from these values, we have computed the mean and the variance both for Gaussian distribution and for a lognormal behaviour using the following formulas [21] to be consistent with the data: We performed a comparative analysis using three different IOP distributions based on the same clinical data, which have been defined above.In particular, we compare a normal, a truncated normal and a lognormal distribution.
Fig. 5 highlights that the normal probability density (blue) is going beyond some physiological constraints for healthy patients such as IOP > 5mmHg.The truncated normal probability density function (green) does not show this issue; however, this IOP distribution presents an abnormal cut in the left tail, making it ineffective for a sensitivity analysis study.For our simulations, with the data provided by Suh et al. [52], the lognormal distribution seems to be the more natural one.Our choice is dictated by the fact that we want to avoid miscalculation due to unphysiological input parameters that may lead to unrealistic discontinuities in the simulation results.This lognormal assumption has also been accepted in other works [56].For this reason we utilized the log normal distribution based on the population based study operated by Suh et al. [52] (Fig. 6).
Retrolaminar tissue pressure.For the RLTp we used a normal distribution (Fig. 7) with mean µ normal = 9.5 mmHg and standard deviation σ normal = 2.2 mmHg [29].In this case, no issues regarding the independence with other input or the sampling of unphysiological values occur.Remark 2. Our choices for the input probability density functions are dictated by clinical measurements reported in the literature, in particular [49] for the normal distribution of MAP, [52] for the lognormal distribution of IOP, and [29] for the normal distribution of RLTp.For the blood pressure and the retrolaminar tissue pressure, all the clinical literature refers to a normal probability density function of these inputs, whereas for what concerns the IOP, the prior knowledge on the uncertainty distribution is still a matter of debate.The rationale behind our specific choices is that we aim to avoid unphysiological input values that may lead to unrealistic discontinuities in the simulation results.It would also be possible in the future to estimate the input probability density functions in a more patient-specific manner, by using given repeated clinical measurements of patient-specific targets, such as systolic and diastolic blood pressure, intracranial pressure, intraocular pressure etc.

Results
In this section, we present two virtual studies using the OMVS described in Sec.2.1 and employing the two methods described in Sec.2.2.The input data assumptions have been discussed in Sec.2.3.In particular, in the first study, we will use uncertainty propagation, while in the second one, we will run a sensitivity analysis on the model.

Uncertainty propagation study
We completed three different sets of N = 10000 evaluations of the model using only the IOP as a stochastic input distribution.For the RTLp we fixed its value for the three sets at 9.5 mmHg (mean value provided by [29]).These three sets of evaluations differ from each other by the blood pressure value imposed; we selected three cases of clinical interest -in the same spirit as in [16]: 1. baseline subjects with a systolic/diastolic blood pressure of 120/80 mmHg; 2. low blood pressure subjects with SP = 100 mmHg and DP = 70 mmHg; 3. high blood pressure subjects with SP/DP = 140/90 mmHg.
Recall that the outputs on which we focus for the UQ study are listed in Tab. 2.
Numerical simulations.Fig. 8 shows the computed probability density functions for the three different locations we selected (CRA, CRV and LC), for three different time instants (peak systole, end systole and end diastole), and for the three populations of clinical interest (low, baseline and high blood pressure).We also report the simulated mean and standard deviation for each computed output in Tab. 3. Thus, we highlight that: • as expected, in each location and for all time instants, blood flow values decrease as we move from low to baseline and high blood pressure populations; • for the simulated CRA blood flow, the three populations (baseline, low and high) have distinct pdfs at all evaluated time instants (peak systolic, end systolic, end diastolic), and present a significant asymmetry; • for the simulated CRV blood flow, we have well established different pdfs at peak and end systole.We remark that within the veins we have a delayed peak of blood flows and CRV es values higher than CRV ps values (Tab.3), in good agreement with experimentally observed patterns in time velocity curves acquired with Doppler imaging.The CRV blood flow pdf at end diastole exhibits a peculiar shape: a peak of frequencies for high values and a plateau in the number of realizations for low values separated by an almost empty frequency area of relatively middle values.Even if the clinical interpretation of these distributions is difficult, it is interesting to observe that these results suggest a different repartition of frequencies between the three cases, in particular, for low blood pressures, the high values peak is narrower and the plateau is wider and with more realizations than the high blood pressure case.This fact is confirmed by the boxplot in Fig. 9a: the tail of CRV blood flow low values is within the first and third quantile range for the low blood pressure case, whereas this tail is composed just by outliers for baseline and high blood pressure populations; • for the LC blood flow, the OMVS suggests a similar analysis to CRV blood flow.The pdfs at peak, and end systole are distinct, whereas at end diastole it shows a peak of frequencies at low values and a more uniform distribution elsewhere.Also, Fig. 8i points out that for only for the low blood pressure virtual population a second peak and a considerable high number of frequencies can be identified at high values.Similarly to the CRV analysis, we propose the boxplot for the output LC ed in Fig. 9b: as predicted by the pdf.The tail of high values in the low blood pressure population is within the first and third quantile range.In contrast, the same tail consists of only outliers for baseline and high blood pressure populations.Finally, the computed LC blood flow variability is considerably lower than for the other two outputs (CRA and CRV) as reported in Tab. 3. [

Sobol' index study
In this study, we utilize as random input variables the IOP, the retrolaminar tissue pressure (RLTp), and the SP and DP to compute the Sobol' indices using the distribution introduced in Sec.2.3.We compute the Sobol' first and total indices using the Saltelli algorithm [46]        Figure 8: Output probability density functions predicted by the OMVS for three virtual populations: baseline subjects with normal SP and DP values (blue), low blood pressure subjects (orange), and high blood pressure subjects (green).and the FAST first and total indices [47].We performed 5 analysis with increasing N = [1000, 2000, 5000, 7500, 10000].The stopping criteria is based on the absolute iterative error: max where O = CRA ps,es,ed , CRV ps,es,ed , LC ps,es,ed , idx are the first or total index for two consecutive choices (n > n * ) in N. The final figures proposed have been obtained with n = 10000 where the error was less than 4 • 10 −2 for all input indices and using both algorithms to compute the results (see Tab. 4).
Numerical simulations.For what concerns the Sobol' indices, Figs.10a and 10b point out that the CRA blood flow at peak and end systole is highly dependent on the value of the MAP, whereas the influence of IOP and RTLp remains minimal.In Fig. 10c, the results provided by the OMVS suggest that the CRA blood flow at end diastole is highly influenced by the IOP, moderately by the MAP, and almost negligibly by the RLTp.Moreover, in this case, we notice that the total order and first order, especially for the MAP, are significantly different, meaning that high order interactions among parameters contribute to the variance of this output.For the CRV blood flow, Fig. 10d illustrates that there is a high dependency of CRV ps on MAP and only a moderate one from IOP. CRV es depends mainly on MAP (Fig. 10e), and CRV ed is highly dependent on the IOP and only mildly on the other two inputs (Fig. 10f).Also, for CRV blood flow at end diastole, we notice that high order interactions occur, especially for MAP and IOP. Figure 9: Boxplots of the end diastolic blood flow in the CRV and the lamina cribrosa predicted by the OMVS for three virtual populations: baseline subjects with normal SP and DP values, low blood pressure subjects, and high blood pressure subjects.

Input
Monte-Carlo approach [46] FAST method [ For the lamina cribrosa blood flow, the results provided by the OMVS sensitivity analysis suggest that the MAP is the dominant factor with a moderate influence of the IOP only at end diastole (Fig. 10g, 10h, and 10i).
The FAST indices (Figs.11) suggest that CRA ps (Fig. 11a), CRA es (Fig. 11b) and LC ed (Fig. 11i) depend mainly on MAP and mildly on IOP and RLTp.For what concerns the CRA ed (Fig. 11c) and CRV ed (Fig. 11f), they are highly influenced by the IOP, moderately by MAP and negligibly little by the RLTp.In this case the first and total order -especially for the IOP -are very different, implying that high order interactions are quite relevant to explain the variability of the output.Fig. 11d suggests that CRV ps has a high dependency on MAP and a moderate one on IOP.For CRV es (Fig. 11e), LC ps (Fig. 11g) and LC es (Fig. 11h), the sensitivity analysis results show that their variability is almost solely due to changes in MAP.  4 Discussion.
Hemodynamics in the ocular posterior tissue vasculature results from the combined effects of different factors.Specifically, ocular blood flow is driven by the difference between arterial and venous blood pressure, is impeded by IOP and RTLp (directly related to cerebrospinal fluid pressure) and is modulated by vascular regulation.Understanding these complex interrelated effects represents a major challenge when interpreting results from various clinical studies in ophthalmology.The current contribution introduces a networkbased model in the framework of the Ocular Mathematical Simulator, that couples retinal blood circulation with a simplified description of the LC hemodynamics.Further on, the model is employed to theoretically assess the relative contribution of IOP, RTLp and MAP stochastic variations on several clinically meaningful outputs characterizing hemodynamics in the CRA, CRV and LC, by means of uncertainty propagation methods and variance-based sensitivity indexes.
The results suggest that the hemodynamic response of the CRA, CRV and LC vasculature to variations in IOP and RTLp presents noticeable differences among individuals with different blood pressures, that strongly influences all the computed outputs.These model predictions are in good agreement with the experimental findings in [54] and several theoretical studies [16,15,4], designed to elucidate how and to which extent blood pressure is influencing the distribution of ocular hemodynamics.Note that in a clinical setting, there is an intrinsic difficulty of evaluating the contribution of each factor and the complex relationships among them.
The use of a lognormal distribution for the IOP translates into a non trivial interpretation of the output, especially at the level of the CRV and LC due to the non-linear character of the model and the complex interplay between factors (as pointed out in Sec.3.2).In this context, mathematical models are crucial to reproduce these mechanisms and help to unveil their interpretation.
We have compared the CRA blood flow simulated results to other data in the literature.Our results -Tab.3, in particular baseline  To compare our results with other clinical and mathematical studies that are more focused on the CRA blood velocities, we set the hypothesis of CRA diameter of about 160 µm [7,23].Using this assumption, our simulations provide similar values for CRA blood flows than the one measured by Harris et al. [19] (CRA ps = 120.6 µl/min, CRA ed = 30.1 µl/min) and the three virtual populations simulated by Guidoboni et al. [16] (baseline: CRA ps = 119.4µl/min, CRA ed = 33.7 µl/min; low: CRA ps = 95.3 µl/min, CRA ed = 28.9µl/min; high: CRA ps = 142.3 µl/min, CRA ed = 41.0 µl/min).This comparison shows the quality of the results, albeit the simplicity of the model we have employed for our study.From the SA (Sec.3.2) we evince that, as expected, the CRA blood flow depends mainly by the MAP, and only at end diastole -when the arterial pressure is at minimum -the IOP is affecting the CRA results.
For the CRV, we compared the output CRV es results (Tab.3) with total venous blood flow measurements reported in the literature.The simulated baseline mean value (= 52.1 µl/min) agrees with the experiments performed by Garcia et al. [12] (64.9 ± 12.8 µl/min) and Feke et al. [9] (80±12 µl/min for age 25-38 group, 73±13 µl/min for age 54-58 group).From a qualitative viewpoint, the analysis of the low values tail proposed in the numerical results paragraph may have an interesting physiological interpretation.As described in Section 2.1, IOP has a non-linear effect on the retinal vasculature, particularly on the venous part.Following this reason, the CRV high values peak may represent the natural state when the IOP is lower than the venous blood pressure, whereas the low values plateau denotes the collapse state.This statement is consistent with the previous analysis where we found a uniform distribution for low blood pressures and an important frequency in the peak for high blood pressures.The SA supports this prediction: the dependency to IOP is significant not only at end diastole, but also at peak systole.In contrast with the CRA, the CRV is a venous vessel, which blood pressure is lower therefore more easily influenced by external pressures (e.g.IOP).These indices show also significant differences between the first-order and the total-order index, which means that there are high-order interactions between these parameters.This fact is not surprising but it appears as a consequence of incorporating the Starling resistor effect, [51], which is a crucial requirement to retrieve clinical data, as reviewed in [14].
For the lamina cribrosa hemodynamics, we highlight the fact that all mathematical results are crucial in the investigation of disease because non-invasive measurements are not available nowadays for this tissue.Following similar consideration made for the CRV, we notice an interesting physiological interpretation following the uncertainty propagation study (Sec.3.1).For the low blood pressure population, which is notably the most at risk for ocular neuropathies, the results may suggest that an overperfusion of the lamina -with respect to the mean values of that population -occurs in more cases than for the other two population.The simplicity of the model, as confirmed by the low variability in the LC blood flows, does not allow us to clarify this fact.The SA (Sec.2.2) supports this analysis.Sobol and FAST indices show a high dependency only due to the MAP, which can be not so intuitive as for the CRA or the CRV, indeed we know that the pressure gradient across the LC (IOP − RLT p) may influence the hemodynamics [44].Further analysis by adding a three-dimensional hemodynamical and biomechanical description of the lamina would certainly help in this investigation, in particular to notice the impact of IOP and RLTp on LC blood flow.
Finally we discuss the estimates of the first and total indices using the Monte-Carlo and FAST approaches.For all considered quantities of interest, the two indices are suggesting similar outcomes, which allows us also to cross-validate the results of our analysis.The main differences concern CRA ed (Figs.10c and 11c) and CRV ed (Figs.10f and 11f), in particular on the high order interactions between IOP and MAP -more emphasized for the FAST estimates.We privilege the results provided by the Monte-Carlo estimates which are unbiased compared to the FAST method.Also, the latter has the initial advantage of being more computationally efficient, but at the cost of extra assumptions of smoothness in the model [57] that we are not currently able to verify for the OMVS.

Conclusion
Thanks to its special connection to the brain and its accessibility to measurements, the eye provides a unique window on the brain, thereby offering non-invasive access to a large set of potential biomarkers that might help in the early diagnosis and clinical care of Neuro-Degenerative Diseases [17].The OMVS has already shown great potentiality to reproduce the ocular biomechanics and the hemodynamics [42].
Pursuing this concept, in this contribution, we have proposed an uncertainty propagation study and a sensitivity analysis to evaluate the impact of uncertainties on this ophthalmological virtual laboratory.First, we have set up a framework to perform a forward UQ analysis, which allowed us to evaluate the effects of the propagation of uncertainty from input to output.Second, we completed a SA study based on Sobol indices to capture the interplay between the different model parameters and their relative importance.Finally, we have assessed qualitatively and quantitatively this computational framework in view of clinical applications.
In the context of computational models, a coupled UQ/SA analysis is crucial for the scientific research, especially in biology and medicine.The use of such mathematical framework may be employed in the process of product development or diseases understanding, decreasing the number of physical tests necessary and therefore the economic cost.The low fidelity model employed already provided useful information for analysis, however this study must be pursued with higher fidelity models such as System II and System III, see Figs. 1b and 1c respectively.The current methodology could thus be further improved, in particular by (i) using a multilevel multifidelity estimator, as for instance the one developed in Dakota toolkit, see [32,30]; (ii) devising a reduced order modelling approach for the 3D elastic and poroelastic models, following the reduced basis framework [31,33], developed in the open source software Feel++ [36,37].Finally, considering sensitivity indices over time, as in [3], as well as characteristics instants, provides a promising perspective of the present work.

Figure 2 :
Figure 2: Schematic of the 0D circuit-based model (System I of the OMVS) employed in this contribution for UQ and SA.The coloured nodes are the input of the model that can be inferred from clinical measurements.Key name Unit Brief description DP mmHg Diastolic blood Pressure SP mmHg Systolic blood Pressure IOP mmHg IntraOcular Pressure RLTp mmHg RetroLaminar Tissue pressure

Figure 3 :
Figure 3: Specific instant during the last cardiac cycle (peak systole, end systole and end diastole) used to compute the clinical outputs of interest.Image edited from https://en.wikipedia.org/wiki/Cardiac_cycle.

Figure 4 :
Figure 4: Normal probability density functions of the input blood pressure variables
Peak systolic CRA blood flow.
End diastolic CRA blood flow.
Peak systolic CRV blood flow.
End systolic CRV blood flow.
End diastolic CRV blood flow.
End systolic LC blood flow.
End diastolic LC blood flow.

Fig. 10
report the Sobol indices using the Saltelli algorithm, while Fig.11report the FAST indices.
End diastolic LC blood flow.
(a) Peak systolic CRA blood flow.(b) End systolic CRA blood flow.(c) End diastolic CRA blood flow.(d) Peak systolic CRV blood flow.(e) End systolic CRV blood flow.(f) End diastolic CRV blood flow.(g) Peak systolic LC blood flow.(h) End systolic LC blood flow.(i) End diastolic LC blood flow.

Figure 10 :
Figure 10: Estimation of Sobol' indices using the Monte-Carlo approach [46].Red squares represent first order indices, while blue squares represent total order indices.

( a )
Peak systolic CRA blood flow.(b) End systolic CRA blood flow.(c) End diastolic CRA blood flow.(d) Peak systolic CRV blood flow.(e) End systolic CRV blood flow.(f) End diastolic CRV blood flow.(g) Peak systolic LC blood flow.(h) End systolic LC blood flow.(i) End diastolic CRV blood flow.

Figure 11 :
Figure11: Estimation of first and total indices using the FAST method[47].Red squares represent first order indices, while blue squares represent total order indices.

Table 1 :
Input parameters for the reduced OMVS.
CRA ps µl/min peak systolic CRA blood flow CRA es µl/min end systolic CRA blood flow CRA ed µl/min end diastolic CRA blood flow CRV ps µl/min peak systolic CRV blood flow CRV es µl/min end systolic CRV blood flow CRV ed µl/min end diastolic CRV blood flow LC ps µl/min peak systolic lamina cribrosa blood flow LC es µl/min end systolic lamina cribrosa blood flow LC ed µl/min end diastolic lamina cribrosa blood flow Table 2: Quantities of interest in the UQ and SA analysis.

Table 3 :
Simulated mean and standard deviation for all the quantities of interest in the uncertainty propagation study.

Table 4 :
Error based on Eq. (3.1) computed between two consecutive choices of n ∈ N.