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

1.   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 ^[1]. 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 ^[2,3,4,5], cardiovascular simulations ^[6,7,8], electrophysiology ^[9], possibly coupled with electromechanical simulations ^[10] and/or hemodynamics ^[11].

To the best of our knowledge, the eye's mathematical and computational modelling is still at its early stages, as recently reviewed in ^[12]. Biomechanical and fluid-dynamical aspects are of particular relevance for several clinical conditions ^[13], 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 ^[14]. 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 ^[15]. A theoretical model to study the effects of intraocular pressure elevation on the central retinal artery hemodynamics was proposed in ^[16] and extended to account for the central retinal venous hemodynamics and the retinal microcirculation in ^[17]. 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 ^[18] and its coupling with ocular hemodynamics in a simplified manner ^[19].

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 ^[20,21]. The clinical relevance of results provided by the OMVS is described in ^[22] or, more extensively, in [20,Chapter 13], and it has been confirmed by an independent population-based study including nearly $10000$ individuals ^[23]. The reduced version of the OMVS model employed in the present study originates from ^[15] 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 ^[24]. 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.

2.   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.

2.1.   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 Figure 1):

1) System Ⅰ (Figure 1a): a circuit-based ($0$D) model for blood flow in the retinal vasculature, central retinal artery (CRA), and central retinal vein (CRV);

2) System Ⅱ (Figure 1b): a three-dimensional ($3$D) porous media model for the perfusion of the lamina cribrosa;

3) System Ⅲ (Figure 1c): a $3$D 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 Ⅰ (Figure 1a), coupled to a $0$D adaptation of System Ⅱ (Figure 1b). For the proposed study, System Ⅲ (Figure 1c) has not been considered. In this manner, the model (ⅰ) provides a multiscale hemodynamics overview of the overall system, while maintaining a relatively accessible mathematical complexity and low computational costs; and (ⅱ) combines information on ocular sites for which quantitative data are available, e.g., blood flow in the central retinal artery, and crucial ocular areas that are not accessible with clinical images, e.g., lamina cribrosa perfusion.

The $0$D reduced version of the OMVS, see Figure 2, exploits the electric analogy to fluid flow in complex vascular network ^[25]. 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 ^[26]. 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 ^[15] and ^[27]. 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 ^[15] and ^[27]. 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.9$ mmHg s cm$^{-3}$, $lcR = 23988.25$ mmHg 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 ^[16] 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 Ⅲ 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 Table 1 and displayed as coloured dots in Figure 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 Table 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 ^[28], 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 ^[29] 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 ^[29].

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 Figure 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 ^[30].

2.2.   Uncertainty quantification and sensitivity analysis approach

2.2.1.   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 ^[31] and further discussion in section 2.3.

