<figure><img src="image_1.png" /></figure>

HAL
open science

Model order reduction and sensitivity analysis for
complex heat transfer simulations inside the human
eyeball

Thomas Saigre, Christophe Prud'Homme, Marcela Szopos

► To cite this version:

Thomas Saigre, Christophe Prud'Homme, Marcela Szopos. Model order reduction and sensitivity analysis for complex heat transfer simulations inside the human eyeball. International Journal for Numerical Methods in Biomedical Engineering, 2024, pp.e3864. 10.1002/cnm.3864. hal-04361954v2

HAL Id: hal-04361954

https://hal.science/hal-04361954v2

Submitted on 9 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Model order reduction and sensitivity analysis for complex heat transfer simulations inside the human eyeball

Thomas Saigre¹†, Christophe Prud'homme¹, Marcela Szopos²

¹Institut de Recherche Mathématique Avancée, UMR 7501 Université de Strasbourg et CNRS

²Université Paris Cité, CNRS, MAP5, F-75006 Paris, France

†Corresponding author contact: saigre@math.unistra.fr

v1.0.0 2024-10-09 (44002a0)

## Abstract

Heat transfer in the human eyeball, a complex organ, is significantly influenced by various pathophysiological and external parameters. Particularly, heat transfer critically affects fluid behavior within the eye and ocular drug delivery processes. Overcoming the challenges of experimental analysis, this study introduces a comprehensive three-dimensional mathematical and computational model to simulate the heat transfer in a realistic geometry. Our work includes an extensive sensitivity analysis to address uncertainties and delineate the impact of different variables on heat distribution in ocular tissues. To manage the model's complexity, we employed a very fast model reduction technique with certified sharp error bounds, ensuring computational efficiency without compromising accuracy. Our results demonstrate remarkable consistency with experimental observations and align closely with existing numerical findings in the literature. Crucially, our findings underscore the significant role of blood flow and environmental conditions, particularly in the eye's internal tissues. Clinically, this model offers a promising tool for examining the temperature-related effects of various therapeutic interventions on the eye. Such insights are invaluable for optimizing treatment strategies in ophthalmology.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Three-dimensional biophysical model</td><td>3</td></tr><tr><td>2.1</td><td>Geometry of the human eyeball</td><td>3</td></tr><tr><td>2.2</td><td>Biomechanical non-linear continuous model and its linearization</td><td>4</td></tr><tr><td>2.3</td><td>Model parameters</td><td>6</td></tr><tr><td>3</td><td>Mathematical and computational framework</td><td>6</td></tr><tr><td>3.1</td><td>Continuous and discrete model</td><td>6</td></tr><tr><td>3.2</td><td>Reduced order modeling with the reduced basis method</td><td>8</td></tr><tr><td>3.2.1</td><td>Error estimates</td><td>10</td></tr><tr><td>3.2.2</td><td>Generation of the reduced basis</td><td>10</td></tr><tr><td>3.3</td><td>Numerical results</td><td>11</td></tr><tr><td>3.3.1</td><td>Mesh convergence</td><td>11</td></tr><tr><td>3.3.2</td><td>Scalability</td><td>11</td></tr><tr><td>3.3.3</td><td>Linearized model</td><td>12</td></tr><tr><td>3.3.4</td><td>Verifications of the reduced basis model</td><td>12</td></tr><tr><td>3.4</td><td>Validation and comparison with previous studies</td><td>14</td></tr></table>
<table><tr><td>4</td><td>Uncertainty quantification</td><td>16</td></tr><tr><td>4.1</td><td>Deterministic sensitivity analysis</td><td>17</td></tr><tr><td>4.2</td><td>Stochastic sensitivity analysis</td><td>17</td></tr><tr><td>4.2.1</td><td>Choice of the distributions</td><td>19</td></tr><tr><td>4.2.2</td><td>Uncertainty propagation</td><td>20</td></tr><tr><td>4.2.3</td><td>Results of the SA</td><td>21</td></tr><tr><td>5</td><td>Conclusion</td><td>22</td></tr><tr><td>A</td><td>Reproductibility of results</td><td>22</td></tr><tr><td></td><td>Acknowledgements</td><td>24</td></tr><tr><td></td><td>References</td><td>26</td></tr></table>

**Keywords**: mathematical and computational ophthalmology, heat transfer, validation, finite element method, real-time model order reduction, uncertainty quantification, sensitivity analysis, Sobol index analysis.

# 1 Introduction

The development of new technologies allows us to simulate more and more complex models in order to apprehend the world we live in. In this study, we will focus on a specific model: heat transfer inside the human eyeball. The temperature of the eyeball may have an impact on the distribution of drugs in the eye, partly due to the aging of the tissues [BBS20]. Hyperthermia is one of the most common treatments for eye tumors [Li+10], and understanding the mechanism of heat transfer could enhance the efficacy of ophthalmic treatments, such as laser therapy of the retina [Mas04].

Heat transfer is also a key factor in the study of the effects of electromagnetic radiation on the eye, as pointed out in [Hir+07; ON09]. The model, originally introduced in [Sco88] to examine temperature rises induced by exposure to infrared radiation, has been expanded upon in subsequent studies [Ng 06; NO07; ON08; Li+10] using diverse methods for computing heat transfer.

While invasive studies on animals have been conducted [PW05], non-invasive measurements on human subjects are scarce, complex to perform, and may yield inaccurate results [RF77]. Most studies focus on temperature measurements at the eye's surface [Map68; EYB89] but report significant differences and identify several sources of uncertainty. Alternatively, numerical simulations can provide complementary information. However, in order to guarantee the reliability of such results, a rigorous validation step is required.

The present contribution aims to contribute to these developments, by means of a mathematical and computational modeling approach, combined with a sensitivity analysis study performed thanks to a model reduction technique. The comparison with data available in the literature, obtained either by measurement on healthy subjects [EYB89] or by other simulations [Ng 06; NO07; Li+10] will ensure the validity of the approach.

In this model, numerous parameters, both biomechanical and geometrical, are involved. The present study concentrates on biomechanical parameters, in a large range that include potential extremal conditions. The variation of these parameters can have a significant impact on the results. To quantify their impact, we set up a framework to perform a forward uncertainty quantification study, complemented by a sensitivity analysis. Deterministic sensitivity analysis has already been performed in [Sco88; Ng 06; NO07; Li+10], using various numerical methods. In this work, we reproduce and extend these results, to incorporate the effect of blood flow, as suggested for instance in [Sco88]. We also run a global sensitivity analysis, that accounts for stochastic effects, and discriminate among different factors by means of Sobol's indices [Sob93]. The combination of deterministic and stochastic methods is an effective practice in the field of uncertainty quantification, and has been successfully used in recent various applications, such as [Dod+14] where the impact of uncertainties on the distribution of the electromagnetic field in the ocular tissue is studied; or more generally in the human head [Suš+22; Šuš+21]. To the best of our knowledge, this is the first time that such a study is performed in the context of bioheat transport in the tissues of the human eyeball.
<figure><img src="image_3.png" /><figcaption>Figure 1: Geometrical model of the human eye.</figcaption></figure>

While Sobol's indices are effective in measuring parameter impact and interactions, the complexity and the significant computational time of our model are very challenging. To overcome this, we adopt the certified reduced basis method [Pru+01; QMN16] to obtain a reduced model, maintaining its 3D nature while significantly reducing computational demands. This method aligns with the paradigm observed in patient-specific mathematical models applied to biomedical problems, ensuring a comprehensive approach involving data integration, model derivation, numerical solving, validation, and uncertainty quantification, as seen in mature research fields like cardiovascular simulations or cerebral hemodynamics. In ophthalmology, a similar paradigm is imperative due to the richness and heterogeneity of available data, requiring innovative approaches for diagnosis and monitoring.

More generally, the present work aims to contribute to the project Eye2Brain [Sal16], that has the ambitious objective to connect the cerebral and ocular environments and contribute in the long term to a better understanding of neurodegenerative diseases [Gui+20]. In this context, a model accounting for the combined effects of ocular blood flow and different ocular tissues was proposed in [Sal+]0 To incorporate inherent uncertainties and variability, an uncertainty propagation and sensitivity analysis on the component simulating the fluid flows in the eye was developed in [PSS21]. We here focus on the heat propagation phenomena, with the perspective of coupling the fluid and thermal contributions in future work.

The structure of the paper is the following. After the introduction, we describe in Section 2 the geometrical model describing the human eyeball, the biophysical model governing the heat transfer, as well as the parameters involved in the equations. Next, we present in Section 3 the methods developed to simulate the full and reduced models, including a step of verification and validation, to ensure that the mathematical and computational framework is correct. We report in Section 4 our results of the sensitivity analysis, using two methods: a deterministic one and a stochastic approach. All the methods are implemented in the open-source software Feel++ [Pru+23] and can be reproduced following guidelines described in Appendix A. Finally, conclusions and perspectives are outlined in Section 5.

# 2 Three-dimensional biophysical model

## 2.1 Geometry of the human eyeball

In this section, we describe the realistic three-dimensional geometry that will be used in the sequel. The model we employ in the present work stems from [Sal+], and was constructed using a CAD (Computer Aided Design) module from SALOME [CAS22]. Figure 1 shows a cut-away view along a vertical plane of the reconstructed eye anatomy.

The eye is composed of several regions, which have different physical properties. The original geometry contained five subdomains: the sclera, the choroid, the retina, the cornea and the lamina cribrosa. To have a better assessment of the thermal properties of each part, we further decompose the geometry as follows: (i) the cornea which allows heat transfer between the eye and the ambient air, (ii) the envelope of the eye composed of the sclera, the optic nerve, and the lamina cribrosa, (iii) the vascular beds namely the choroid and the retina, mostly composed of blood vessels, (iv) the anterior and posterior chambers, filled with aqueous humor, (v) the lens and (vi) the vitreous body filled with the vitreous humor, a transparent
<figure><img src="image_4.png" /><figcaption>Figure 2: Featured geometrical locations for the output of interest (pointwise temperature).</figcaption></figure>

liquid allowing the light to reach the retina. In the present model, the optic nerve domain is assumed to be homogeneous, the contribution of the inner vessels is not directly taken into account in heat transfer.

Several more simplified geometrical descriptions were already utilized in the literature to study heat transport in the eye; mostly in 2D [Sco88; Ng 06] or in 3D [NO07; Li+10]. In particular, the 3D model developed in [NO07] did not incorporate a detailed description of the vascular beds, although previous studies [Sco88] and our further sensitivity analysis pointed out the importance of the influence of the blood temperature on the heat distribution.

In order to compare in a first stage our results with previously reported findings [EYB89], we define on the front part of the cornea the geometrical center of the cornea (GCC), see Figure 1, which is an imaginary line “cutting” the cornea horizontally. This region is interesting because this part of the eye is easily accessible and the temperature can be measured non-invasively.

We focus on outputs of interest that are studied in the literature [Sco88; Ng 06]. These outputs are the temperature values at given locations or the mean temperature on a given domain. Precisely, we select on points present at the interface of two regions of the eye, as well as the mean temperature over the cornea. For a precise description of these locations, see Figure 2.

