The ocular mathematical virtual simulator: A validated multiscale model for hemodynamics and biomechanics in the human eye

We present our continuous efforts from a modeling and numerical viewpoint to develop a powerful and flexible mathematical and computational framework called Ocular Mathematical Virtual Simulator (OMVS). The OMVS aims to solve problems arising in biomechanics and hemodynamics within the human eye. We discuss our contribution towards improving the reliability and reproducibility of computational studies by performing a thorough validation of the numerical predictions against experimental data. The OMVS proved capable of simulating complex multiphysics and multiscale scenarios motivated by the study of glaucoma. Furthermore, its modular design allows the continuous integration of new models and methods as the research moves forward, and supports the utilization of the OMVS as a promising non‐invasive clinical investigation tool for personalized research in ophthalmology.


| INTRODUCTION
The interest in patient-specific mathematical models applied to biomedical problems has greatly increased in the last years.To provide meaningful and complementary insights to traditional research, these models should build upon clinical measures inputs combined with physiological knowledge.In order to fully exploit their quantitative predictive capabilities, a complete pipeline should next be developed and should incorporate: data integration, model derivation, numerical solving, validation, and uncertainty quantification.Significant research progress was accomplished towards these challenging goals in specific biomedical applications, for example, cardiovascular simulations 1,2 or cerebral hemodynamics, 3,4 which are now quite mature and well-established research fields.Moreover, efficient computational frameworks 5,6 oriented towards medical applications successfully tackled such complex multiphysics and multiscale problems and provided useful insights to interesting clinical questions.
In the field of ophthalmology, a similar paradigm is needed: the availability of rich, heterogeneous data published in the literature, possibly making contradictory statements and the lack of understanding of the underlying mechanisms of several ocular diseases, calls for innovative approaches to help diagnosis and monitoring of these clinical conditions.For a recent review of the state-of-the-art and the open questions, see for instance 7,8 and the references therein.However, computational and mathematical ophthalmology is still an emergent field and to the best of our knowledge, only few studies focused on this particular topic, as reviewed for example in References 9,10 or, in the context of uncertainty quantification, in References 11,12.The present work aims to contribute to this growing area of research by proposing a novel modeling and simulation environment, called the Ocular Mathematical Virtual Simulator (OMVS).From a clinical viewpoint, the focus of the present work is on glaucoma.The reasons are manifold: it is a leading cause of irreversible blindness worldwide, its prevalence increases with age, it poses a significant public health burden, it currently lacks a cure, existing treatments focus on managing the condition and slowing its progression.Our collaborative research with the specialist center at the Icahn School of Medicine, Mount Sinai (NY) provides a unique opportunity to develop a multidisciplinary approach to tackle these issues.Specifically, the OMVS has the purpose to provide a window on both ocular biomechanics and hemodynamics from a macroscale viewpoint.Moreover, this innovative framework means to account for the complexity stemming from the availability of data and their heterogeneity due different measurement techniques. 13,14The multiscale and multiphysics characteristics dictated by data acquisition and the application itself have been considered from the beginning in the design and accounted thanks to a modular structure.On one hand, the combined effects of ocular blood flow and different ocular tissues are described by a coupled hemodynamics and biomechanics model.On the other hand, the multiscale aspect, essential to properly account for systemic effects of blood circulation coupled with local effects on tissues of interest, is represented by a coupled partial and ordinary differential equations for fluid flow.More details about these multiphysics and multiscale features and the associated developments are provided in Section 2, while Section 4 discusses information about geometrical and functional parameters utilized as OMVS input data.
The next challenge within the OMVS development is to have a robust and efficient numerical strategy able to preserve the physical properties at the discrete level.To address this ambitious goal, suitable innovative methods are adopted: a Hybridazable Discontinuous Galerkin (HDG) method with an original integral boundary condition 15 detailed in Section 3.2 and a novel effective multiscale coupling solving method described in Section 3.1.
A preliminary version of the OMVS has been presented in Reference 16 and subsequently used in Reference 17 to target specific clinical applications related to glaucoma.A detailed uncertainty quantification and global sensitivity analysis study was conducted for a reduced version of the OMVS framework in Reference 18.However, until now, the full OMVS framework and the specific question of its validation have not been tackled directly yet.We therefore describe in Section 5 the results of a thorough validation analysis.Finally, a comprehensive scheme of the different connections between data and models within the OMVS framework can be found in Figure 1, illustrating our paradigm and validation strategy.F I G U R E 1 General design and strategy of the ocular mathematical virtual simulator (OMVS) framework.A MULTISCALE MODEL The OMVS framework is based on a multiscale mathematical description of blood circulation (or hemodynamics) and tissue deformation (or biomechanics) in various parts of the eye.Specifically, the anatomical regions considered in this work are sclera, cornea, lamina cribrosa, retina and choroid, see Figure 2. The multiscale approach is adopted to couple models of different complexity for the hemodynamic and/or biomechanic behaviors in different anatomical regions depending on their relevance for the study at hand, as detailed below.
The current development of the OMVS is motivated by the study of glaucoma, for which the lamina cribrosa is considered to be a site of particular interest due to its crucial role in the balance of pressures between the eye and the brain. 20However, direct measurements within this tissue are difficult or even impossible to obtain in clinical practice, thus alternative approaches are investigated.
2][23][24][25][26] On the other hand, few studies have been completed on the hemodynamics of the lamina. 27For this reason, the lamina cribrosa is modeled here as a three-dimensional (3D) poroelastic medium, with the goal of capturing the interaction between the blood flow through the tissue and the deformation of the tissue itself (see Section 2.2).However, blood flow and tissue deformation occurring within the lamina are also influenced by non-local factors, such as systemic blood pressure, that need to be taken into account in order to obtain physiologically-reasonable solutions (see Section 2.1).
Tissue deformations in the lamina cribrosa are mainly due to the load induced by intraocular pressure (IOP) and cerebrospinal fluid pressure (CSFp) acting on its inner and outer surfaces, respectively.However, these surfaces also experience a transversal tension due to the overall IOP-related inflation of the eye ball.This transversal tension can be incorporated in the model in several ways: either as given boundary conditions similarly to the approach proposed in Reference 28, or as interface conditions between the lamina cribrosa and the sclera.In the present work, we employ the latter strategy, which requires the lamina cribrosa to be coupled with models for other ocular tissues via interface conditions expressing the continuity of displacements and stresses (see Section 2.4).The other ocular tissues we account for in the present contribution are sclera, cornea, retinal nerve fiber layer and choroid, described as elastic media (see Section 2.3) This approach has the advantage of providing a global description of the biomechanics of main tissues composing the eye ball.
The blood supply to the lamina cribrosa is provided by the posterior ciliary arteries, which branch out from the ophthalmic artery before the central retina artery (CRA), whereas the blood drainage occurs via the central retinal vein (CRV), which also drains the retina (see Figure 3).Thus, the hemodynamic conditions within the lamina cribrosa are inherently coupled with the hemodynamics in the central retinal vessels and in the retinal vasculature.To capture this interaction, the circulation in the CRA, CRV and retina are modeled via a zero-dimensional (0D) model based on a nonlinear system of ordinary differential equations (ODEs) describing blood flow as the analogous of an electric current flowing through an electric circuit (see Section 2.1).Next, we employ a multiscale connection between the 0D network F I G U R E 2 Anatomy of the eye.Image created with BioRender. 19nd the 3D description of the hemodynamics in the lamina cribrosa.In particular, the idea is to have a first approximation of the contribution of the lamina cribrosa perfusion incorporated in the retinal vasculature via the 0D model described in Section 2.1.Then a more detailed 3D lamina cribrosa representation is introduced by means of the coupling between the hemodynamics and the biomechanics via the poroelastic model presented in Section 2.2.A suitable exchange of information between the 0D model for the circulation in CRA/CRV/retina and the 3D poroelastic medium for the circulation in the lamina cribrosa must be prescribed in order to ensure the continuity of mass and pressure (see Section 2.4).Finally, the description of the geometry, computational mesh, and parameter values utilized in the models is postponed to Section 4.