2.2.2.   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\in \{1, \ldots 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 section 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 ^[32], see also ^[33]. More precisely, we define the first-order indices as

where (i) ${\mathrm{var}}$ denotes the variance and $\mathbb{E}$ the expected value; (ii) ${\mathrm{var}}[Y]$ corresponds to the variability of $Y$ with respect to the overall uncertainty including non-linear effects; (iii) ${\mathrm{var}}[\mathbb{E}(Y|X_j)]$, the variance of the conditional expectation $\mathbb{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$, $\mathbb{E}(Y|X_j)$ is likely to vary a lot and hence ${\mathrm{var}}[\mathbb{E}(Y|X_j)]$ as well.

Another useful index is the total Sobol' index, defined as

where $X_{(-j)} = (X_1, \ldots, X_{j-1}, X_{j+1}, \ldots, 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\neq j$) can be measured by second-order Sobol' indices, and so on for high-order interactions, see for more details ^[33].

Several approaches have been proposed to numerically compute these sensitivity indices, as reviewed for instance in ^[33]. In the present work, we adopted the following two strategies: (ⅰ) a Monte Carlo-type approach and an estimator proposed in ^[34] on the basis of a combinatoric argument, and (ⅱ) a Fourier amplitude sensitivity test (FAST) ^[35], 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 ^[34], where $d$ is the input space dimension and $n$ is the sample size; as for the FAST method the computational cost for first and total order indices is of $d \; n$ model evaluations ^[35].

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 OpenTURNS ^[36].

2.3.   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.

2.3.1.   Blood pressures

For SP and DP data we refer to the paper of Sesso and co-authors ^[37] where these two quantities showed a normal distribution. In particular SP has a mean of $124.1 \, {\rm{mmHg}}$ and a standard deviation of $11.1 \, {\rm{mmHg}}$ (Figure 4a), whereas DP has a mean of $77.5 \, {\rm{mmHg}}$ and a standard deviation of $7.1 \, {\rm{mmHg}}$ (Figure 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 section 2.2. For IOP, RLTp, we can make this assumption reasonably, however, this is not valid for SP and DP ^[38]. 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 \, {\rm{mmHg}}$ and standard deviation of $7.6 \, {\rm{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 (Figure 4c).

From Eq (2.3) and from the correlation assumption between SP and DP made in ^[38], 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.

2.3.2.   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 ^[39].

To explain our modelling choices we consider different IOP distributions using the data recovered from the clinical work of Suh and collaborators ^[40] (mean $\mu_{normal} = 14.7 \, {\rm{mmHg}}$, standard deviation $\sigma_{normal} = 2.8\, {\rm{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 ^[41] 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.

Figure 5 highlights that the normal probability density (blue) is going beyond some physiological constraints for healthy patients such as $IOP > 5$ mmHg. 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. ^[40], 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 ^[42]. For this reason we utilized the log normal distribution based on the population based study operated by Suh et al. ^[40] (Figure 6).

2.3.3.   Retrolaminar tissue pressure

For the RLTp we used a normal distribution (Figure 7) with mean $\mu_{normal} = 9.5 \, {\rm{mmHg}}$ and standard deviation $\sigma_{normal} = 2.2\, {\rm{mmHg}}$ ^[43]. 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 ^[37] for the normal distribution of MAP, ^[40] for the lognormal distribution of IOP, and ^[43] 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.

3.   Results

In this section, we present two virtual studies using the OMVS described in section 2.1 and employing the two methods described in section 2.2. The input data assumptions have been discussed in section 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.

3.1.   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 \, {\rm{mmHg}}$ (mean value provided by ^[43]). 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 ^[15]:

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 Table 2.

Numerical simulations Figure 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 Table 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 (Table 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 Figure 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, Figure 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 Figure 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 Table 3.

3.2.   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 section 2.3. We compute the Sobol' first and total indices using the Saltelli algorithm ^[34] and the FAST first and total indices ^[35]. We performed 5 analysis with increasing $\mathcal{N} = \left[ 1000, \: 2000, \: 5000, \: 7500, \:10000\right]$. The stopping criteria is based on the absolute iterative error:

where $\mathcal{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 $\mathcal{N}$. The final figures proposed have been obtained with $n = 10000$ where the error was less than $4 \cdot 10^{-2}$ for all input indices and using both algorithms to compute the results (see Table 4).

Numerical simulations Figure 10 report the Sobol indices using the Saltelli algorithm, while Figure 11 report the FAST indices.

For what concerns the Sobol' indices, Figure 10a, b 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 Figure 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, Figure 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 (Figure 10e), and CRV_ed is highly dependent on the IOP and only mildly on the other two inputs (Figure 10f). Also, for CRV blood flow at end diastole, we notice that high order interactions occur, especially for MAP and IOP. 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 (Figure 10g–i).

The FAST indices (Figure 11) suggest that CRA_ps (Figure 11a), CRA_es (Figure 11b) and LC_ed (Figure 11i) depend mainly on MAP and mildly on IOP and RLTp.

For what concerns the CRA_ed (Figure 11c) and CRV_ed (Figure 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.

Figure 11d suggests that CRV_ps has a high dependency on MAP and a moderate one on IOP.

For CRV_es (Figure 11e), LC_ps (Figure 11g) and LC_es (Figure 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 network-based 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 ^[23] and several theoretical studies ^[15,16,17], 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 section 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—Table 3, in particular baseline where CRA_ps $= 73.3 \pm 1.0\, \mu l/min$, CRA_es $= 32.8 \pm 0.6\, \mu l/min$, CRA_ed $= 21.7 \pm 1.7\, \mu l/min$—are in agreement with the experimental results of Dorner et al. ^[44] (CRA mean blood flow of $38.1 \pm 9.1 \, \mu l/min$) and Riva et al. ^[45,46] (CRA mean blood flow of $33 \pm 9.6 \, \mu l/min$).

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 \, \mu m$ ^[44,47]. Using this assumption, our simulations provide similar values for CRA blood flows than the one measured by Harris et al. ^[30] (CRA_ps $= 120.6\, \mu l/min$, CRA_ed $= 30.1\, \mu l/min$) and the three virtual populations simulated by Guidoboni et al. ^[15] (baseline: CRA_ps $= 119.4\, \mu l/min$, CRA_ed $= 33.7\, \mu l/min$; low: CRA_ps $= 95.3\, \mu l/min$, CRA_ed $= 28.9\, \mu l/min$; high: CRA_ps $= 142.3\, \mu l/min$, CRA_ed $= 41.0\, \mu 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 (section 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 (Table 3) with total venous blood flow measurements reported in the literature. The simulated baseline mean value ($= 52.1 \, \mu l/min$) agrees with the experiments performed by Garcia et al. ^[48] ($64.9 \pm 12.8 \, \mu l/min$) and Feke et al. ^[49] ($80 \pm 12 \, \mu l/min$ for age 25–38 group, $73 \pm 13 \, \mu 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 ^[26], which is a crucial requirement to retrieve clinical data, as reviewed in ^[12].

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 (section 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 (section 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-RLTp$) may influence the hemodynamics ^[50]. 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 (Figures 10c and 11c) and CRV_ed (Figures 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 ^[51] that we are not currently able to verify for the OMVS.

5.   Conclusions

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 ^[52]. The OMVS has already shown great potentiality to reproduce the ocular biomechanics and the hemodynamics ^[20].

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 Ⅱ and System Ⅲ, see Figure 1b, 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 ^[6,53]; (ii) devising a reduced order modelling approach for the 3D elastic and poroelastic models, following the reduced basis framework ^[54,55], developed in the open source software Feel++ ^[56,57]. Finally, considering sensitivity indices over time, as in ^[11], as well as characteristics instants, provides a promising perspective of the present work.

Conflict of interest

All authors declare no conflicts of interest in this paper.