## 2.2 Biomechanical non-linear continuous model and its linearization

Based on Section 2.1, the geometry of the eye can be written as a disjoint union of different regions: $\Omega = \cup_{i=1}^{10} \Omega_i$, where $i$ is the index of the subdomain and $\Omega_i$ corresponds to the following regions: cornea, vitreous humor, aqueous humor, retina, iris, choroid, lens, sclera, lamina cribrosa, and optic nerve.

We focus on stationary heat transfer in this domain. Following [Sco88; Ng 06] the steady-state condition of the heat transfer in the human eye can be described by the following system

$$ \nabla \cdot (k_i \nabla T) = 0 \quad \text{in } \Omega = \cup_{i=1}^{10} \Omega_i \qquad (1a) $$

where:

* $i$ is the volume index (cornea, vitreous humor...),

* $T_i$ [K] is the temperature in the domain $\Omega_i$,

* $k_i$ [W m$^{-1}$ K$^{-1}$] is the thermal conductivity of $\Omega_i$.

We set the global thermal conductivity $k$ [W m$^{-1}$ K$^{-1}$] as a discontinuous piece-wise constant function: $k = k_i$ on $\Omega_i$. The boundary $\partial\Omega$ is decomposed as: $\partial\Omega = \Gamma_{amb} \cup \Gamma_{body}$ (see Figure 3), where $\Gamma_{amb}$ corresponds to the boundary region exposed to the ambient environment and $\Gamma_{body}$ the boundary of the internal domain. Denote by $n$ the outward normal vector to the domain $\Omega$. The following boundary conditions are adopted:

* To model the exchange between the eye and the ambient air, and incorporate radiative heat transfer we impose the following non-linear Neumann condition:

$$ -k \frac{\partial T}{\partial n} = h_{amb}(T - T_{amb}) + \sigma \varepsilon (T^4 - T_{amb}^4) + E \quad \text{on } \Gamma_{amb} \qquad (1b) $$
<figure><img src="image_5.png" /><figcaption>Figure 3: Description of the physical boundaries and interfaces of the domain $\Omega$.</figcaption></figure>

Three terms are present in this condition to describe different heat loss mechanisms occurring on the cornea: (i) The first term in the equation represents the convective heat transfer between the surface of the eye and the surrounding air. The parameter $h_{amb}$ [W m$^{-2}$ K$^{-1}$] is the air convective coefficient, and $T_{amb}$ [K] is the ambient temperature; (ii) the second term represents the radiative heat transfer between the surface of the eye and the surrounding environment, where the parameter $\sigma$ is the Stefan-Boltzmann constant ($\sigma = 5.67 \times 10^{-8}$ W m$^{-2}$ K$^{-1}$), and $\varepsilon$ [–] is the emissivity of the surface; (iii) the third term represents the heat loss due to tear evaporation. The parameter $E$ [W m$^{-2}$] represents the heat transfer rate due to evaporation, which depends on the environmental conditions and the tear film characteristics. This process causes a cooling effect on the surface of the eye, which can be significant in dry environments or cases of reduced tear production.

* To model the thermal exchanges between the eye and the body, we impose:

$$ -k \frac{\partial T}{\partial \boldsymbol{n}} = h_{\text{bl}}(T - T_{\text{bl}}) \quad \text{on } \Gamma_{\text{body}} \qquad (1c) $$

where the parameter $h_{bl}$ [W m$^{-2}$ K$^{-1}$] is the blood convection coefficient and $T_{bl}$ [K] is the blood temperature.

Finally, to ensure a continuous flow of heat flux and no temperature jump, we impose at the interface between two adjacent regions $\Omega_i$ and $\Omega_j$ the following condition:

$$ \left\{ \begin{array}{l} T_i = T_j \\ k_i(\nabla T_i \cdot \boldsymbol{n}_i) = -k_j(\nabla T_j \cdot \boldsymbol{n}_j) \quad \text{on } \partial\Omega_i \cap \partial\Omega_j \end{array} \right. \qquad (1d) $$

where $\boldsymbol{n}_i$ (resp. $\boldsymbol{n}_j$) denotes the outward normal vector to the domain $\Omega_i$ (resp. $\Omega_j$).

System (1a) - (1d) defines a non-linear problem, denoted $\mathcal{E}_{NL}$ in the sequel.

**Remark 2.1.** Note that the condition (1b) modeling radiative transfer is non-linear, because of the term in $T^4$, which requires a more complex treatment, both from the mathematical standpoint, for the reduced basis method; and from the numerical standpoint, due to extra computational cost. As an alternative, a linearization of the condition (1b) was proposed in [Sco88]:

$$ \sigma\varepsilon(T^4 - T_{\text{amb}}^4) = (T - T_{\text{amb}}) \underbrace{\sigma\varepsilon(T^2 + T_{\text{amb}}^2)(T + T_{\text{amb}})}_{=:h_r}, $$

which leads to a linear Robin condition. The value $h_r$ stands for the radiation heat transfer coefficient and is approximately equal to $6 \text{ W m}^{-2} \text{ K}^{-1}$ [Sco88].

Condition (1b) can hence be rewritten as:

$$ -k \frac{\partial T_i}{\partial \boldsymbol{n}} = h_{\text{amb}}(T - T_{\text{amb}}) + h_r(T - T_{\text{amb}}) + E \quad \text{on } \Gamma_{\text{amb}} \qquad (2) $$

The model described by Equations (1a)-(2)-(1c)-(1d) is further denoted $\mathcal{E}_L$.
<table><caption>Table 1: Parameters involved in the model, baseline values and ranges used in the sensitivity analysis.</caption><thead><tr><th>Symbol</th><th>Name</th><th>Dimension</th><th>Baseline value</th><th>Range</th></tr></thead><tbody><tr><td>T<sub>amb</sub></td><td>Ambient temperature</td><td>[K]</td><td>298</td><td>[283.15, 303.15]</td></tr><tr><td>T<sub>bl</sub></td><td>Blood temperature</td><td>[K]</td><td>310</td><td>[308.3, 312]</td></tr><tr><td>h<sub>amb</sub></td><td>Ambient air convection coefficient</td><td>[W m<sup>-2</sup> K<sup>-1</sup>]</td><td>10</td><td>[8, 100]</td></tr><tr><td>h<sub>bl</sub></td><td>Blood convection coefficient</td><td>[W m<sup>-2</sup> K<sup>-1</sup>]</td><td>65</td><td>[50, 110]</td></tr><tr><td>E</td><td>Evaporation rate</td><td>[W m<sup>-2</sup>]</td><td>40</td><td>[20, 320]</td></tr><tr><td>k<sub>lens</sub></td><td>Lens conductivity</td><td>[W m<sup>-1</sup> K<sup>-1</sup>]</td><td>0.4</td><td>[0.21, 0.544]</td></tr><tr><td>k<sub>cornea</sub></td><td>Cornea conductivity</td><td>[W m<sup>-1</sup> K<sup>-1</sup>]</td><td>0.58</td><td>-</td></tr><tr><td>k<sub>sclera</sub> = k<sub>iris</sub> =<br>k<sub>lamina</sub> = k<sub>opticNerve</sub></td><td>Eye envelope components conductivity</td><td>[W m<sup>-1</sup> K<sup>-1</sup>]</td><td>1.0042</td><td>-</td></tr><tr><td>k<sub>aqueousHumor</sub></td><td>Aqueous humor conductivity</td><td>[W m<sup>-1</sup> K<sup>-1</sup>]</td><td>0.28</td><td>-</td></tr><tr><td>k<sub>vitreousHumor</sub></td><td>Vitreous humor conductivity</td><td>[W m<sup>-1</sup> K<sup>-1</sup>]</td><td>0.603</td><td>-</td></tr><tr><td>k<sub>choroid</sub> = k<sub>retina</sub></td><td>Vascular beds conductivity</td><td>[W m<sup>-1</sup> K<sup>-1</sup>]</td><td>0.52</td><td>-</td></tr><tr><td>ε</td><td>Emissivity of the cornea</td><td>[-]</td><td>0.975</td><td>-</td></tr></tbody></table>

## 2.3 Model parameters

In the model presented in the previous section, many parameters are involved, but not all of them are directly measurable. Moreover, inherent uncertainties due to noise and individual variability must be taken into account in the modeling process. We therefore fixed in a first stage a set of baseline values, corresponding to the nominal values for the human body, according to the literature [Sco88; Ng 06] (see Table 1). In a second step, we split the total set of parameters into two subsets: a first part kept fixed to baseline values, and a second part that varies in a certain range (see Table 1). The aim is to perform a refined sensitivity analysis, that encompasses previously published studies [Sco88; Ng 06; NO07], and extends the analysis to a larger parameter space.

Specifically, we set the varying parameter space $D^\mu \subset \mathbb{R}^6$ as the Cartesian product of the intervals defined in the last column of Table 1. For the purpose of the sensitivity analysis, an element $\mu = \{T_{amb}, T_{bl}, h_{amb}, h_{bl}, E, k_{lens}\} \in D^\mu$ is called a parameter, and we denote $\bar{\mu}$ the baseline parameter, extracted from the corresponding column in Table 1. The dependence of the model concerning the parameter $\mu$ is emphasized by the notation $\mathcal{E}_L(\mu)$ and $\mathcal{E}_{NL}(\mu)$.

# 3 Mathematical and computational framework

This section outlines the mathematical and computational framework, including the variational formulation derivation, the high fidelity finite element method (FEM) resolution technique, and the construction of reduced basis metamodel. It is followed by the presentation of numerical results, verification and validation steps.

## 3.1 Continuous and discrete model

We first compute the variational formulation of the linearized model $\mathcal{E}_L(\mu)$ described in Section 2. Let $v \in H^1(\Omega)$ be a test function. As the union $\Omega = \bigsqcup_i \Omega_i$ is disjoint, we have:

$$ \int_{\Omega} -\nabla \cdot (k \nabla T) v \, dx = \sum_{i} \int_{\Omega_i} -\nabla \cdot (k_i \nabla T) v \, dx \quad (3) $$

Hence, using Green's theorem:

$$ \sum_i \int_{\Omega_i} -\nabla \cdot (k_i \nabla T) v \, dx = 0 \Leftrightarrow \sum_i \int_{\Omega_i} k_i \nabla T \cdot \nabla v \, dx - \int_{\partial \Omega_i} k_i \frac{\partial T}{\partial n_i} v \, d\sigma = 0 \quad (4a) $$

with boundary and interface conditions Equations (1c), (1d) and (2), we obtain
$$ \sum_i k_i \int_{\Omega_i} \nabla T \cdot \nabla v \, dx + \int_{\Gamma_{amb}} [h_{amb}T + h_rT] v \, d\sigma + \int_{\Gamma_{body}} h_{bl}Tv \, d\sigma = \\ \int_{\Gamma_{amb}} [h_{amb}T_{amb} + h_rT_{amb} - E] v \, d\sigma + \int_{\Gamma_{body}} h_{bl}T_{bl}v \, d\sigma \quad (4b) $$

The previous equation is equivalent to:

$$ a_L(T, v; \mu) = f_L(v; \mu) \quad (5a) $$

with:

$$ a_L(T, v; \mu) := k_{\text{lens}} \int_{\Omega_{\text{lens}}} \nabla T \cdot \nabla v \, dx + \sum_{i \neq \text{lens}} k_i \int_{\Omega_i} \nabla T \cdot \nabla v \, dx + \int_{\Gamma_{\text{amb}}} [h_{\text{amb}}T + h_rT] v \, d\sigma + \int_{\Gamma_{\text{body}}} h_{\text{bl}}Tv \, d\sigma \quad (5b) $$

$$ f_L(v; \mu) := \int_{\Gamma_{amb}} [h_{amb}T_{amb} + h_rT_{amb} - E] v \, d\sigma + \int_{\Gamma_{body}} h_{bl}T_{bl}v \, d\sigma \quad (5c) $$

The problem statement is therefore: for $\mu \in D^\mu$ given, find the output of interest $s(\mu) \in \mathbb{R}$ given by

$$ s(\mu) = \ell(T(\mu)), \quad (6) $$

where $T(\mu) \in H^1(\Omega)$ is solution to

$$ a_L(T(\mu), v; \mu) = f_L(v; \mu) \quad \forall v \in H^1(\Omega). \quad (7) $$

The functional $\ell$ returns the desired output of interest, which can be the mean temperature in a selected region e.g. $\ell(T(\mu)) = \frac{1}{|\Omega_{\text{cornea}}|} \int_{\Omega_{\text{cornea}}} T(\mu) \, dx$, or the temperature at a fixed point e.g. $\ell(T(\mu)) = \langle \delta_O, T(\mu) \rangle$.

**Theorem 3.1.** Let $\mu \in D^\mu$ fixed. The problem (5) is well-posed for $v \in H^1(\Omega)$: there exists a unique $T(\mu) \in H^1(\Omega)$ such that $a_L(T(\mu), v; \mu) = f_L(v; \mu)$ for all $v \in H^1(\Omega)$. If $T(\mu) \in C^1(\bar{\Omega}) \cap C^2(\Omega)$, then $T(\mu)$ is solution to problem $\mathcal{E}_L(\mu)$.

*Proof.* The result is a straightforward application of the Lax-Milgram theorem [EG21], and of the regularity of $T$. $\square$

**Remark 3.2.** The well-posedness of the fully non-linear problem $\mathcal{E}_{\text{NL}}(\mu)$ can also be obtained by the mean of a variational approach, in the spirit of [Mil93].

**High fidelity FEM resolution.** We present here the discretization approach and briefly describe the in-house computational framework we developed.

In the sequel, we set $V := H^1(\Omega)$ and focus on outputs of interest, $s_k(\mu)$, for $k \in [1, n_{\text{output}}]$ given by the formula $s_k(\mu) = \ell_k(u(\mu); \mu)$, where $\ell$ is a bounded linear form and $u(\mu)$ is the solution of Equation (7).

Denote by $V_h \subset V$ the approximate functional space of dimension $\mathcal{N}$, $h$ standing for the discretization of the space, for a finite element approach. The previous problem is equivalent to:

$$ \mathbf{A}_L(\mu)\mathbf{T}^{\text{fem}}(\mu) = \mathbf{f}_L(\mu) \quad (8a) $$

$$ s_k(\mu) = \mathbf{L}_k(\mu)^T \mathbf{T}^{\text{fem}}(\mu) \quad (8b) $$

with $\mathbf{A}(\mu) \in \mathbb{R}^{\mathcal{N} \times \mathcal{N}}$, $\mathbf{f}(\mu) \in \mathbb{R}^{\mathcal{N}}$, $\mathbf{L}_k(\mu) \in \mathbb{R}^{\mathcal{N}}$, and $k$ is the index of the output. The vector $\mathbf{T}^{\text{fem}}(\mu) \in V_h \simeq \mathbb{R}^{\mathcal{N}}$ is the solution, and $s(\mu) \in \mathbb{R}$ is the computed output.

More precisely, the steps run during resolution are given in Algorithm 1.
**Algorithm 1:** High fidelity resolution.

Input: $\mu \in D^\mu$

Construct $\mathbf{A}(\mu)$, $\mathbf{f}(\mu)$, $\mathbf{L}_k(\mu)$;

Solve $\mathbf{A}(\mu)\mathbf{T}^{\text{fem}}(\mu) = \mathbf{f}(\mu)$;

Compute outputs $s_k(\mu) = \mathbf{L}_k(\mu)^T \mathbf{T}^{\text{fem}}(\mu)$;

Output: Numerical solution $\mathbf{T}^{\text{fem}}(\mu)$ and outputs $s_k(\mu)$

To establish numerical results, we implement Algorithm 1 in the framework of the open-source library Feel++ [Pru+23]<fnref n="1" />, and specifically the heat toolbox<fnref n="2" /> where both models $\mathcal{E}_{\text{NL}}(\mu)$ and $\mathcal{E}_L(\mu)$ can be simulated, with both $\mathbb{P}_1$ and $\mathbb{P}_2$ piece-wise polynomials [EG21]. The solution strategy uses conjugate gradient method solver preconditioned by an algebraic multigrid method; in the context of $\mathcal{E}_{\text{NL}}$, non-linear iterations are required.

All the results presented in this document are available and can be reproduced, refer to Appendix A for more details. All subsequent computational simulations are performed on the same machine equipped with the following hardware: AMD EPYC 7552 48-Core Processor.

## 3.2 Reduced order modeling with the reduced basis method

We now introduce the reduced basis metamodel [Pru+01; RHP08; QMN16]. The goal of the reduced basis method (RBM) is to approximate the solution of the parametrized-PDE described by Equations (1a)-(2)-(1c)-(1d). For complex geometries and biomechanical problems, such as the one described in Section 2.1, numerical solving has a prohibitive cost, especially for studies of uncertainty quantification, requiring the resolution of the system for many parameters. We briefly present the implemented strategy, following [Pru+01].

Recall that $\mathbf{T}^{\text{fem}}(\mu) \in V_h$ can be written as $\mathbf{T}^{\text{fem}}(\mu) = \sum_{n=1}^{N} T_n^{\text{fem}}(\mu)\phi_{h,n}$, where $(\phi_{h,n})_{n \in [1,N]}$ is a basis of $V_h$.

The main idea of RBM is to construct a low-dimension subspace $V_N \subset V_h$, of dimension $N$ with $N \ll N$, such that the approximation error is small: $\|\mathbf{T}^{\text{fem}}(\mu) - \mathbf{T}^{\text{rbm}}(\mu)\| \le \varepsilon_{\text{tol}}$, while the procedure to compute $\mathbf{T}^{\text{rbm}}(\mu)$ is efficient and stable.

The reduced equivalent of the variational form Equation (7) is: given $\mu \in D^\mu$, find $\mathbf{T}^{\text{rbm},N}(\mu) \in V_N$ such that:

$$a_L(\mathbf{T}^{\text{rbm},N}(\mu), \mathbf{v}; \mu) = f_L(\mathbf{v}; \mu) \quad \forall \mathbf{v} \in V_N \qquad (9)$$

The reduced space $V_N$ is constructed from snapshots, which are high fidelity solutions. The RBM consists of two main phases: (i) the offline stage, where the reduced space is constructed, and (ii) the online stage, where the reduced space is used to compute the solution of the system. The first step is performed only once and can be costly, while the second step is performed for each parameter $\mu$ and is efficient.

During the offline stage, snapshots are computed for a set of parameters $\{\mu_i\}_{i=1}^N$. This gives a family of vectors $(\mathbf{T}^{\text{fem}}(\mu_i))_{1 \le i \le N} \subset V_h$. The reduced space is defined by $V_N := \text{span}(\xi_i)$, where $(\xi_i)_{1 \le i \le N}$ is an orthonormal family of vectors, obtained by the Gram-Schmidt process applied to the snapshots $\{\mathbf{T}^{\text{fem}}(\mu_i)\}_{1 \le i \le N}$. We define the snapshots matrix $\mathbf{Z}_N = [\xi_1, \dots, \xi_N] \in \mathbb{R}^{N,N}$.

The snapshots can be selected in different ways. The first approach is to select the snapshots randomly in the parameter space, but this could lead to a poor approximation of the solution [Buf+12]. Another approach is to select the snapshots greedily, by selecting the parameter that maximizes the error between the reduced solution and the high fidelity solution, see Section 3.2.2.

Setting $\mathbf{A}_N(\mu) = \mathbf{Z}_N^T \mathbf{A}_L(\mu) \mathbf{Z}_N \in \mathbb{R}^{N \times N}$ and $\mathbf{f}_N(\mu) = \mathbf{Z}_N^T \mathbf{f}_L(\mu) \in \mathbb{R}^N$, we obtain the reduced algebraic system of size $N$:

$$\mathbf{A}_N(\mu)\mathbf{T}^{\text{rbm},N}(\mu) = \mathbf{f}_N(\mu) \qquad (10a)$$

$$s_{k,N}(\mu) = \mathbf{L}_{k,N}(\mu)^T \mathbf{T}^{\text{rbm},N}(\mu) \qquad (10b)$$

<footnote n="1">¹Source code: https://github.com/feelpp/feelpp/</footnote>

<footnote n="2">²See documentation: https://docs.feelpp.org/toolboxes/latest/heat/toolbox.html</footnote>
the same process applying for the outputs $\mathbf{L}_k$.

Thanks to the linearity of the model $\mathcal{E}_L(\mu)$, we can further write the following affine decomposition: for $T, v \in V$,

$$a_L(T, v; \mu) = \sum_{q=1}^{Q_a} \beta_A^q(\mu) a_L^q(T, v) \qquad (11a)$$

with

$$\beta_A^1(\mu) = k_{\text{lens}} \qquad a_L^1(T, v) = \int_{\Omega_{\text{lens}}} \nabla T \cdot \nabla v \, dx \qquad (11b)$$

$$\beta_A^2(\mu) = h_{\text{amb}} \qquad a_L^2(T, v) = \int_{\Gamma_{\text{amb}}} T v \, d\sigma \qquad (11c)$$

$$\beta_A^3(\mu) = h_{\text{bl}} \qquad a_L^3(T, v) = \int_{\Gamma_{\text{body}}} T v \, d\sigma \qquad (11d)$$

$$\beta_A^4(\mu) = 1 \qquad a_L^4(T, v) = \int_{\Gamma_{\text{amb}}} h_r T v \, d\sigma + \sum_{i \neq \text{lens}} k_i \int_{\Omega_i} \nabla T \cdot \nabla v \, dx \qquad (11e)$$

and

$$f_L(v; \mu) = \sum_{p=1}^{Q_f} \beta_F^p(\mu) f_L^p(v) \qquad (12a)$$

with

$$\beta_F^1(\mu) = h_{\text{amb}} T_{\text{amb}} + h_r T_{\text{amb}} - E \qquad f^1(v) = \int_{\Gamma_{\text{amb}}} v \, d\sigma \qquad (12b)$$

$$\beta_F^2(\mu) = h_{\text{bl}} T_{\text{bl}} \qquad f^2(v) = \int_{\Gamma_{\text{body}}} v \, d\sigma \qquad (12c)$$

where $Q_a = 4$ and $Q_f = 2$. We furthermore define the algebraic matrices $\mathbf{A}_L^q \in \mathbb{R}^{N \times N}$ and vectors $\mathbf{f}_L^p \in \mathbb{R}^N$, so the following equality holds:

$$\mathbf{A}_L(\mu) = \sum_{q=1}^{Q_a} \beta_A^q(\mu) \mathbf{A}_L^q, \quad \mathbf{f}_L(\mu) = \sum_{p=1}^{Q_f} \beta_F^p(\mu) \mathbf{f}_L^p \qquad (13)$$

From this decomposition and Equation (10a), we obtain the following algebraic system:

$$\mathbf{A}_N(\mu) = \sum_{q=1}^{Q_a} \beta_A^q(\mu) \underbrace{\mathbf{Z}_N^T \mathbf{A}_L^q \mathbf{Z}_N}_{\mathbf{A}_N^q} \qquad (14)$$

We set $\mathbf{A}_N^q := \mathbf{Z}_N^T \mathbf{A}_L^q \mathbf{Z}_N \in \mathbb{R}^{N \times N}$. The matrices $\mathbf{A}_N^q \in \mathbb{R}^{N \times N}$ are independent of $\mu$ and can be computed only once and stored. The same process applies to $\mathbf{f}_N(\mu)$ and $\mathbf{L}_{k,N}(\mu)$:

$$\mathbf{f}_N(\mu) = \sum_{q=1}^{Q_f} \beta_F^q(\mu) \mathbf{f}_N^q, \quad \mathbf{L}_{k,N}(\mu) = \sum_{q=1}^{Q_q} \beta_\ell^q(\mu) \mathbf{L}_{k,N}^q \qquad (15)$$

For the outputs we study in this work, the decomposition of $\mathbf{L}_{k,N}$ has only one term since the output does not depend on the parameters.

This decomposition allows implementing an offline/online procedure. During the offline phase, the basis of $V_N$ is constructed from the snapshots, as well as the matrices $\mathbf{A}_N^q$, $\mathbf{f}_N^q$, and $\mathbf{L}_{k,N}^q$ are computed and stored. More details about this construction are given in Section 3.2.2. This procedure is costly and is performed only once for the problem. During the online phase, the reduced system Equation (10) is solved for any parameter $\mu$. The entire procedure is synthesized in Algorithm 2.

During the offline stage, two approaches can be used to select the size of the reduced basis $N$: (i) an approach where we set the size of the reduced basis $N$ to a fixed value, and (ii) an approach where we set a tolerance $\epsilon_{\text{tol}}$ on the error committed on the output. The second approach is more interesting since it allows having a reduced basis of size $N$ that is adapted to the desired tolerance.
<caption>**Algorithm 2:** Offline and online stages of the RBM.</caption>

**Offline procedure:**
**Input:** Parameters $\mu_1, \dots, \mu_N \in D^\mu$
Compute snapshots $\mathbf{T}^{\text{fem}}(\mu_1), \dots, \mathbf{T}^{\text{fem}}(\mu_N)$;
Construct $\mathbf{Z}_N \leftarrow [\xi_1, \dots, \xi_N]$ (orthonormal);
Construct the reduced matrices $(\mathbf{A}_N^q)_{1 \le q \le Q_a}, (\mathbf{F}_N^p)_{1 \le p \le Q_f}, (\mathbf{L}_{k,N})_{1 \le k \le n_{\text{output}}}$;
**Output:** Reduced basis and reduced matrices, stored.

**Online procedure:**
**Input:** $\mu \in D^\mu$
Assemble $\mathbf{A}_N(\mu), \mathbf{F}_N(\mu), \mathbf{L}_{k,N}(\mu)$ using the saved matrices and the affine decomposition;
Solve $\mathbf{A}_N(\mu)\mathbf{u}_N(\mu) = \mathbf{F}_N(\mu)$;
Compute the output $s_{k,N}(\mu) = \mathbf{L}_{k,N}(\mu)^T\mathbf{u}_N(\mu)$;
**Output:** $\mathbf{u}_N(\mu), s_{k,N}(\mu)$.

### 3.2.1 Error estimates

**Error bound** Given the reduced basis approximation $\mathbf{T}^{\text{rbm},N}(\mu)$ of the high fidelity resolution solution $\mathbf{T}^{\text{fem}}(\mu)$ for $\mu \in D^\mu$, we defined the field error

$$ \mathbf{e}(\mu) := \mathbf{T}^{\text{fem}}(\mu) - \mathbf{T}^{\text{rbm},N}(\mu) \qquad (16) $$

We want to construct quantities $\Delta_N(\mu)$ and $\Delta_N^s(\mu)$ such that

$$ \|\mathbf{e}(\mu)\|_V \le \Delta_N(\mu) \quad \text{and} \quad s(\mu) - s_N(\mu) \le \Delta_N^s(\mu) \qquad (17) $$

Those quantities are named a *posterior error bounds* [Pru+01; RHP07]. To quantify the sharpness and *rigor* properties of the error bound, we introduce the *effectivity*:

$$ \eta_N(\mu) := \frac{\Delta_N(\mu)}{\|\mathbf{e}(\mu)\|_V} \qquad \eta_N^s(\mu) := \frac{\Delta_N^s(\mu)}{s(\mu) - s_N(\mu)} \qquad (18) $$

It has been proven in [Pru+01] that the error bound is *rigorous*, namely that it is always greater than the error; and *sharp*, namely that it is as close as possible to the actual error.

These properties can be summarized as:

$$ 1 \le \eta_N(\mu) \le \eta_{ub}(\mu) \quad \forall \mu \in D^\mu \qquad (19) $$

where $\eta_{ub}(\mu)$ is a the sharpness of the bound and is proven to be bounded [Pru+01] when $N$ increases.

Finally, to construct the reduced space, we require the error bound to be *efficient*, that is, its evaluation is independent of the size of the high fidelity space $\mathcal{N}$. This is critical when heuristic algorithms are used to construct the reduced space, such as the Greedy algorithm discussed in Section 3.2.2.

Such an error bound can be constructed efficiently from the *residual* $r$ of the problem (5), a lower bound $\alpha_{\text{lb}}(\mu)$ of the coercivity constant $\alpha(\mu)$ of $a_L(\cdot, \cdot; \mu)$, and the affine decomposition of $a_L$ and $f_L$:

$$ \alpha(\mu) = \inf_{v \in V} \frac{a(v, v; \mu)}{\|v\|_V^2}, \quad r(v, \mu) := \ell(v; \mu) - a(u_N(\mu), v; \mu) \quad \forall v \in V \quad \text{and} \quad \Delta_N := \frac{\|r(\cdot, \mu)\|_{V'}}{\alpha_{\text{lb}}(\mu)} \qquad (20) $$

For more details, refer to [Pru+01].

### 3.2.2 Generation of the reduced basis

The a posteriori error estimator introduced earlier provides an efficient criterion to select the desired dimension of the reduced space $N$, in the offline phase. Given a fixed tolerance $\varepsilon_{\text{tol}}$, we can greedily select the greatest $N$ such that the error bound is smaller than the tolerance. In this section, we describe an algorithm to generate the reduced basis.

For this algorithm, a large set of parameters $\Xi_{\text{train}} \subset D^\mu$ is required. This set is called the *training set*, and is generated log-randomly. A first snapshot is computed for a given parameter $\mu_0 \in D^\mu$. To get the $N+1$-th snapshot to be inserted in the basis, we select the parameter $\mu^*$ that maximizes the error bound $\Delta_N(\mu)$, for $\mu \in \Xi_{\text{train}}$. This step is performed until a selected tolerance for the maximal error bound is reached. The greedy algorithm is presented in Algorithm 3.
**Algorithm 3:** Greedy algorithm to construct the reduced basis.

Input: $\mu_0 \in D^\mu$, $\Xi_{\text{train}} \subset D^\mu$ and $\varepsilon_{\text{tol}} > 0$

S $\leftarrow$ [$\mu_0$];

while $\Delta_N^{\max} > \varepsilon_{\text{tol}}$ do

$\mu^* \leftarrow \arg \max_{\mu \in \Xi_{\text{train}}} \Delta_N(\mu)$ (and $\Delta_N^{\max} \leftarrow \max_{\mu \in \Xi_{\text{train}}} \Delta_N(\mu)$);

$V_{N+1} \leftarrow \{\boldsymbol{\xi} = \mathbf{T}^{\text{fem}}(\mu^*)\} \cup V_N$;

Append $\mu^*$ to S;

$N \leftarrow N + 1$;

end

Output: Sample S, reduced basis $V_N$

<table><caption>Table 2: Characteristics of the meshes.</caption><thead><tr><th>Mesh</th><th>h</th><th>nDof P<sub>1</sub></th><th>nDof P<sub>2</sub></th></tr></thead><tbody><tr><td>M0</td><td>0.86</td><td>47,284</td><td>3.27 &middot; 10<sup>5</sup></td></tr><tr><td>M1</td><td>0.74</td><td>68,993</td><td>4.73 &middot; 10<sup>5</sup></td></tr><tr><td>M2</td><td>0.51</td><td>1.21 &middot; 10<sup>5</sup></td><td>8.83 &middot; 10<sup>5</sup></td></tr><tr><td>M3</td><td>0.47</td><td>2.08 &middot; 10<sup>5</sup></td><td>1.58 &middot; 10<sup>6</sup></td></tr><tr><td>M4</td><td>0.29</td><td>9.96 &middot; 10<sup>5</sup></td><td>7.87 &middot; 10<sup>6</sup></td></tr><tr><td>M5</td><td>0.15</td><td>7.36 &middot; 10<sup>6</sup></td><td>5.87 &middot; 10<sup>7</sup></td></tr></tbody></table>

## 3.3 Numerical results

### 3.3.1 Mesh convergence

In Section 2.1, we detailed the geometry of the eyeball, derived from computer-aided design (CAD) data. As illustrated in Figure 1, certain regions exhibit greater complexity than others. For instance, the lamina cribrosa is notably thinner, while the iris presents a less uniform structure. Achieving an effective mesh requires a well-distributed arrangement of elements. This is attainable through the application of a specialized meshing algorithm designed to tailor the mesh according to the geometric intricacies. Utilizing the MMG library [too22], we have generated a family of meshes with varying levels of refinement. These meshes are used to our subsequent simulation processes.

Table 2 displays the characteristics of the meshes, such as their characteristic size $h$ and the number of degrees of freedom (nDof) for both $\mathbb{P}_1$ and $\mathbb{P}_2$ finite element discretizations.

In this section, we detail the outcomes of our mesh convergence analysis. This analysis involves solving the given problem across various mesh configurations and subsequently comparing the resultant data. To conduct this study, we first solve the model denoted as $\mathcal{E}_{\text{NL}}(\bar{\mu})$. Following this, we compute the output $T_{\text{NL,O}}^{\text{fem}}$ representing the temperature at the cornea's center as determined by the high-fidelity model. The primary objective of this analysis is to ascertain whether the obtained temperature values demonstrate convergence towards a consistent value. Figure 4 illustrates the results of our mesh convergence study, clearly indicating a pattern of satisfactory convergence. We select for further comparisons the values obtained for M3 and for $\mathbb{P}_2$.

### 3.3.2 Scalability

In this section, we explore the scalability of our computational framework. This involves measuring the time required to solve the model in relation to the number of MPI parallel processes utilized. The time measured pertains to the duration necessary for assembling the algebraic system and solving the problem, as per Algorithm 1. Our experiments utilized mesh M3, with both $\mathbb{P}_1$ and $\mathbb{P}_2$ discretizations. The results, presented in Figure 5, demonstrate satisfactory scalability: the execution time decreases as the number of parallel processes increases. However, we observed that beyond 12 processes, the reduction in execution time becomes less significant. Consequently, for optimal efficiency, we have selected 12 processes for our subsequent analyses. This study sets the stage for a subsequent comparison with a reduced-order
<figure><img src="image_6.png" /><figcaption>Figure 4: Temperature at the center of the cornea computed with the high-fidelity model $\mathcal{E}_{\text{NL}}(\bar{\mu})$, depending on the level of refinement of the mesh.</figcaption></figure>

model, which employs a reduced basis with reliable, certified output bounds derived from the high-fidelity solutions. This comparison aims to highlight that, while parallel computing can accelerate the high-fidelity computation, the reduced-order approach offers even more substantial computational gains.

### 3.3.3 Linearized model

We now compare the results obtained after solving $\mathcal{E}_{\text{NL}}(\bar{\mu})$ against $\mathcal{E}_{\text{L}}(\bar{\mu})$. We denote the solution of the nonlinear model $\mathcal{E}_{\text{NL}}(\bar{\mu})$ by $T_{\text{NL}}(\bar{\mu})$, and by $T_{\text{L}}(\bar{\mu})$ the solution of the linearized model $\mathcal{E}_{\text{L}}(\bar{\mu})$, and compute the relative error:

$$e_{\text{lin}}(\mu) = \frac{\|T_{\text{NL}}(\mu) - T_{\text{L}}(\mu)\|_{L^2(\Omega)}}{\|T_{\text{NL}}(\mu)\|_{L^2(\Omega)}} \quad (21)$$

For $\mu = \bar{\mu}$, we get $e_{\text{lin}} = 4.176 \cdot 10^{-7}$. In Figure 6, we plot the difference between the two solutions $|T_{\text{NL}}(\mathbf{x}) - T_{\text{L}}(\mathbf{x})|$ for $\mathbf{x} \in \Omega$. We notice that the difference is the largest on the front of the eye, where the boundary condition has been changed. Whereas at the back of the eye, the solutions are superposed. We also compute the maximal difference: $2.1667 \cdot 10^{-3}$ K. Therefore, we consider in the sequel that the linearized model does not induce a significant error in the results.

Figure 7 displays the results of the simulation of the linear model $\mathcal{E}_{\text{L}}(\mu)$, for three parameters $\mu$: $\bar{\mu}$ the baseline value parameters, $\mu_{\text{min}}$ (resp. $\mu_{\text{max}}$) where each component is the lowest (resp. highest) bound of its range of values.

### 3.3.4 Verifications of the reduced basis model

We compare the results of the reduced basis method with the output of the high fidelity FEM model. We generate a sample $\Xi_{\text{test}}$ of 100 parameters in $D^\mu$. For $\mu \in \Xi_{\text{test}}$, we compute on the one hand $T_{\text{O}}^{\text{fem}}(\mu)$, the value of the temperature at point O from the model $\mathcal{E}_{\text{L}}(\mu)$, and on the other hand $T_{\text{O}}^{\text{rbm},N}(\mu)$, the value of the temperature for the reduced basis model, with a basis of size N. In Figure 8(a), the value of the error $|T_{\text{O}}^{\text{fem}}(\mu) - T_{\text{O}}^{\text{rbm},N}(\mu)|$ is plotted for each $\mu \in \Xi_{\text{test}}$, for various reduced basis sizes N. Statistics on the error committed over the sample $\Xi_{\text{test}}$ are displayed in Figure 8(b), as well as the effectivity $\eta_N(\mu)$ in Figure 8(c).

We observe that even for small values of N, the error on the output is remarkably small: an error of $10^{-4}$ is reached for $N = 6$. On the other hand, we find that the convergence on the output is twice as fast as the convergence on the field, as predicted by the theoretical error estimate [Pru+01, Eq. (36)].
<figure><img src="image_7.png" /><figcaption>Figure 5: Time of execution to run $\mathcal{E}_L^N(\bar{\mu})$ and corresponding speed-up, for an increasing number of parallel processes. Simulations are performed on the M3.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Difference of the temperature between the full model and the linearized model, computed on the mesh M3 with $P_2$ elements, and the baseline values $\bar{\mu}$ for the parameters.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Distribution of the temperature [K] in the eyeball from the linear model $\mathcal{E}_L(\mu)$.</figcaption></figure>

<table><caption>Table 3: Times of execution, using mesh M3 for high fidelity simulations.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Finite element resolution</th><th>Reduced model</th></tr><tr><th>P<sub>1</sub></th><th>P<sub>2</sub> (np=1)</th><th>P<sub>2</sub> (np=12)</th><th>T<sup>rbm,N</sup>(μ), Δ<sub>N</sub>(μ)</th></tr></thead><tbody><tr><td>Problem size</td><td>N = 207 845</td><td>N = 1 580 932</td><td></td><td>N = 10</td></tr><tr><td>t<sub>exec</sub></td><td>5.534 s</td><td>62.432 s</td><td>10.76 s</td><td>2.88 × 10<sup>-4</sup> s</td></tr><tr><td>speed-up</td><td>11.69</td><td>1</td><td>5.80</td><td><b>2.17 × 10<sup>5</sup></b></td></tr></tbody></table>

Note that the anticipated error behavior aligns with theory when the output functional maintains continuity. In this work, we deviate from the standard case, attributed to the utilization of the Dirac functional in output computation used for pointwise evaluation. Nevertheless, a similar behavior is observed, and additional insights into this phenomenon will be provided in future work.

Table 3 offers a comparative analysis of execution times for solving the heat transfer problem. We first discuss the execution times for the high-fidelity solution, encompassing both $P_1$ and $P_2$ finite-element discretizations. The measured time, denoted as $t_{exec}$, includes assembling and solving the problem. In contrast, we also evaluate the execution time of the online phase of our certified reliable reduced basis model. This comparison highlights a significant reduction in the time required to assemble and solve the problem using our advanced reduced basis approach. Importantly, this efficiency does not compromise accuracy; the results from the reduced basis model are effectively exact with respect to the high fidelity model. As anticipated in our earlier scalability analysis, we achieve remarkable computational gains with our model, reinforcing the benefits of our approach in both precision and performance.

The reduced bases constructed for the various outputs of interest are generated with the greedy algorithm (Algorithm 3), using a maximal tolerance $\epsilon_{tol} = 1 \cdot 10^{-6}$ and a maximal size for the basis $N = 20$. In practice, the tolerance is reached for $N = 10$ to 12.

## 3.4 Validation and comparison with previous studies

We present in this section a thorough comparison between the results of this work and previously published data on the temperature of the eye, obtained either by experimental procedures or via computational modeling. Note that only scarce data are available for the entire human eyeball, since most of the measurement techniques estimated only the surface temperature of the cornea. In particular, [EYB89] gathers the outputs of 19 studies conducted with various instruments (mercury bulbs, liquid crystal thermometers or infrared thermometers), and the mean value reported, according to [Ng 06], is $T_O^{exp} = 307.15$ K. The temperature at the center of the cornea computed with baseline value from our model is $T_O^{fem}(\bar{\mu}) = 306.02$ K, which lies in the interval of results found in the literature (see [EYB89,
<figure><img src="image_10.png" /><figcaption>(a) Error on RBM for various reduced basis sizes with error bound $\Delta_N^s(\mu)$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Convergence of the errors on the field and the output on (c) Stability of the effectivity $\eta_N^s(\mu) = \frac{\Delta_N^s(\mu)}{|T_O^{\text{fem}}(\mu) - T_O^{\text{rbm,N}}(\mu)|}$ point O, for $\mu \in \Xi_{\text{test}}$, for various reduced basis sizes. The maximal, minimal, and mean values are represented. and $\eta_N^{\text{en}}(\mu) = \frac{\Delta_N^{\text{en}}(\mu)}{\|T^{\text{fem}}(\mu) - T^{\text{rbm,N}}(\mu)\|_{\mu}}$ for $\mu \in \Xi_{\text{test}}$, for various reduced basis sizes. The full red line represents the theoretical lower bound of the effectivity, i.e. 1.</figcaption></figure>

Figure 8: Comparison of the temperature between the full order model and the reduced basis model, tested over a sample $\Xi_{\text{test}} \subset D^\mu$ of 100 parameters.
<figure><img src="image_12.png" /><figcaption>Figure 9: Temperature on the GCC: experimental data (mean and standard deviation) vs. numerical results. In this analysis, we cannot distinguish graphical difference between the linear and non-linear models.</figcaption></figure>

Table 1] and [Ng 06, Table 9]).

Additionally, in [EYB89], the temperature is measured along an imaginary horizontal line, the *Geometrical Center of the Cornea* (GCC), as described in Figure 1, on a panel of 21 subjects. The experimental data are displayed in Figure 9, together with the findings of the present work. On the horizontal axis, the distance to the center of the eye is represented, and on the vertical axis is the temperature difference to the central one (mean value and standard deviation). Note that as the geometry of the simulated eye is not the same as the one used in the experiment, we scaled the results over the x-axis. The result shows that the high fidelity model is able to closely replicate the same behavior as the one experimentally measured, and the model $E_L(\bar{\mu})$ provides very close values (see Section 3.3.3). Moreover, thanks to the error bound introduced in Section 3.2.1 for the RBM, the approach is considered to be valid for the sensitivity analysis procedure hereafter.

In Figure 10, we present a comparative analysis between the results of our current study and various numerical findings reported in existing literature. This comparison features temperatures calculated along a line traversing the eye's center, the specific location of which is depicted in Figure 2. This comparative approach is crucial as it verifies the accuracy of our computed values, encompassing not just the corneal surface but also the eye's internal tissue structures. It is noteworthy that our analysis includes a mix of both 2D and 3D results, derived from both non-linear and linearized models.

The results show a very good agreement between the findings of the present study and previously reported temperature results, along the different locations in the eyeball.

# 4 Uncertainty quantification

The *uncertainty quantification* (UQ) allows quantifying the uncertainty of the model parameters on various outputs. In the present work, we focus on forward UQ, that is, we want to quantify the uncertainty and the sensitivity of the output of the model, given the uncertainty of the input parameters. More precisely, two studies are performed: (i) an uncertainty propagation, to understand how the uncertainties of the inputs of the model are propagated to the output via the computational model, and (ii) a *sensitivity analysis* (SA) to assess the impact of varying selected parameters on several outputs of interest, namely temperature at specific locations in the eye. Their locations are detailed in Figure 2.

The SA is conducted in two different approaches. First, to recapitulate findings from the literature, we performed a deterministic SA, where for each simulation, only one parameter is allowed to vary in a given range, whereas the others are fixed to their baseline value. In a second stage, we extended the SA to a stochastic framework, where each selected parameter follows a given random distribution and
<figure><img src="image_13.png" /><figcaption>Figure 10: Temperature on a line going through the center of the eye (see Figure 2), comparison with numerical results from literature.</figcaption></figure>

the impact on the quantity of interest is assessed via sensitivity indices. The advantage of the latter is the global perspective provided by this method and its ability to capture high-order interactions among several input parameters.

## 4.1 Deterministic sensitivity analysis

Our initial investigation of the impact of varying selected parameters is conducted through a *deterministic sensitivity analysis* (DSA). Specifically, we choose a parameter among the ones defined in Table 1, and we set the other to their baseline value. Next, we vary the selected parameter among pre-defined values and compute the outputs of the high fidelity model. Similar studies were performed in [Sco88; Ng 06; Li+10]. We gather in the present study information about several parameters of interest for the heat transfer model from these studies, namely baseline values and ranges. These variations correspond not only to physiological conditions but also include some extreme situations. We postpone a more in-depth discussion on this topic to Section 4.2, where the random distributions characterizing these parameters are set up. Note that in this case, we do not need to use the reduced model, since only a relatively small number of simulations is required.

In Figure 11, we show the results of the DSA for the parameter $\mu = \{E, h_{amb}, h_{bl}, k_{lens}, T_{amb}, T_{bl}\}$, on point O which is at the surface of the cornea. The plain-line curves correspond to the results reported in the literature, which we compare with the results of our simulations; the vertical dashed line corresponds to the baseline value for each parameter. The results are in very good agreement with previous findings and show that temperature at the level of the cornea is strongly influenced by $h_{amb}$, $T_{amb}$, $E$, and $T_{bl}$, whereas the influence of $h_{bl}$ and $k_{lens}$ is less significant. For instance, high air conductivity can result in a temperature 7 K lower than the baseline value, while the difference obtained for $h_{bl}$ and $k_{lens}$ in the computed temperature is at most of 1 K.

## 4.2 Stochastic sensitivity analysis (SSA)

We consider an output quantity Y depending on a set of input parameters $\mu \in D^\mu$, and we estimate the sensitivity of Y to each parameter $\mu_i$ for $i \in [1, d]$, where d is the dimension of the parametric space. To this end, we compute the *Sobol' sensitivity indices* introduced in [Sob93] as follows. We assume that each component $\mu_i$ of $\mu$ follows a random variable $X_i$, independent of the others. The first-order indices are defined as:

$$S_i := \frac{\mathrm{var}(\mathbb{E}[Y|X_i])}{\mathrm{var}(Y)} \qquad (22)$$
<figure><img src="image_14.png" /><figcaption>(a) Evaporation rate $E$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Ambient convective coefficient $h_{amb}$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) Blood convection coefficient $h_{bl}$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) Conductivity of the lens $k_{lens}$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(e) Ambient temperature $T_{amb}$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(f) Blood temperature $T_{bl}$</figcaption></figure>

Figure 11: Results of the DSA for the 6 parameters studied, among previous studies from the literature (markers x $E_{NL}(\mu)$, □ [Ng 06], + [Sco88], * [Li+10]). The vertical dashed line corresponds to the baseline value of the parameter.
where var(Y) corresponds to the variance of Y including the eventual non-linearity effect of the coefficient on the output, and var (E [Y | X$_{i}$]) is the variance of the conditional expectation of Y given X$_{i}$, corresponding to the first order effect of the parameter µ$_{i}$ on the output: if the parameter modeled by the distribution X$_{j}$ has a great impact on the output Y, then E [Y | X$_{j}$] will vary has well, and so its variance.

We also define the total Sobol’ index:

$$
S_i^{\text{tot}} := \frac{\operatorname{var}(\mathbb{E}[Y|X_{(-i)}])}{\operatorname{var}(Y)} = 1 - S_{-i} \quad (23)
$$

where $X_{(-i)} = (X_1, \cdots, X_{i-1}, X_{i+1}, \cdots, X_d)$ is the set of parameters without the parameter $X_i$, and $S_{-i}$ is the sum of the indices where $X_i$ is not present.

To compute the Sobol' indices, we use an algorithm of functional chaos, implemented in the library OpenTURNS [Bau+16] by the class FunctionalChaosAlgorithm, using a bootstrap method<fnref n="3" /> for the confidence intervals.

4.2.1 Choice of the distributions

We discuss now the prior distributions for each parameter. Each parameter does not depend on the
others, resulting in a family of 6 random independent variables. Figure 12 shows the probability density
function (PDF) of distributions of the parameters, associated with the baseline values (see Table 1),
where the parameters used in the literature for the deterministic sensitivity analysis are represented with
a vertical line. We present hereafter some details on how the random distributions were constructed.

*   **Evaporation rate E:** According to [Sco88], the evaporation rate's range is of 40 to 100 W m⁻², using data from literature [Adl53]. The value E = 40 W m⁻² is chosen as the baseline value. Some high values are also considered, to study the impact of important evaporation rates. The values used in the literature run from 20 to 320 W m⁻². As this parameter varies by several orders of magnitude, we decided to use a log-normal distribution to represent it. More precisely we set E ~ log-N(µE, σE, γE), with σE = 0.7, µE = log(40) - 0.15²/2 and γE = 20, restricted to [20, 130]. The distribution is presented in Figure 12(a). This choice of the distribution leads to a mean value of Ē = 55.8 W m⁻².

*   **Ambient air convection coefficient h_amb:** In [Sco88], the sole value given for the ambient air convection coefficient is 10 W m⁻² K⁻¹, and similar values are used to run the DSA, from 8 to 15 W m⁻² K⁻¹. Other results in the literature coroborate this value: [Kos+13, Table 12.2] reports a range of 2.5 to 25 W m⁻² K⁻¹ for a free convection, and 10 to 500 W m⁻² K⁻¹ for a forced convection. [Edg] proposes a range of 10 to 100 W m⁻² K⁻¹ for the air. In their DSA, [Ng 06] and [Li+10] use higher values of h_amb, up to 100 W m⁻² K⁻¹ to simulate a forced convection condition. As high values are not a common case, such a coefficient should not have a high frequency in the distribution. We chose to use a log-normal distribution: h_amb ~ log-N(log(10) - 1/2, 1, 8). In Remark 2.1, we discussed the linearization process of the model, inducing the usage of a fixed parameter h_r chosen to fit temperature in usual ambient room conditions, which leads to a restriction of the distribution to the interval [8, 100] W m⁻² K⁻¹. The distribution is presented in Figure 12(b). The mean value of the distribution is h̄_amb = 17.6 W m⁻² K⁻¹.

*   **Ambient air convection coefficient $h_{amb}$:** In [Sco88], the sole value given for the ambient air convection coefficient is $10 \text{ W m}^{-2} \text{ K}^{-1}$, and similar values are used to run the DSA, from 8 to 15 $\text{W m}^{-2} \text{ K}^{-1}$. Other results in the literature coroborate this value: [Kos+13, Table 12.2] reports a range of 2.5 to 25 $\text{W m}^{-2} \text{ K}^{-1}$ for a free convection, and 10 to 500 $\text{W m}^{-2} \text{ K}^{-1}$ for a forced convection. [Edg] proposes a range of 10 to 100 $\text{W m}^{-2} \text{ K}^{-1}$ for the air. In their DSA, [Ng 06] and [Li+10] use higher values of $h_{amb}$, up to 100 $\text{W m}^{-2} \text{ K}^{-1}$ to simulate a forced convection condition. As high values are not a common case, such a coefficient should not have a high frequency in the distribution. We chose to use a log-normal distribution: $h_{amb} \sim \log-N(\log(10) - \frac{1}{2}, 1, 8)$. In Remark 2.1, we discussed the linearization process of the model, inducing the usage of a fixed parameter $h_r$ chosen to fit temperature in usual ambient room conditions, which leads to a restriction of the distribution to the interval [8, 100] $\text{W m}^{-2} \text{ K}^{-1}$. The distribution is presented in Figure 12(b). The mean value of the distribution is $\bar{h}_{amb} = 17.6 \text{ W m}^{-2} \text{ K}^{-1}$.

*   **Blood convection coefficient $h_{bl}$:** A control value of 65 $\text{W m}^{-2} \text{ K}^{-1}$, derived from experimental data [Lag82] is provided in [Sco88]. For the DSA, the values used run from 50 to 120 $\text{W m}^{-2} \text{ K}^{-1}$. This leads us to the following assumption for the distribution of the parameter: $h_{bl} \sim \log-N(\log(65) - \frac{0.15^2}{2}, 0.15, 0)$, restricted over [50, 120], see Figure 12(c). The mean of this distribution is $\bar{h}_{bl} = 65.8 \text{ W m}^{-2} \text{ K}^{-1}$.

*   **Lens conductivity $k_{lens}$:** This parameter is chosen among all the conductivities since the water content of the lens varies with aging [Sco88]. [Sco88] and [Ng 06] run the DSA with this parameter, using values from 0.21 to 0.544 $\text{Wm}^{-1}\text{K}^{-1}$. As the range of values is not very large, it seems reasonable to use a uniform distribution for this parameter: $k_{lens} \sim U(0.21, 0.544)$, see Figure 12(d).

*   **Blood convection coefficient $h_{bl}$:** A control value of 65 W m⁻² K⁻¹, derived from experimental data [Lag82] is provided in [Sco88]. For the DSA, the values used run from 50 to 120 W m⁻² K⁻¹. This leads us to the following assumption for the distribution of the parameter: $h_{bl} \sim \log-N(\log(65) - \frac{0.15^2}{2}, 0.15, 0)$, restricted over [50, 120], see Figure 12(c). The mean of this distribution is $\bar{h}_{bl} = 65.8 \text{ W m}^{-2} \text{ K}^{-1}$.

*   **Lens conductivity $k_{lens}$:** This parameter is chosen among all the conductivities since the water content of the lens varies with aging [Sco88]. [Sco88] and [Ng 06] run the DSA with this parameter, using values from 0.21 to 0.544 Wm⁻¹K⁻¹. As the range of values is not very large, it seems reasonable to use a uniform distribution for this parameter: $k_{lens} \sim U(0.21, 0.544)$, see Figure 12(d).

<footnote n="3">³See documentation https://openturns.github.io/openturns/latest/auto_meta_modeling/p polynomial_chaos_m etam
odel/plot_chaos_sobol_confidence.html</footnote>
<figure><img src="image_20.png" /><figcaption>(a) Distribution of $E$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Distribution of $h_{amb}$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) Distribution of $h_{bl}$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(d) Distribution of $k_{lens}$</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(e) Distribution of $T_{amb}$</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(f) Distribution of $T_{bl}$</figcaption></figure>