| A 0D network model for the circulation in the retina and the central retinal vessels
A nonlinear system of ODEs is utilized as a simplified representation of the blood circulation in the retina, CRA and CRV.The blood circulation is modeled via a lumped-parameter model, which exploits the electric analogy to fluid flow 29 Chapter 15.Within this analogy, the flow of a fluid through a hydraulic network corresponds to the flow of an electric current through an electric circuit.Thus, volume, flow rate, fluid velocity, pressure correspond to electric charge, current, current density, potential, respectively.
The ODEs model illustrated in Figure 4 stems from the circuit described and validated in References 30 and 31.The vasculature is divided into six main compartments: lamina cribrosa (lc), CRA (cra), arterioles (r,a), capillaries (r,c), venules (r,v), and CRV (crv ð Þ are exposed to the IOP and the retrobulbar segments R cra,1 , R cra,2 , R crv,3 , R crv,4 ð Þ are exposed to the pressure within the optic nerve tissue, also referred to as retrolaminar tissue pressure (RLTp).IOP is mainly due to the pressure of aqueous humor in the anterior chamber, whereas RLTp is mainly due to the pressure of the cerebrospinal fluid in the subarachnoid space. 32The lamina cribrosa compartment, with inner pressure P lc , corresponds to a 0D simplified description of the hemodynamics in this tissue.
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 the simplified FSI (fluid-structure interaction) model studied in Reference 33 on the CRA/CRV interaction with the lamina cribrosa.In the present lumped-parameter model (Figure 4), this contribution has been simplified and the external pressure felt by the CRA/CRV is only due to the IOP in the intraocular and translaminar segments, and to the RLTp in the retrobulbar segments.
Similarly to the work presented in References 29,31,34, we employ the following modeling choices for the nonlinear resistors in the circuit: where L is the length of the vessel, b p is the average pressure inside the vessel, p e is the external pressure to the vessel, and k 0 , k L and K p are univocally characterizing the nonlinear resistive behavior.Equation (1) assumes that the cross-section remains circular in presence of external pressures; this situation applies to R cra,i , i ¼ 1, …, 4. Equation (2) models a Starling resistor behavior, which is based on experiments suggested by References 35,36.In particular, this formula expresses the fact that, in case of compressible tubes, the circular cross-section is preserved for dilation, whereas when the transmural pressure becomes negative the cross-section shape reflects the physiological high collapsibility of the venous segments.R r,v1 , R r,v2 , R crv,i , i ¼ 1, …,4 follow this equation.

| A 3D poroelastic model for the lamina cribrosa
In this section we introduce the 3D model for the lamina cribrosa.We describe it as a poroelastic medium, where hemodynamics and biomechanics are combined with the goal of studying their interactions, as they are essential for the physiological functions of this tissue.To this end, let us consider the poroelastic system introduced by Biot 37 and subsequently studied by References 38-41: where Ω & ℝ d , d ¼ 2, 3 is the computational domain, u is the solid displacement, σ is the stress tensor of the mixture, also known as total stress, F el is the volumetric force term, p is the fluid pressure, α is the Biot coefficient, j is the discharge velocity, F fl is the volumetric fluid source term, and ζ is the fluid content.The constitutive laws for ζ, σ, and j are: where M is the Biot modulus, μ and λ are the Lamé parameters of the elastic matrix, and K is the permeability tensor.
In the context of this study, Equations (3a) and (3b) represent the lamina biomechanics and hemodynamics, respectively.Guided by the particular application at hand, we make the following assumptions and considerations: • The volumetric source term in Equation (3b) is assumed to be zero, that is, F el ¼ 0. This implies that stresses and strains inside the lamina are merely due to its boundary conditions, that is, IOP and CSFp, and interface conditions, that is, coupling with other ocular tissues (see Section 2.4).• The volumetric source term in Equation (3a) is assumed to be zero, that is, F fl ¼ 0. This assumption implies that blood circulating in the lamina is solely provided by a multiscale effective connection between the lamina and the posterior ciliary arteries on the arterial side, and between the lamina and the CRV on the venous side (see Section 2.4).• Changes in fluid content are assumed to be mostly due to changes in pressure, so that the term α ∂ ∂t r Á u ð Þ can be neglected in Equation (4a).Thus, we obtain the following simplified system of PDEs: Referring to the same constitutive laws described in (4), system (5a) and (5b) is accompanied by appropriate boundary (BCs) and initial conditions (ICs), which are discussed in Section 2.4.
Remark 1.The well posedness of the poroelastic system (3) is proved in References 42,43.
Remark 2. In the context of lamina cribrosa modeling, various extensions of the Biot system have been considered, such as the use of a nonlinear model to describe its biomechanics 44 and the use of a poroviscoelastic model to characterize the interaction between the blood flow and the elastic tissue within the optic nerve head. 43,45,46In this work, we opted to consider the simplified system (5) since the focus of our study is the coupling between the lamina biomechanics and hemodynamics and the circulation upstream and downstream of the lamina (Section 2.1).Thus, the simplified system (5) should be seen as a first step to create the "big picture" of an ocular system in which each mathematical component could be replaced with more sophisticated models adapted as the research evolves.
To simplify the presentation, we detail hereafter the hemodynamics and the biomechanics of the model proposed in Equation ( 5) separately as sketched in Figure 5.

| Hemodynamics of the lamina cribrosa
We consider here just the contribution of Equation (5b) and its mixed formulation: find j, p such that endowed with appropriate initial conditions, where Ω lc corresponds to the lamina cribrosa domain.The choice of the boundary conditions for the fluid flow through the lamina is of particular importance, given that blood enters the lamina from the branches of the posterior ciliary arteries located at its outer boundary and it is drained by the central retinal vein running through the lamina transversally in its course from the retina to the cavernous sinus.These concepts are schematized in Figure 6A.We adopt the isotropic hypothesis (Figure 6B) for the permeability tensor K which accounts for capillaries within the tissue of the lamina.Thus, we assume that K ¼ κ I with a given permeability κ > 0. This approximation hypothesis has been also supported by experimental studies. 47

| Biomechanics of the lamina cribrosa
For what concerns the biomechanics, we consider the contribution of Equation (5a) and its mixed formulation: find σ el,lc ,u lc such that HEMODYNAMICS BIOMECHANICS  together with appropriate initial conditions.Specifically, in the biomechanics part, we consider the contribution of the extracellular matrix, the collagen beams and nerve bundles as described in Figure 6B, but we neglect the biomechanical behavior of the blood vessels.

| A 3D biomechanical model for the sclera, choroid, retina, and cornea
In this section we discuss the model that represents the biomechanics within the sclera (Figure 7A), choroid (Figure 7B), retina (Figure 7C) and cornea (Figure 7D).In this case we adopt a mixed formulation of the classic linear elastic problem: find σ el,i , u i defined on endowed with appropriate initial conditions.
Ocular tissue domains introduced in the biomechanical description of the ocular mathematical virtual simulator.

| Boundary and interface conditions in the OMVS framework
This section illustrates the boundary conditions that are adopted for the PDE systems and the interface conditions that describe the coupling among all the model blocks within the multiscale system.

| Hemodynamics within the lamina cribrosa
The conditions to be imposed at the surfaces delimiting the lamina cribrosa are depicted in Figure 8. Specifically, we have: where is the boundary delimiting the lamina cribrosa and n is the outward unit normal vector.
The notations Q I t ð Þ and P I t ð Þ refer to the variables obtained from the 0D model (Section 2.1).These conditions replicate the physiology of the lamina cribrosa.This tissue is nourished through the lateral boundaries, and finally the blood is gathered by the CRV at the central opening (Figure 6A).Equation (9a) describes the physiological drainage through the CRV in a simplified manner, via a Dirichlet boundary condition, where the value of p int,hemo is given (see Section 4).Equation (9b) corresponds to the physiological situation where there are no outward fluxes from the top or the bottom surfaces of the lamina.Equations (9c) represents the nourishment of the lamina cribrosa.This condition describes the feedback from the ODEs circuit to the hemodynamics part of the poroelastic model.The main advantages are: •this framework provides a zoom in a region of interest for the lamina cribrosa, while accounting for systemic feedback via the 0D circuit; •the value P I is constant on Γ ext thanks to the integral condition we have adopted (Equation (9c)): •the mixed formulation in Equation ( 6) provides the natural spatial multiscale connection between the 3D flux R Γ ext j Á n and its 0D representation in the circuit Q I ð Þ.
F I G U R E 8 Hemodynamics boundary conditions for the lamina cribrosa in the ocular mathematical virtual simulator.

| Biomechanics of the lamina cribrosa
The interface and boundary conditions pertaining to the lamina biomechanics are summarized in Figure 9.
Interface conditions are imposed on the external lateral boundary Γ ext .Here the lamina connects with the sclera and we require the continuity of displacement and stress as follows: The BCs are depicted in Figure 9, that is, we have: • on the internal boundary of the lamina Γ int in this case due to lack of biological information we assumed that all the internal boundaries feel the same pressure (more details on the value of this pressure in Section 4) that is compressing (negative sign) 8t 0, T ½ ; • on the top boundary Γ top the negative sign is to indicate the fact that IOP is compressing the tissue; • on the bottom boundary Γ bottom F I G U R E 9 Biomechanical boundary and interface conditions for the lamina cribrosa in the ocular mathematical virtual simulator.
also in this case the negative sign is to indicate the compressing effect of the CSF behind the eye on the tissue.

| Biomechanics of sclera, choroid, retina and cornea
First we list the boundary conditions associated with the description given in system (8) for the biomechanical behavior of sclera, choroid, retina and cornea, where.10 and 11).
Thus, we have for the sclera: For the choroid, BCs read: Boundary condition labels of the OMVS geometry.The original geometry has been clipped along the z-plane in order to view within the eye.For each subfigure, the domain of interest is highlighted in red.
For the retina we have: Finally, we impose the following BCs for the cornea: Second, we impose the continuity of stresses and displacements on the interfaces between these ocular tissues, namely sclera, choroid, retina and cornea:

| NUMERICAL METHODS AND COMPUTATIONAL FRAMEWORK
In this section we present the discretization methods adopted in order to numerically solve the multiscale model introduced in Section 2. First, we describe the time discretization approach, in particular an innovative effective treatment Interface condition labels of the OMVS geometry.The original geometry has been clipped along the z-plane in order to view within the eye.For each subfigure, the domain of interest is highlighted in red.
of the multiscale description of the hemodynamics in the OMVS framework (see Section 3.1).Next, we introduce the spatial discretization approach, in particular the Hybridizable Discontinuous Galerkin (HDG) method with an original non-standard Integral Boundary Condition (IBC-see Section 3.2).Finally, we briefly present the open-source computational framework in which the numerical methods are implemented.

| Time discretization
First, for what concerns the time approximation of the biomechanics problem, namely systems ( 7) and ( 8), we adopt a standard BDF2 approach, see for instance 48 Chapter 5.
Next, we focus on the hemodynamics problem and provide details about an innovative method we adopted to address the time discretization of this multiscale system.The main challenge is to numerically solve in an efficient manner the system describing the 3D lamina cribrosa hemodynamics in the domain Ω lc and the ocular posterior vasculature hemodynamics described by the network model represented in Figure 4 and denoted hereafter Υ.Note that the problem involves several nonlinearities, for instance the interaction between IOP/RLTp and the blood pressure within the CRA and CRV (Equations ( 1) and ( 2)).
On the one hand, recall that the system of PDEs (6) in Ω lc is endowed with the boundary conditions specified in Equations (9a), (9b) and (9c), and accounts for the condition expressed in Equation (10).On the other hand, the circuit Υ can be uniquely described in a generic manner by the following nonlinear system of ODEs: with the initial conditions y t ¼ 0 ð Þ¼y 0 .If d is the number of unknowns of the circuit, that is, the number of ODEs to be solved, y is a d-dimensional vector representing the state variables of the circuit Υ, A is a d Â d tensor including the topology and physics of the connections among the circuit nodes, and r y,t is composed by sources and sinks within the circuit, and by the contribution due to the coupling with the PDE region Ω lc .In addition, the conservation of flux between the 3D and the 0D is given by where R Γ ext j Á n is the outward flux from the domain Ω lc , P I is the pressure at the interface Γ ext and P in is the unknown pressure at the entrance of the lamina cribrosa in the network Υ (see Figure 4).As such, P in is part of the unknown vector y.
Let t n ¼ nΔt with Δt > 0 the time-discretization step.We introduce the notation for the pressure in the 3D domain Ω lc and the vector of unknowns in the circuit Υ: In particular, we rewrite y as to highlight the crucial role played by the unknown P I t ð Þ in the subsequent algorithm: Note that Steps 1 and 2 are defined on the discrete time interval t n ,t nþ1 ð Þ, but the differential operators have yet to be fully discretized in time and space.
In this work, we adopted a partitioned approach for the numerical solving of the multiscale problem, in order to preserve the modular structure of the OMVS framework and in the perspective of potential extensions that would imply changes limited to each building block.Thus, the communication between the two steps is only realized through the initial conditions and the proposed effective strategy allows notable flexibility in the choice of solution methods for each sub-problem resulting from this specific algorithm.However, the design of the two steps and in particular solving the PDE system jointly with the condition of the conservation of the fluxes (Equation ( 13)) in Step 2, is meant to preserve an implicit treatment of the multiscale connection.
Remark 3. Other partitioned strategies could be alternatively employed to numerically solve the multiscale system under consideration in this contribution.One possibility, referred to as functional iterations, or fixed-point (Picard) iteration is to design a strategy based on sub-iterations between PDE and ODE solvers at each time step.Another approach addressing the need to solve in separate sub-steps the PDE and the ODE problem in based on operator splitting.An extensive discussion about the attractive features of these strategies and their limitations in the context of the coupling between a poroelastic medium and a lumped hydraulic circuit can be found in Reference 49.The implementation of such approaches and a thorough analysis of their potential advantages in the present context raises interesting questions and might be considered as future research directions.
The computational framework relies on the open-source finite element library Feel++, 50 in which the global discretization procedure was implemented.Feel++ is a Finite Element Embedded Library in C++ that allows using a very wide range of Galerkin methods, as well as other advanced numerical methods such as domain decomposition methods or certified reduced basis.The ingredients of the software include a very expressive embedded language, seamless interpolation, mesh adaption and seamless parallelization.The first step was implemented within OpenModelica, 51 an open-source Modelica-based modeling and simulation environment intended for industrial and academic studies of complex dynamic systems.The multiscale connection between the different components of the solver is handled thanks to the capabilities of the Feel++ library.As for the spatial discretization, the splitting approach provided by Algorithm I is very flexible, allowing for different methods.Our specific choice in the present work is the Hybridizable Discontinuous Galerkin method, as described in the next section.