Figure 12: Distributions of the parameters. The vertical lines represent the values chosen in literature for the DSA.

* **Ambient temperature** $T_{amb}$: The baseline value of this parameter is taken to a usual room temperature of 294 K (20 °C). The values used for the DSA vary from extreme conditions of 273 K (0 °C) to 308 K (35 °C). As these extreme values are not very common, we choose to restrict the values taken by $T_{amb}$ from 283.15 K (−35 °C) to 303.15 K (30 °C): $T_{amb} \sim U(283.15, 303.15)$. The distribution is presented in Figure 12(e).

* **Blood temperature** $T_{bl}$: The temperature of human blood is commonly accepted to be 310 K (37 °C). For the DSA, cases of hypothermia and hyperthermia are considered, with a range from 308 K to 312.15 K (35 °C to 39 °C). We therefore take $T_{bl} \sim U(308, 312.15)$, see Figure 12(f).

### 4.2.2 Uncertainty propagation

We focus on the distribution of outputs of interest, from a random sample of the input parameters of get from the distributions presented 10,000 points, which leads to a number of 10,000 simulations. The computational cost of the high fidelity simulations becomes in this case prohibitive and therefore we employ the reduced basis metamodel developed in Section 3.2. Figure 13 presents the distribution of three outputs, namely the mean of the temperature over the cornea $T_{cornea}$, and the temperature on points *O* and *G* respectively are the front and the back of the eyeball. Note that $T_O$ and $T_{cornea}$ display a Gaussian distribution, whereas $T_G$ is more difficult to interpret, but could correspond to a uniform or bi-modal distribution.

We provide in Table 4 results about mean values and standard deviation for the same quantities. We note that the mean values of $T_O$ and $T_{cornea}$ are of the same order of magnitude as the experimental data in the validation section (Section 3.4): the difference of temperature is about 2 K, and standard deviations are in the same ranges. The mean value of $T_G$ is very close to results reported in Figure 10 from the literature with a small standard deviation.
<table><caption>Table 4: Statistics of the outputs.</caption><thead><tr><th></th><th><i>T</i><sub>cornea</sub></th><th><i>T</i><sub>O</sub></th><th><i>T</i><sub>G</sub></th></tr></thead><tbody><tr><td>Mean</td><td>305.590082</td><td>303.185187</td><td>310.028526</td></tr><tr><td>Standard deviation</td><td>1.788358</td><td>2.457063</td><td>1.055978</td></tr></tbody></table>

<figure><img src="image_26.png" /><figcaption>Figure 13: Distribution of the output, from the composed input distribution.</figcaption></figure>

### 4.2.3 Results of the SA

We perform a sensitivity analysis to compute Sobol's indices, a convergence analysis varying the sampling size $N_{\text{param}}$. Table 5 reports the maximal deviation of these indices, and the time taken by the application to compute the 6 sets of Sobol's indices.

Additionally, we compute the predictivity factor $Q_2$ for the polynomial chaos metamodel is defined as:

$$Q_2 := 1 - \frac{\sum_{l=1}^{N} (Y_l - \hat{f}(X_l))^2}{\text{var}(Y)}, \qquad (24)$$

measuring how accurate the metamodel $\hat{f}$ is at predicting the output $Y$ from the input $X$. The closer $Q_2$ is to 1, the better the metamodel is. In the context of the Sobol's indices experiment, the metamodel $\hat{f}$ is the polynomial chaos expansion of the output $Y$. The test of convergence is performed using the temperature on point $O$ as the output. The convergence of Sobol's indices is reached for $N_{\text{param}} = 200$ with a $10^{-2}$ accuracy, which is a threshold used in the sequel.

Figure 14 shows the results of the Sobol analysis for different outputs of interest. Recall that Figure 2 shows where the points are in the eye.

<table><caption>Table 5: Convergence of the Sobol's indices.</caption><thead><tr><th><i>N</i><sub>param</sub></th><th>Max deviation</th><th><i>t</i><sub>exec</sub></th><th><i>Q</i><sub>2</sub></th></tr></thead><tbody><tr><td>60</td><td>0.18102</td><td>0.75609 s</td><td>0.999153</td></tr><tr><td>100</td><td>0.03698</td><td>1.99651 s</td><td>0.992648</td></tr><tr><td>150</td><td>0.02969</td><td>2.83743 s</td><td>0.99986</td></tr><tr><td>200</td><td>0.02923</td><td>4.16046 s</td><td>0.998926</td></tr><tr><td>400</td><td>0.00739</td><td>8.36701 s</td><td>0.999931</td></tr><tr><td>600</td><td>0.00496</td><td>15.7947 s</td><td>0.9998</td></tr><tr><td>1000</td><td>0.00248</td><td>22.364 s</td><td>0.999904</td></tr></tbody></table>
In the deterministic sensitivity analysis conducted in Section 4.1, the impact of the variation of a
sole parameter on the temperature at point *O* was studied. Using Sobol's indices, we are now able to
measure the impact when all of them are varying. The results of Sobol analysis at point *O* presented in
Figure 14(a) are in very good agreement with the deterministic findings: the temperature at the level of
the cornea is strongly influenced by external factors such as *h*<sub>amb</sub>, as well subject-specific parameters such
as *T*<sub>amb</sub>, *E*, and *T*<sub>bl</sub>. Moreover, it is minimally influenced by the lens conductivity *k*<sub>lens</sub> and the blood
convection coefficient *h*<sub>bl</sub>.