| Spatial discretization
In this section we briefly present the Hybridizable Discontinuous Galerkin (HDG) method applied for the spatial discretization of the model, which was originally introduced in this context in Reference 15.The rationale for adopting this discretization approach is the high accuracy that all the model solvers require to reproduce physiological behaviors.Specifically, the advantages of HDG are that: i. the equation are enforced element-by-element leading to local conservation properties; ii. the optimal convergence properties for both primal (potential and displacements) and dual (stresses and fluxes) variables are ensured; in particular fluxes are crucial in the communication between the micro-and macroscale.
However, these features induce a significant computational cost, which can be largely mitigated using static condensation as described in Reference 15 Section 4.2.
Remark 4. The use of HDG formulation combined with IBC R Γ lc,I j Á n ¼ Q I supports the direct solution of the pressure by means of a natural coupling between the 3D and the 0D parts of the model, without any sub-iteration.
The main technical details of this approach are provided in Appendix A, where we report the formulation for the elasticity and the porous media model, which were not detailed in the original paper. 15The final forms of the hemodynamics and biomechanics systems, obtained from the static condensation algorithm, are expressed in terms of the potential or the displacement, respectively.Potential and flux-or displacement and stresses for elasticity-are then computed element wise at higher order in a post-processing step.For the present application, this property is fundamental, since it provides the same level of accuracy for the primal variables (potential and displacements) and for the dual unknowns, fluxes and stresses of notable interest here.
The overall discretization procedure is again implemented in the Feel++ library. 50Numerically the problem is solved by means of domain decomposition methods, in parallel, using an algebraic multigrid solver for both Darcy and elasticity systems; more details are available in Reference 15.

| OMVS INPUT DATA
The purpose of this section is twofold: (i) first, we describe the pipeline we developed to generate a realistic computational mesh of the entire human eye; (ii) second, we provide a thorough description of the input parameters of the OMVS and the main rationale behind our choices.We emphasize the important role of the information described in this section in view of the reproducibility of our work.

| Geometry and mesh
In this subsection we present the process we developed to generate a realistic geometry and the associated computational mesh used for the 3D biomechanics and hemodynamics simulations of the tissues represented in the OMVS.This geometry is parameterized on patient-specific inputs, such as the distance between the cornea and the lamina cribrosa and some lamina cribrosa features as described in Table 1.Patient-specific geometrical inputs could be relevant in several clinical contexts, such as: (i) higher risk factors for glaucoma among people of African descent (AD) when compared to European descent (ED), due to differences in their ocular structures (thinner corneas 52 and thinner sclera and lamina cribrosa 53 for the former compared to the latter); (ii) dynamic changes taking place in the anterior segment with aging, suggesting prevalence of angle closure glaucoma increased in females after middle age 54 or (iii) structural changes in patients affected by myopia. 55he generation of the geometry and the mesh is realized via the software Salome, 56 which is an integration platform for numerical simulations with pre-and post-processing tools.We specifically utilized the CAD (Computer Aided Design) module, which is based on OpenCascade 57 and its constructive solid geometry features.In the following we briefly describe the main steps of the pipeline we designed to obtain a computational mesh of the eye starting from a CAD image (as schematically illustrated in Figure 12).

| Geometry
The starting CAD image of the eye is illustrated in Figure 12A (courtesy of CEMOSIS 1 ).This preliminary image is imported in Salome, where the following 3D ocular structures are identified: cornea, vitreous humor, iris and ciliary body, ligament, lens, sclera, choroid, retina, central retinal vein and artery.Furthermore, we added the lamina cribrosa representation, motivated by the special role that this tissue plays in the onset and progression of glaucoma. 20On the basis of the ocular anatomy description extracted from, 58 we created a lamina cribrosa geometry within the imported ocular framework.In the generation of the geometry, we made some specific parameters accessible and changeable by the user, as described in Table 1.A cut of the newly created geometry of the eye is displayed in Figure 12B.