Sobol's indices for several other locations are gathered in Figure 14(b–f). From these results, we
can infer the following ranking of the influential parameters: T<sub>amb</sub>, h<sub>amb</sub>, E, and T<sub>bl</sub>. In particular, the
dependence of the ambient temperature T<sub>amb</sub> decreases when we go deeper inside the eye. Precisely, the
impact of T<sub>amb</sub> is still significant for the mean temperature of the cornea, but the other parameters are
equally influential. These behaviors are coherent with physiological conditions. Moreover, regardless of
the output studied, the parameters k<sub>lens</sub> and h<sub>bl</sub> are minimally influencing the output. Consequently in
future simulations, their value can be set at baseline. Surprisingly, the temperature at B<sub>1</sub>, on the lens, is
minimally influenced by k<sub>lens</sub>, but this parameter has a minimal role in the modeling process. On the
other hand, T<sub>bl</sub> is very influential at D<sub>1</sub> and G, close to vascular beds, again in a coherent manner with
the physiological situation. Finally, we can notice a slight difference between the first-order and total-order
indices, mostly for h<sub>amb</sub> and T<sub>amb</sub>, implying that there are high-order interactions among these selected
parameters. To measure the impact of coupled parameters, second-order Sobol's indices computation is
required, but the polynomial chaos expansion does not directly provide these values. Alternatively, a
Monte-Carlo based method could be implemented which is very costly to the computational viewpoint.

5 Conclusion

We have successfully developed a numerical model that accurately simulates heat transfer within
the complex three-dimensional structure of the human eyeball, enabling us to calculate the temperature
distribution across various ocular tissues. This model has undergone rigorous validation against both
experimental data and numerical results from existing literature. A key advancement in our study is the
implementation of a certified reduced basis method. This method significantly accelerates the simulations
of our complex model while maintaining high accuracy, making it highly efficient for many-queries
computations essential in uncertainty quantification studies. Our sensitivity analysis pinpointed four
main physiological parameters as most influential in affecting the results: blood temperature, ambient
temperature, the ambient air convection coefficient, and the evaporation rate. These findings build
upon and enrich prior studies, such as those highlighted in [Sco88; Ng 06; Li+10], underscoring the
vital role of blood flow characteristics and environmental conditions, particularly in the inner ocular
tissues. Additionally, through Sobol' indices analysis, we identified the significant impact of parameter
interactions, particularly those related to ambient temperature. From a clinical standpoint, our insights
into heat transport in the human eye could inform studies on the effects of electromagnetic wave radiation,
as explored in [Hir+07; NO07; ON09] and related references. As a next step, we plan to couple the heat
transfer model with models describing the aqueous humor flow, as in [ON08; Abd+21]. This work is
crucial to understand the formation of Krukenberg's spindle, to better understand the pathophysiology of
pigment dispersion syndrome and to enhance drug delivery efficiency in the anterior chamber [Wan+16;
BBS20]. Our preliminary findings in this direction can be found in [Sai+24]. Ultimately, our work, in
conjunction with previous initiatives such as the Ocular Mathematical Virtual Simulator [Sal+], lays
the groundwork for a comprehensive, multi-physics, multiscale framework in ophthalmology, tailored for
personalized medical applications.

A Reproductibility of results

All the codes used to run the simulations are available on the Feel++ [Pru+23] GitHub repository^4).
Details are given in the README.md file of the repository: https://github.com/feelpp/feelpp/tree/
develop/mor/examples/eye2brain

<footnote n="4">⁴https://github.com/feelpp/feelpp</footnote>
<figure><img src="image_27.png" /><figcaption>Figure 14: Sobol's indices for the SSA.</figcaption></figure>
# Acknowledgements

The authors would like to acknowledge the support of the platform Cemosis at University of Strasbourg and the French Ministry of Higher Education, Research and Innovation. This work has benefited from a national grant managed by the French National Research Agency (Agence Nationale de la Recherche) attributed to the Exa-MA project of the NumPEx PEPR program, under the reference ANR-22-EXNU-0002.

# References

[Abd+21] Farah Abdelhafid et al. "Operator Splitting for the Simulation of Aqueous Humor Thermo-Fluid-Dynamics in the Anterior Chamber". In: *Recent Developments in Mathematical, Statistical and Computational Sciences*. Ed. by D. Marc Kilgour et al. Cham: Springer International Publishing, 2021, pp. 489-499. ISBN: 978-3-030-63591-6.

[Adl53] Francis Heed Adler. *Physiology of the eye, clinical applications*. St. Louis Mosby, 1953.

[Bau+16] Michaël Baudin et al. "OpenTURNS: An Industrial Software for Uncertainty Quantification in Simulation". In: *Handbook of Uncertainty Quantification*. Ed. by Roger Ghanem, David Higdon, and Houman Owhadi. Cham: Springer International Publishing, 2016, pp. 1-38. ISBN: 978-3-319-11259-6. DOI: 10.1007/978-3-319-11259-6_64-1. URL: https://doi.org/10.1007/978-3-319-11259-6_64-1.

[BBS20] Ajay Bhandari, Ankit Bansal, and Niraj Sinha. "Effect of aging on heat transfer, fluid flow and drug transport in anterior human eye: A computational study". In: *Journal of Controlled Release* 328 (2020), pp. 286-303. ISSN: 0168-3659. DOI: https://doi.org/10.1016/j.jconrel.2020.08.044. URL: https://www.sciencedirect.com/science/article/pii/S0168365920304855.