| Mesh
The mesh module in Salome is used to generate 1D/2D/3D meshes.In particular we employed the NETGEN meshing plug-in 59 to specify the parameters of the mesh and eventually build the mesh.We set markers on every volume or surface of the structured mesh allowing simulations that can be run only on single parts of the eye or on the whole geometry.Figure 13A displays the cut along the z-plane of the overall computational mesh, while Figure 13B shows the generated computational mesh for the lamina cribrosa.In summary, from a CAD of the whole eye we developed a parameterized ocular geometry in order to generate a complex computational mesh, with the following characteristics: • 20 boundaries and interfaces, with conforming elements at edges and at interfaces; • specific discretization for the lamina cribrosa due to its different geometrical dimensions with respect to the other components.
This generating process has addressed the goal of having a full realistic description of the ocular architecture, which is used by the OMVS for patient-specific numerical simulations.For the entire process and in the case of default values, the code takes a total time of 272:6s, of which 186:1 for the geometry and 84:7 for the mesh.

| Model parameters
In this subsection we report all the baseline values for the parameters employed in the OMVS.The selection was made in such way that, at baseline, the model predictions agree with values reported in the literature.Some of the parameters were directly measured or extracted from literature, whereas others were not readily available and thus are indirectly inferred from accessible data, under clinical hypothesis, which will be indicated when relevant.In addition, whenever pertinent, a brief discussion on the measurement instruments to obtain or estimate these parameters is presented.
We remark that parameters choices and calibration is well-recognized as a very delicate problem in the context of mathematical modeling in medicine.On the one hand, as explained above, suitable data to feed the models might be difficult or even impossible to be measured directly.On the other hand, they usually suffer from high variability due to (i) operator-dependent and instrument-dependent measurements, which can be different across clinical centers and (ii) inherent differences among individuals.To take into account these sources of uncertainty and variability, we refer to the preliminary sensitivity analysis we carried out in Reference 18 for the 0D components of the OMVS.
We divide the set of parameters of the proposed ocular mathematical model in three groups: (i) parameters common to both the ODE-based and the PDE-based systems, called hereafter global OMVS parameters, (ii) specific parameters employed only in the 3D hemodynamical and biomechanical description of the model, and (iii) specific parameters related to the 0D hemodynamical part.
(A) Cut along the z-plane.
(B) Zoom on the lamina cribrosa.
F I G U R E 1 3 Mesh of the eye realized with Salome using NETGEN algorithms.

| Global OMVS parameters
First, we focus on the systemic parameters, namely the heart rate (HR), the systolic blood pressure (SP), and the diastolic blood pressures (DP), that can be easily measured with standard non-invasive devices.Specifically, blood pressure is measured with a sphygmomanometer, a clinical instrument which typically consists of an inflatable rubber cuff that is applied to the arm and connected to a column of mercury next to a graduated scale, enabling the determination of SP and DP by increasing and gradually the pressure in the cuff.With the help of a stethoscope, the reading of the HR is easy and noninvasive.Second, the parameters pertaining to the ocular region described in the model are the IOP, and the RLTp.Recall that IOP is the pressure of the fluids inside the eye and it is determined by the balance between the production and the drainage of aqueous humor.Its value can be obtained non-invasively in a standard clinical setting with the Goldmann Application Tonometer, which measures the force necessary to flatten an area of the cornea and computes the pressure exploiting the Imbert-Fick law. 60The RLTp is the pressure within the optic nerve tissue, behind the lamina cribrosa.This pressure is largely determined by the pressure within the subarachnoid space, that is, the intracranial pressure (ICP). 32ICP is usually measured with a lumbar puncture, an invasive technique in which a needle is inserted into the spinal canal to collect cerebrospinal fluid (CSF).However, some recent clinical works 61,62 and mathematical models, 63 propose innovative methods in order to estimate ICP noninvasively.The baseline values for the global parameters are reported in Table 2.

| 3D hemodynamics and biomechanics parameters for the OMVS
Here we describe the baseline parameters values related to the 3D modeling of the OMVS.
The baseline permeability coefficient of the lamina cribrosa and the baseline Lamé coefficients for the different ocular tissues are reported in Table 3.In particular we recall that for the lamina cribrosa we have adopted the isotropic hypothesis (Figure 6B) for the permeability tensor K ¼ κI which takes into account for capillaries within the tissue of the lamina.This hypothesis provides a satisfactory approximation of physiological conditions as confirmed by experimental studies. 47The majority of the Lamé coefficients values have to be considered valid in case of small stresses 0:0 À 8:0kPa ≈ 0 À 60 mmHg ð Þ , see Reference 68, which includes the physiological range of pressures involved in the realistic framework developed in the present work.
To complete the set of parameters involved in the 3D model, we consider p int,hemo ¼ 19 mmHg as baseline value in Equation (9a) and p int,mech ¼ 40:88 mmHg as baseline value in Equation (11c).These values have been retrieved from Reference 30 in the case of baseline systemic blood pressure of SP ¼ 120 mmHg and DP ¼ 80 mmHg, as in Table 2.For the different validation cases for which we report simulations in the next section we have calibrated the values of p int,hemo and p int,mech with a similar approach, namely directly proportional to the values of SP and DP.Thus, every time the values of SP and DP are modified from their baseline values, also p int,hemo and p int,mech are changed accordingly.
Regarding available measurements in this context, the utilization of Optical Coherence Tomography Angiography offers valuable insights about crucial dimensions, thickness, and vessel structure of the lamina and retina. 74,75Additionally, the pachymeter serves as a precise tool for obtaining Central Corneal Thickness measurements.In terms of material properties, our focus turns to ultrasound elastography, a technique that holds promise in this domain.Moreover, we point to a rich reservoir of ocular biomechanics literature, particularly highlighting ex-vivo measurements and traction tests supported by advanced mathematical analysis. 76

| 0D hemodynamics parameters for the OMVS
In this paragraph we define the baseline parameter values for the 0D component of the OMVS.The resistances and the capacitances baseline values are reported in Tables 4 and 5, respectively, together with the corresponding references from which they were extracted.In the particular case of nonlinear resistors, the values are computed using Equations ( 1) and ( 2), respectively.Finally, we characterize the blood pressure source and sink within the circuit in a consistent manner with experiments.For P eye,in , we employed the same approach as in Reference 30: we reconstruct the pressure profile at the entry of the circuit (see Figure 4) from typical Color Doppler Imaging (CDI) measurement of the blood velocity in the CRA and we impose it as a time-dependent pressure source (see Figure 14).The CDI employs the Doppler effect to generate imaging of the movement of tissues and body fluids-usually blood-and their relative velocity to the probe by calculating the frequency shift of a particular sample volume.We parametrize the reconstructed signal using HR, SP and DP as inputs; thus, we have divided this blood pressure time profile into six parts for each cardiac cycle, as follows:  where b t ¼ mod t, 60 HR À Á .Regarding P Eye,out , we set its baseline value to 14 mmHg, as reported in Reference 30.For the retinal blood flow another quantity that the ophthalmologist can measure is the oxygenation of the blood by means of the retinal oxymetry, 79 however for the current implementation of the OMVS, this measurement is not exploited.