[Buf+12] Buffa, Annalisa et al. "A priori convergence of the Greedy algorithm for the parametrized reduced basis method". In: *ESAIM: M2AN* 46.3 (2012), pp. 595-603. DOI: 10.1051/m2an/2011056. URL: https://doi.org/10.1051/m2an/2011056.

[CAS22] Open CASCADE. *SALOME: The Open Source Integration Platform for Numerical Simulation*. 2022. URL: https://salome-project.org.

[Dod+14] H. Dodig et al. "Stochastic sensitivity of the electromagnetic distributions inside a human eye modeled with a 3D hybrid BEM/FEM edge element method". In: *Engineering Analysis with Boundary Elements* 49 (2014). Special issue on Bioelectromagnetics, pp. 48-62. ISSN: 0955-7997. DOI: https://doi.org/10.1016/j.enganabound.2014.04.005. URL: https://www.sciencedirect.com/science/article/pii/S0955799714000708.

[Edg] Engineers Edge. *Convective Heat Transfer Coefficients Table Chart*. https://www.enginearsedge.com/heat_transfer/convective_heat_transfer Coefficients_13378.htm. Accessed: 2022-12-28.

[EG21] Alexandre Ern and Jean-Luc Guermond. *Finite elements II*. Texts in applied mathematics. Cham: Springer International Publishing, 2021.

[EYB89] Nathan Efron, Graeme Young, and Noel A Brennan. "Ocular surface temperature." In: *Current eye research* 8 9 (1989), pp. 901-6.

[Gui+20] Giovanna Guidoboni et al. "Neurodegenerative disorders of the eye and of the brain: A perspective on their fluid-dynamical connections and the potential of mechanism-driven modeling". en. In: *Front. Neurosci.* 14 (Nov. 2020), p. 566428.

[Hir+07] A Hirata et al. "Temperature elevation in the eye of anatomically based human head models for plane-wave exposures". In: *Physics in Medicine and Biology* 52.21 (Oct. 2007), pp. 6389-6399. ISSN: 1361-6560. DOI: 10.1088/0031-9155/52/21/003. URL: http://dx.doi.org/10.1088/0031-9155/52/21/003.
[Kos+13] Philip Kosky et al. "Chapter 12 - Mechanical Engineering". In: *Exploring Engineering (Third Edition)*. Ed. by Philip Kosky et al. Third Edition. Boston: Academic Press, 2013, pp. 259-281.
ISBN: 978-0-12-415891-7. DOI: https://doi.org/10.1016/B978-0-12-415891-7.00012-1.
URL: https://www.sciencedirect.com/science/article/pii/B9780124158917000121.

[Lag82] J J W Lagendijk. "A mathematical model to calculate temperature distributions in human and rabbit eyes during hyperthermic treatment". In: *Physics in Medicine and Biology* 27.11 (Nov. 1982), pp. 1301-1311. DOI: 10.1088/0031-9155/27/11/001. URL: https://doi.org/10.1088/0031-9155/27/11/001.

[Li+10] Eric Li et al. "Modeling and simulation of bioheat transfer in the human eye using the 3D alpha finite element method (αFEM)". In: *International Journal for Numerical Methods in Biomedical Engineering* 26.8 (2010), pp. 955-976. DOI: https://doi.org/10.1002/cnm.1372. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.1372.

[Map68] R. Mapstone. "Measurement of corneal temperature". In: *Experimental Eye Research* 7.2 (1968), 237-IN29. ISSN: 0014-4835. DOI: https://doi.org/10.1016/S0014-4835(68)80073-9. URL: https://www.sciencedirect.com/science/article/pii/S0014483568800739.

[Mas04] Barry R. Masters. "Franz Frankhauser and Sylwia Kwasniewska (eds): Lasers in Ophthalmology—Basic, Diagnostic and Surgical Aspects: A Review: Kugler Publications, The Hague, The Netherlands, 2003, pp. 450 (ISBN 90-6299-189-0)". In: *Graefe's Archive for Clinical and Experimental Ophthalmology* 242.5 (Apr. 2004), pp. 446-447. ISSN: 1435-702X. DOI: 10.1007/s00417-004-0873-3. URL: http://dx.doi.org/10.1007/s00417-004-0873-3.

[Mil93] Zdeněk Milka. "Finite element solution of a stationary heat conduction equation with the radiation boundary condition". eng. In: *Applications of Mathematics* 38.1 (1993), pp. 67-79. URL: http://eudml.org/doc/15737.

[Ng 06] Ng, E.Y.K. and Ooi, E.H. "FEM simulation of the eye structure with bioheat analysis". In: *Computer Methods and Programs in Biomedicine* 82.3 (2006), pp. 268-276. ISSN: 0169-2607. DOI: https://doi.org/10.1016/j.cmpb.2006.04.001. URL: https://www.sciencedirect.com/science/article/pii/S0169260706000708.

[NO07] E.Y.K. Ng and E.H. Ooi. "Ocular surface temperature: A 3D FEM prediction using bioheat equation". In: *Computers in Biology and Medicine* 37.6 (2007), pp. 829-835. ISSN: 0010-4825. DOI: https://doi.org/10.1016/j.compbiomed.2006.08.023. URL: https://www.sciencedirect.com/science/article/pii/S0010482506001557.

[ON08] Ean-Hin Ooi and Eddie Yin-Kwee Ng. "Simulation of aqueous humor hydrodynamics in human eye heat transfer". In: *Computers in Biology and Medicine* 38.2 (2008), pp. 252-262. ISSN: 0010-4825. DOI: https://doi.org/10.1016/j.compbiomed.2007.10.007. URL: https://www.sciencedirect.com/science/article/pii/S001048250700176X.

[ON09] E. H. Ooi and E. Y. K. Ng. "Ocular Temperature Distribution: a Mathematical Perspective". In: *Journal of Mechanics in Medicine and Biology* 09.02 (2009), pp. 199-227. DOI: 10.1142/S0219519409002936.

[Pru+01] C. Prud'homme et al. "Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods". In: *Journal of Fluids Engineering* 124.1 (Nov. 2001), pp. 70-80. ISSN: 0098-2202. DOI: 10.1115/1.1448332. eprint: https://asmedigitalcollection.asme.org/fluidsengineering/article-pdf/124/1/70/5902394/70\_1.pdf. URL: https://doi.org/10.1115/1.1448332.

[Pru+23] Christophe Prud'homme et al. *feelpp/feelpp: Feel++ Release V111 alpha.5*. Version v0.111.0-alpha.5. Aug. 2023. DOI: 10.5281/zenodo.8272196. URL: https://doi.org/10.5281/zenodo.8272196.

[PSS21] Christophe Prud'homme, Lorenzo Sala, and Marcela Szopos. "Uncertainty propagation and sensitivity analysis: results from the Ocular Mathematical Virtual Simulator". In: *Mathematical Biosciences and Engineering* 18.3 (2021), pp. 2010-2032. ISSN: 1551-0018. DOI: 10.3934/mbe.2021105. URL: https://www.aimpress.com/article/doi/10.3934/mbe.2021105.

[PW05] Christine Purslow and James S Wolffsohn. "Ocular surface temperature: a review". en. In: *Eye Contact Lens* 31.3 (May 2005), pp. 117-123.
[QMN16] Alfio Quarteroni, Andrea Manzoni, and Federico Negri. *Reduced Basis Methods for Partial Differential Equations*. Springer International Publishing, 2016. DOI: 10.1007/978-3-319-15431-2. URL: https://doi.org/10.1007/978-3-319-15431-2.

[RF77] Robert F. Rosenbluth and Irving Fatt. "Temperature measurements in the eye". In: *Experimental Eye Research* 25.4 (1977), pp. 325-341. ISSN: 0014-4835. DOI: https://doi.org/10.1016/0014-4835(77)90100-2. URL: https://www.sciencedirect.com/science/article/pii/0014483577901002.

[RHP07] G. Rozza, D. B. P. Huynh, and A. T. Patera. "Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations". In: *Archives of Computational Methods in Engineering* 15.3 (Sept. 2007), pp. 1-47. ISSN: 1886-1784. DOI: 10.1007/BF03024948. URL: https://doi.org/10.1007/BF03024948.

[RHP08] G Rozza, D B P Huynh, and A T Patera. "Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations". en. In: *Arch. Comput. Methods Eng.* 15.3 (Sept. 2008), pp. 229-275.

[Sai+24] Thomas Saigre et al. "A coupled fluid-dynamics-heat transfer model for 3D simulations of the aqueous humor flow in the human eye". In: *8th International Conference on Computational and Mathematical Biomedical Engineering - CMBE2024 Proceedings*. Ed. by P. Nithiarasu and R. Löhner. Arlington (Virginia), United States, June 2024. URL: https://www.compbio.med.net/2024/cmbe-proceedings.htm.

[Sal+] Lorenzo Sala et al. "The ocular mathematical virtual simulator: A validated multiscale model for hemodynamics and biomechanics in the human eye". In: *International Journal for Numerical Methods in Biomedical Engineering* (), e3791. DOI: https://doi.org/10.1002/cnm.3791. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cnm.3791. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.3791.

[Sal16] Lorenzo Sala. *Eye2brain*. http://www.cemosis.fr/projects/eye2brain/. 2016.

[Sco88] J.A. Scott. "A finite element model of heat transport in the human eye". In: *Physics in Medicine and Biology* 33.2 (1988), pp. 227-242. DOI: 10.1088/0031-9155/33/2/003. URL: https://doi.org/10.1088/0031-9155/33/2/003.

[Sob93] Ilya M. Sobol. "Sensitivity Estimates for Nonlinear Mathematical Models". In: 1993.

[Šuš+21] Anna Šušnjara et al. "Stochastic-Deterministic Thermal Dosimetry Below 6 GHz for 5G Mobile Communication Systems". In: *IEEE Transactions on Electromagnetic Compatibility* 63.5 (2021), pp. 1667-1679. DOI: 10.1109/TEMIC.2021.3098431.

[Šuš+22] Anna Šušnjara et al. "Uncertainty quantification and sensitivity analysis of transcranial electric stimulation for 9-subdomain human head model". In: *Engineering Analysis with Boundary Elements* 135 (2022), pp. 1-11. ISSN: 0955-7997. DOI: https://doi.org/10.1016/j.enganabound.2021.10.026. URL: https://www.sciencedirect.com/science/article/pii/S0955799721003192.

[too22] MMG tools. *MMG: a mesh-based finite element library*. 2022. URL: https://www.mmgtools.org/.

[Wan+16] Wenjia Wang et al. "Fluid and structure coupling analysis of the interaction between aqueous humor and iris". In: *BioMedical Engineering OnLine* 15.2 (Dec. 2016), p. 133. ISSN: 1475-925X. DOI: 10.1186/s12938-016-0261-3. URL: https://doi.org/10.1186/s12938-016-0261-3.