| OMVS VALIDATION
In this section we discuss several significant outcomes of the OMVS and compare our simulations with clinical or experimental data whenever available.The aim is to provide a thorough validation for the OMVS framework and to complement predictive results previously published in References 17,80, The verification of the numerical methods presented in Section 3 has been already extensively performed in Reference 15.
Figure 15 displays all the output available within the OMVS and lists the main sources where we retrieved the experimental data.The section is organized as follows: • validation of the hemodynamics simulated by the OMVS; • validation of the biomechanics of the lamina cribrosa simulated by the OMVS; • validation of the biomechanics of other ocular tissues represented in the OMVS.
Note that the intrinsic physiological variability makes validation a very challenging task; thus we strive to report both qualitative and quantitative comparisons with experimental data whenever the latter ones are possible.Furthermore, we highlight the fact that for the hemodynamics within the lamina cribrosa, to the best of our knowledge, there are no experimental or clinical data, thus we cannot directly validate the simulation results; however, since all the system components are strongly coupled we consider that if all the other parts reproduce a physiological behavior, also the lamina cribrosa perfusion computed by the OMVS should attain realistic values.

| Hemodynamics of the ocular posterior segment
In this part we discuss the outcome concerning the central retinal vessels hemodynamics, which are part of the circuit presented in Section 2.1 and in the lamina cribrosa, as described in Section 2.2.First, we obtained similar numerical results as reported in Reference 30.Despite the numerous simplifying assumption in the simple 0D scheme, indeed, the proposed mathematical model is able to capture the mechanical action of IOP on clinically measurable hemodynamics quantities, in particular the CRA blood velocity.In particular Guidoboni et al. report a peak systolic CRA flow of 119.4 μ=L and an end distolic CRA flow of 33.7 μ=L for a patient with normal IOP and high blood pressure, values which are in good agreement with our results for the HBP-Normal IOP case (110.58 and 28.45 μ=L, respectively).Therefore their conclusions can be extended also to the present model.The model predicts that the steep decay in retinal blood flow that occurs due to high IOP would shift towards higher values as the blood pressure of the subject increases.This outcome also agrees with the clinical study of He et al., 81 who found that a higher IOP was needed to attenuate ocular blood flow in Long-Evans rats with higher systemic blood pressure.Second, in another set of simulations, we utilize the values of IOP, HR, SP, and DP reported in Reference 82 as input data of the OMVS and compare the predictions of the simulation (Figure 16) with the clinical results provided in Reference 82.We have simulated three test cases of clinical interest: i. systemic hypertensive subjects (HBP) with high IOP; ii.system hypertensive subjects (HBP) with normal IOP; iii.normal systemic blood pressure subjects (NBP) with high IOP.
Input values are listed in Table 6. Figure 16A shows the blood flow profile during two cardiac cycles computed by the OMVS within the CRA and CRV.Note that the CRV flow has a negative sign to indicate its opposite direction with respect CRA flow, the same convention used in clinical settings for CDI (see Figure 14).Systemic hypertensive virtual subjects share a similar blood flow waveform within the CRA (red and blue solid lines), whereas the model predicts a lower CRA blood flow for NBP-High IOP subjects (black solid line).For what concerns the CRV, the simulation results for virtual subjects with elevated IOP (blue and black dashed lines) show a drop in the blood flow with respect to HBP-Normal IOP (red dashed line).This reduction lasts longer in the case of NBP-High IOP (black dashed line).
The top panel of Figure 16B displays the lamina cribrosa pressure drop between its external and internal boundaries during one cardiac cycle.The bottom panel of Figure 16B reports the lamina cribrosa perfusion value on Γ ext during one cardiac cycle.Figure 16B suggests that HBP-Normal IOP subjects (red line) experience a similar hemodynamic behavior as HBP-High IOP subjects (blue line), whereas NBP-High IOP subjects (black line) exhibit a significant drop in the lamina cribrosa pressure and perfusion (up to 28% in both cases).These outcomes are easily identifiable also qualitatively from Figure 16C where we display the 3D spatial distribution of the blood pressure accompanied by the intensity of the blood perfusion (color-coded arrows within the domain).The clinical results reported by Costa et al. 82 suggest that only individuals with high IOP and normal blood pressure may be at higher risk for glaucomatous damage.The model predictions of the OMVS conjecture that elevated IOP has a significant impact on the lamina cribrosa hemodynamics when combined with normal systemic blood pressure, whereas it is remarkably less noticeable in case of high systemic blood pressure.Considering that a deficit in the lamina cribrosa perfusion as a risk factor for glaucoma, which is an increasingly supported idea in ophthalmology, 83,84 the simulation results obtained employing our mathematical model are consistent with the clinical analysis presented above.In addition, we also compared our results with experimental data from a quantitative viewpoint.In particular, we refer to the values of the CRA blood flow at peak systole, end diastole and mean reported by Harris et al., 85 Riva et al., 86 and Dorner et al., 87 which are in good agreement with the present simulation outcomes (see Table 7).
Third, these theoretical predictions on the relationship between intraocular pressure, blood pressure, ocular perfusion and glaucoma have been confirmed by the Singapore Epidemiology of Eye Diseases study, an independent population-based study including nearly 10000 individuals. 88

| Biomechanics of the lamina cribrosa
We focus in this part on the OMVS outcomes related to the biomechanics of the lamina cribrosa.
First, Guidoboni et al. 33 compared the simulation results with experimental measurements 89,90 suggesting that the biomechanics in the lamina has a significant impact also on the hemodynamics in the central retinal vessels.In a previous study 91 (results not reported in detail here) we performed a similar virtual experiment that investigates the effect of IOP on the biomechanics and hemodynamics within the lamina cribrosa.The results obtained employing the OMVS were consistent with the findings reported in Reference 33.
Second, Causin et al. 65 reported results on two different modeling case studies and comparisons with experiments conducted by Yan et al., 92 who mounted three enucleated human eyes on a specially designed experimental apparatus, which allowed to sequentially increase the IOP.We replicated this experiment virtually using the OMVS, completing five tests with increasing IOP values of 5, 15, 25, 35 and 50 mmHg, respectively.According to data reported in Reference 65 and 92, for all virtual tests we considered the following inputs SP ¼ 126:3 mmHg DP ¼ 84:2 mmHg HR ¼ 60 beats= min: Figure 17A shows the deformed 3D geometry of the lamina cribrosa for the five different input values for the IOP, using the one corresponding to IOP ¼ 5mmHg as the reference value (in gray).On this basis, we next computed the LC volume below reference, namely ΔV , and we compared our five simulations with the measurements using three enucleated human eyes (eyeA, eyeB, and eyeC) 92 and the results provided by Reference 65 in the two modeling test cases (Causin-case1 and Causin-case2).The results displayed in Figure 17C show that the values of ΔV predicted by the OMVS fall within the range of values measured experimentally in Reference 92 and theoretically predicted in Reference 65.In addition, the LC maximum displacement at baseline (109 μm with IOP ¼ 15 mmHg) is within the range of values reported by the experimental studies of Yan et al. 93 59-138 μm ð Þ .Figure 17B is an the example of the 3D spatial distribution LC displacement at baseline IOP ¼ 15 mmHg ð Þ , showing that the LC biomechanics is only mildly influenced by the presence of the CRA/CRV opening.Interestingly, this effect has been reported in a previous optic nerve head biomechanics study. 94

| Biomechanics of cornea, sclera, choroid, and retina
In this subsection we discuss the biomechanical simulations provided by the OMVS for other ocular tissues than the lamina cribrosa, and their relevance when compared with real data.
As a premise, we emphasize the fact that few works have been done in this direction.Indeed, the possibility to measure the stresses and the displacements experienced by retina, choroid, sclera and cornea in vivo is very challenging; on the other hand, many studies tried to replicate the behavior of such tissues with phantoms or ex vivo experiments.For this reason, we have compared the simulations results of the sclera and the cornea, which are the external and more accessible tissues, whereas for the biomechanics of the retina and choroid we plan to conclude this exhaustive analysis when more data will be available in the future.For all test cases we have used the 3D OMVS parameters presented in Table 3 and the following input data SP ¼ 126:3 mmHg DP ¼ 84:2 mmHg HR ¼ 60beats=min: First, we refer to the work by Myers et al., 95 where the authors are performing an inflation test on a posterior bovine sclera.This in vitro analysis exhibits a nonlinear response to controlled pressurization.We completed three virtual experiments with increasing IOP values of 14:78, 23:3 and 30 mmHg, respectively.Myers et al. reported that the sclera can be modeled as a quasi-linear elastic material in the physiological pressure range of 2 À 6kPa ≈ 15 À 45 mmHg ð Þ , behavior which is well replicated by our virtual experiments resumed in Figure 18A.The colormap of the 3D overview of the scleral displacement in Figure 18B shows a larger deformability of the optic nerve head area, coherently with the work of Myers et al. 95 Second, we refer to the correlation between axial length and intraocular pressure discussed in the work by Detorakis and Pallikaris. 96The axial length is the distance between the cornea and the back of the retina.It is the combination of the depth of the anterior chamber, the lens and the vitreous humor chamber.In adults, it measures between 22 and 25 mm. 97In this case, we have performed six virtual experiments with increasing IOP values of 1:5, 3:83, 10:5, 14:78, 23:3 and 30 mmHg, respectively.The authors infer a linear dependency between the measurements of intraocular pressure and axial length.The results provided by the OMVS (Figure 18C) are in good agreement with such a correlation between these two quantities and within the range of measurements in adults.
Third, we compare the model predictions of the OMVS with the findings of Boyce et al. 98 on the corneal displacement.In their study, the authors have conducted experiments on a bovine cornea to measure T A B L E 7 CRA flow comparison with experimental data using the hypothesis of CRA diameter of 160 μm. 87

References
Peak  The main goal of the current study was to design the Ocular Mathematical Virtual Simulator, a mathematical and computational framework able to simulate and predict the biomechanics and hemodynamics within the main tissues of the eye.Advanced and innovative numerical methods were developed to take into account (i) the combined effects of flows and different structures from a multiphysics perspective, and (ii) the interplay between local and systemic effects on the blood flow inherent to the multiscale structure of the problem.The proposed results provide an extensive validation of the model, that paves the way for its utilization in the context of clinical research for glaucoma.By comprehensively modeling the lamina cribrosa's and the surrounding tissues behavior under varying conditions, 16,17,91,100 our framework could aid in predicting and understanding the hemodynamical and biomechanical changes associated with glaucomatous damage.This, in turn, may contribute to improved diagnostic and treatment strategies, enhancing patient care.Some hypotheses, however, had been incorporated in the development of the OMVS, which lead to some limitations for its use.First, the retina has been considered as an elastic tissue from the biomechanics standpoint and represented by a 0D model for the hemodynamics.Alternatively, a porous media model could have been used to describe blood flow dynamics in this tissue, similarly to the lamina cribrosa.However, the internal retinal vasculature is more complex than the lamina cribrosa and it is not fully understood yet. 101Moreover, the images that can be obtained on the retina are mostly reduced to a small portion of the total tissue, usually the fundus, which is the part close to the optic nerve head.Second, the present configuration of the OMVS does not incorporate the effect of the vitreous humor and only accounts for a steady description of the aqueous humor, in which the IOP value is considered constant.As a further step, and thanks to the modularity of the OMVS design, other models developed to describe the vitreous humor, 102 the aqueous humor and the IOP modeling 103 could be integrated in the proposed framework.Third, further investigation into uncertainty quantification should be undertaken to comprehensively gauge the level of confidence associated with the obtained results and understand the robustness of the OMVS framework.Fourth, although our primary emphasis remains the investigation of glaucoma, it is conceivable that a comparable approach could be extended to provide insights into other ocular pathologies, including Age-Related Macular Degeneration, Diabetic Retinopathy, or other neurodegenerative diseases. 80Nevertheless, for the model to offer substantial insights into these diverse conditions, it would necessitate overcoming certain limitations outlined earlier.Furthermore, its clinical significance in those areas would necessitate dedicated research efforts and validation specific to those conditions.In conjunction with this, the incorporation of machine learning techniques also presents considerable promise in the advancement of sophisticated models for ocular diseases.It has the potential to enrich our comprehension of diverse ocular pathologies by meticulously analyzing extensive datasets, as recently reviewed for instance in Reference 104, see also References 105,106.Furthermore, the combined use of mathematical modeling and artificial intelligence may help progressing personalized research in ophthalmology.

APPENDIX A
In this appendix we describe the main steps of the HDG method utilized in the OMVS framework, based on Reference 15.We first introduce the general framework and its notations, then specify our approach to the hemodynamics in the lamina cribrosa (Equation ( 6)), and finally to the biomechanics described by systems (7) and (8).

A.1. | Introductory notation
In view of the finite element approximation, we partition Ω i with i ¼ sclera, choroid, retina, cornea, lc into the union of closed straight tetrahedra K, and we denote by T i,h -called triangulation-the collection of elements K such that Ω i ¼ [ K T i,h K.For each K T i,h , we denote by h K the diameter of K. We let h ≔ max K T i,h h K and we consider a family of conforming, regular triangulations of Ω i , T i,h f g h > 0 (see Reference 48 Chapter 3).For each element K T i,h we indicate by ∂K the boundary of K and by n ∂K the associated outward unit normal vector.The 3-dimensional measure of the element K is indicated with j K j while the 2-dimensional measure of each face of ∂K is denoted by j F j.
We let ℱ i,h denote the collection of all the faces of T i,h , whose union forms the skeleton of the decomposition T i,h .The set ℱ i,h naturally splits into the subset ℱ ∂Ω i i,h of faces belonging to ∂Ω i , and into the subset of faces belonging to the interior of Ω i , denoted by ℱ 0 i,h .Finally, assuming that the decomposition T i,h is such that for all faces F in ℱ i,h with with ϕ ?L 2 ℱ i,h ð Þdefined as and with where Þ denotes the space of polynomials of degrees less or equal that k on K (resp.F).Definitions (A1a), (A1b) and (A2) imply that functions belonging to V i,h and W i,h are, in general, discontinuous across elements faces of T i,h , while functions in M i,h are discontinuous across element edges of ℱ i,h ∖ ℱ Γ i,I i,h , single-valued on each face F ℱ i,h of the skeleton of T i,h and constant on Γ i,I .The functions of M i,h play the role of "connectors" between adjacent elements enabling the possibility to use static condensation to reduce the computational cost.To simplify the presentation, we introduce the compact notations for the integrals:

. | Porous media model
We now consider the elliptic boundary value problem of second order described in Equation (6).At the discrete level, in Ω lc,h the numerical normal flux j h is defined as where the quantity τ is a non-negative stabilization parameter.The HDG formulation then reads: find j h V lc,h , p h W lc,h and b p h M lc,h such that. where and h N and h D are the associated Neumann and Dirichlet boundary conditions described in Section 2.4.The dependent variables j h and p h are the approximations of j and p in the interior of each element K T lc,h , whereas the dependent variable b p h is the approximation of the trace of p on each face of ℱ lc,h .The numerical normal flux (A3) is characteristic of a particular class of HDG methods, the so-called Local Discontinuous Galerkin Hybridizable (LDG-H) methods proposed and investigated in a series of seminal papers (see Reference 107 and references therein).

A.3. | Linear elasticity model
In this paragraph we apply the HDG method to the linear elasticity equations (Equations ( 7) and ( 8)).
In addition to the notation introduced in the previous part, we define on Ω i for i ¼ sclera, choroid, retina, cornea, lc the stabilization parameter τ and the L 2 -orthogonal projection from L 2 E h ð Þ onto M h , denoted P M such that: The use of P M is needed in order to obtain the optimal accuracy; without it, the formulation might result slightly less accurate. 108Thus, we write the HDG formulation for linear elasticity: find σ el,i,h V i,h , u i,h W i,h and u i,h M i,h such that 8i ¼ sclera, choroid, retina, cornea, lc and 8 v i , w i , μ i,1 ,μ i,2 , μ i,3 Aσ el,i,h , v where n o with S ¼ set of symmetric tensors in ℝ 3Â3 È É , A defined as and g N and g D are the associated Neumann and Dirichlet boundary conditions described in Section 2.4.Observe that (i) the spaces introduced in this paragraph are the natural extension of the spaces needed for the porous media system, and (ii) the symmetry for elements v i V i,h is enforced in an essential manner in their definition. 108

3
Scheme of the blood supply and drainage for the lamina cribrosa and retina employed in our model.F I G U R E 4 Circuit representing the circulation in the lamina cribrosa, the retina and the central retinal vessels.

F I G U R E 5
Poroelastic model framework for the lamina cribrosa.
Micro-structure F I G U R E 6 Hemodynamic model framework for the lamina cribrosa.

2
Process to generate a computational mesh from a CAD drawing using the Salome platform.T A B L E 1 Geometrical parameters that can be modified by the user during the generation of the computational mesh.max size of the mesh grid h for the lamina cribrosa sub-domain hsize_eye 1.0 max size of the mesh grid h for all the eye except the lamina cribrosa subthickness proportion with respect to the original CAD thickness

F I G U R E 1 4
Color Doppler Imaging (CDI) blood velocity profile within the CRA.Courtesy of the Eugene and Marilyn Glick Eye Institute (Indianapolis, USA).F I G U R E 1 5 Overview of the ocular mathematicalvirtual simulator validation: comparison with experimental data.

(
A) CRA/CRV blood flow.(B) Lamina cribrosa perfusion and pressure computed on Γ .(C) 3D hemodynamics in the lamina cribrosa.F I G U R E 1 6 Ocular mathematicalvirtual simulator validation study: hemodynamics of the CRA/CRV and lamina cribrosa.Input data retrieved from Reference 82.T A B L E 6 Input data for the comparative study with Reference 82.
the displacement under constrained inflation conditions.We conducted six virtual experiments with increasing IOP values of 1:5, 3:83, 10:5, 14:78, 23:3 and 30 mmHg, respectively.Boyce et al. found that, if the IOP applied on the cornea is in physiological ranges, the displacement has a linear behavior (Figures 7, 9 and 10 of Reference 98).Our model predicts a similar behavior for the human cornea, as illustrated in Figure18D.This linear behavior is well captured by the OMVS, despite the presence of several nonlinearities in the biomechanical interaction between the ocular tissues and in the hemodynamics of the ocular posterior segment.Furthermore, Figure18Dreports a maximum corneal displacement up to 0:03 cm, which is consistent with the values reported by Elsheikh et al.99 3D view of the displacement spatial distribution with IOP= 15mmHg.(C) Comparison among the IOP-induced increments of LC volume below reference (Δ ) determined experimentally by 88 , numerically by 65 , and the numerical predictions of the OMVS.(A) Deformed 3D geometry of the LC using 5 different values of IOP.F I G U R E 1 7 Ocular mathematicalvirtual simulator validation study: biomechanics of the lamina cribrosa.

(A )
Scleral displacement (B) Sclera 3D overview (C)Axial length displacement (D) Corneal displacement F I G U R E 1 8 Ocular mathematicalvirtual simulator validation study: biomechanics of the sclera and cornea.Dotted lines represent the ideal linear behavior.
parameters for the ocular mathematical virtual simulator.Summary of the data for the resistors in the 0D hemodynamics model.
Summary of the data for the capacitors in the 0D hemodynamics model.