Modeling the Origin of the Ocular Pulse and Its Impact on the Optic Nerve Head

Abstract

Purpose: To use finite element (FE) analysis to understand the origin of the ocular pulse and predict its biomechanical impact on the optic nerve head (ONH).

Methods: An FE model of a healthy eye was reconstructed. The choroid was biphasic and consisted of a solid phase (connective tissues) and a fluid phase (blood). We applied arterial blood pressure at 18 entry sites (posterior ciliary arteries) and venous blood pressure at 4 exit sites (vortex veins). For one cardiac cycle, we reported the resulting pulse volume, the ocular pulse amplitude (OPA), and diastole-to-systole ONH deformations. We also studied the effect of a change in scleral stiffness, in arterial pressure, and in baseline IOP.

Results: During the cardiac cycle, a change in arterial pressure resulted in choroidal expansion, which in turn induced a change in IOP (the OPA) and ONH deformations. From diastole to systole, we found that choroidal expansion made the peripapillary retina move anteriorly, but both choroidal expansion and the OPA made the prelamina and LC move posteriorly. The net result was shearing of neural tissues in the neuroretinal rim. Both a stiffer sclera and a higher IOP resulted in a higher OPA, smaller pulse volume, larger diastole-to-systole ONH strains, and neural tissue shear in the neuroretinal rim. Increasing the arterial pressure had the same effect, except that it increased the pulse volume.

Conclusions: Our models indicate that, during the cardiac cycle, the OPA and choroidal expansion can deform the ONH with a net shearing of neural tissues within the neuroretinal rim.

Intraocular pressure (IOP) is the only manageable risk factor for glaucoma. However, IOP is not a constant value, but instead is pulsatile in nature. The difference between systolic and diastolic IOP has been defined as the ocular pulse amplitude (OPA) that ranges between 0.9 and 7.2 mm Hg in healthy subjects.¹ Generally, it has been well accepted that the OPA is mainly caused by acute choroidal expansion due to the pulsatile blood flow. 

Although the OPA has been suggested to play a role in glaucoma,^2–4 the exact relationship between OPA and visual field loss and with other glaucoma risk factors has remained elusive. On the one hand, the vascular mechanisms suggest that a smaller OPA may indicate abnormal or insufficient ocular blood supply.⁵ This has been suspected as one of the underlying mechanisms of glaucoma, and in fact a decreased OPA has been observed in normal-tension glaucoma subjects.^3,4 On the other hand, a larger OPA should yield larger stress/strain fluctuations within the optic nerve head (ONH), particularly within the lamina cribrosa (LC), and this may have an impact on retinal ganglion cell axons. LC displacements (between 2 and 9 μm) induced by the ocular pulse have been measured in vivo by low-coherence tissue interferometry and phase-sensitive optical coherence tomography.^6–8 In addition, studies have shown that a higher OPA was usually observed at a higher IOP level^1,9 and with older age in glaucoma patients.¹⁰ This may be consistent with the fact that a higher OPA has been reported to be associated with some types of glaucoma including chronic angle-closure glaucoma and suspected open-angle glaucoma subjects,¹⁰ although it should be emphasized that some of these trends (OPA versus age) were not observed in healthy subjects in another study.¹ 

Because of conflicting evidence, it has remained unclear how the measurement of the OPA (currently assessed with dynamic contour tonometry¹ and pneumotonometry¹¹) could become useful in the clinical management of glaucoma, and whether a higher OPA would have a negative impact on the ONH. To better understand the origin and effects of the OPA on the ONH, we suggest that it is critical to develop models of the eye. While several computational models of the eye have been proposed^12–18 and one considered choroidal expansion (Feola A, et al. IOVS 2017;58:ARVO E-Abstract 3153), none have yet tried to describe the origin of the OPA. Other analytical models^19,20 that can describe the OPA are typically too simple to capture the complex interactions between the OPA and ONH deformations. 

The aim of this study was to use finite element (FE) modeling to better understand the origin of the OPA and its biomechanical impact on the ONH. In addition, FE sensitivity studies were performed to better understand how the OPA is affected by a change in scleral stiffness, ophthalmic artery pressure, and IOP, as those parameters have been shown to affect the OPA clinically. 

Methods

In this study, we used FE modeling to mimic choroidal expansion (by applying arterial and venous blood pressures) and estimated the resulting OPA, ONH deformations, and pulse volume. We further studied the effect of a change in scleral stiffness, ophthalmic artery pressure, and IOP. 

Three-Dimensional Geometry of the Ocular and Orbital Tissues

Our three-dimensional (3D) eye model was adapted and modified from our previous study (Fig. 1).^12,13 In brief, the optic nerve and eye globe were reconstructed from magnetic resonance imaging images of a healthy subject. The corneoscleral shell was assumed to be spherical (outer diameter: 24 mm; thickness: 1 mm), and the optic nerve (circular cross section) consisted of three tissues: the nerve tissue (diameter range, 3.0–3.88 mm; length, 24.8 mm), the pia mater (thickness, 0.06 mm), and the dura mater (thickness, 0.3 mm). We used a generic ONH geometry that was embedded within the corneoscleral shell and that incorporated the scleral flange (length: 0.4 mm; thickness: 0.45 mm), Bruch's membrane (thickness: 5 μm; uniform over the globe), the choroid (thickness: 134 μm; uniform over the globe; initial volume: 213 μL), the LC (central thickness: 0.28 mm; anterior and posterior radii of 2.8 and 2.6 mm, respectively), the prelaminar tissues (thickness: 0.2 mm), and the border tissues of Elschnig and Jacoby as extensions of the pia matter (thickness: 0.06 mm). 

Figure 1

 

Reconstructed geometry and FE mesh of the whole eye model. Detailed ONH structures were included (LC, sclera, neural tissues, pia, dura, choroid, BM, and border tissue) using average measurements from the literature. Bottom right: a zoomed image of the ONH region showing BM, the choroid, and the border tissue.

The reconstructed model was then discretized into a hexahedron-dominant mesh with 67,584 eight-node hexahedra and 3024 six-node pentahedra using ICEM CFD (ANSYS, Inc., Canonsburg, PA, USA; Fig. 1). The mesh density was numerically validated through a convergence test. No symmetry conditions were applied. 

Biomechanical Properties of the Reconstructed Eye Tissues

The sclera was modeled as a fiber-reinforced composite as described in our previous paper.²¹ The LC, neural tissue, and Bruch's membrane (BM) were modeled as isotropic elastic materials and thus described with a single stiffness value.¹⁴ The pia and dura were modeled as Yeoh materials, derived from experimental data in porcine eyes.^12,13 The choroid was modeled as a vascular biphasic structure consisting of a solid phase (connective tissues) and a fluid phase (blood). Since the elastic modulus of the BM–choroidal complex varies between 0.32 and 0.88 MPa (average: 0.45 MPa)²² and that of BM alone was found to be 10.79 MPa (Chan WH, et al. IOVS 2007;48:ARVO E-Abstract 2187), we estimated that the elastic modulus of the choroid was 0.06 MPa (for a given thickness of 134 μm). The Poisson's ratio of the choroid was set to 0 to allow it to expand during the cardiac cycle. We chose a solid volume fraction of 0.55 for the choroid as the choroidal blood volume fraction has been estimated between 0.30 and 0.60.²³ The hydraulic permeability of the choroid was 45,037 mm²/MPa.s, as calculated from the choroidal vascular resistance,^24,25 for an average choroidal vessel diameter²⁶ of 70 μm and a blood viscosity²⁷ of 3.4 mPa.s. All biomechanical parameters are listed in Table 1. 

Table 1

 

Tissue Biomechanical Properties Used for the Baseline Model

 

Boundary Conditions

For our baseline FE model, we used two sets of boundary conditions to ensure numerical stability. First, the orbital apex of the optic nerve was fixed to mimic its connections to the optic canal by fibrous adhesions.¹² Second, the corneoscleral shell was fixed near the equator on two opposite sides using four nodes on each side to secure the eye globe (Fig. 2a). 

Figure 2

 

(a) Half FE geometry and boundary conditions: fixed optic nerve at the orbital apex and fixed sclera nodes near the equator on two opposite sides. (b) Choroidal layer illustrating the nodes where arterial and venous blood pressures were applied (one node per artery and per vein). In particular, the venous blood pressure was applied at the vortex veins (one in each quadrant) and the ophthalmic artery pressure was applied at the short and long PCAs. The color map represents the choroidal blood pressure (higher near the posterior ciliary arteries and lower near the vortex veins). (c) In the choroid, blood can flow from the posterior ciliary arteries (entry of blood in the choroid) to the vortex veins (exit of blood). The direction of blood flow within the choroidal layer is illustrated with blue arrows. (d) ONH region illustrating the different pressure fields (IOP, CSFP) applied in each model and (e) enlarged view of the ONH illustrating the measurement of LC depth.

Loading Conditions

We applied a baseline IOP of 15 mm Hg to the inner limiting membrane, and a baseline cerebrospinal fluid pressure (CSFP) of 11.3 mm Hg within the arachnoid space.²⁸ For simplicity, blood vessels and blood flow were not modeled within the sclera; that is, the sclera did not contain holes to let the blood vessels pass through. However, we assumed that blood could still enter and exit the choroid by applying prescribed fluid (blood) pressure at nodes directly located on the posterior surface of the choroid (Figs. 2b, 3). The nodes where blood can enter the choroid represent the posterior ciliary arteries, and the nodes where blood can exit the choroid represent the vortex veins. We used 16 nodes approximately 2.3 mm away from the dural sheath to represent the entry sites of the short posterior ciliary arteries (PCAs). Two additional nodes, one temporally and one nasally, were used to represent the entry sites of the long PCAs. Note that each node is meant to represent one artery. The long PCAs mainly supply a segment of the iris, ciliary body, and peripheral choroid on each side, respectively, whereas the short PCAs provide blood supply for both the ONH and the choroid.^29,30 We applied the same arterial blood pressure at these 18 nodes, and such pressure was allowed to vary between 70.8 and 93 mm Hg to mimic diastolic and systolic changes of the ophthalmic artery pressure.³¹ In addition, a constant blood pressure of 15 mm Hg was applied at four nodes where the vortex veins leave the choroid (typically one in each quadrant) to drain the blood out of the choroid.³² It should be emphasized that it is the pressure difference between the PCA entry sites/nodes and the vortex vein exit sites/nodes that drives the entire choroidal blood flow. An illustration of choroidal blood flow streamlines originating from the PCA entry sites/nodes and draining into the vortex vein exit sites/nodes is shown in Figure 2c. Each cardiac cycle lasted 1 second, and we ran two cycles for each model. 

Figure 3

 

3D schematic showing the posterior portion of the eye globe and the blood vessels (vortex veins and posterior ciliary arteries) that pass through the posterior sclera and connect with the choroid (as described in the literature³⁰).

Modeling the Origins of the OPA

To better understand how IOP was applied and changed during our simulations, one may consider IOP as having two components: (1) a baseline IOP that was applied to the inner limiting membrane as an input in an initial FE step (as in any other FE simulation of the eye), plus (2) a fluctuating IOP component as an output (i.e., the OPA) that was the direct result of choroidal expansion in a subsequent FE step. In our model, the total IOP at any moment during the cardiac cycle is the sum of these two components. 

How is the second IOP component (i.e., the OPA) generated? First, it is important to realize that we assumed that the vitreous body (filling the entire eye in our model) was incompressible. This is a reasonable assumption because the vitreous mostly consists of water, which is also incompressible. Therefore, the volume of the vitreous body was constrained to remain constant during choroidal expansion in all our simulations. Since our model aims to reproduce choroidal expansion during the cardiac cycle (due to the pulsatile choroidal blood flow), such an expansion will try to act to deform and change the volume of the vitreous body. Since the vitreous body was constrained as incompressible, an internal pressure needs to be applied to maintain the volume of the vitreous body. This pressure term is an output of our model and can be understood as the OPA. 

Our approach can be simply compared to a scenario in which a water-filled balloon is being poked. When the water-filled balloon is poked from outside, the enclosed water pressure increases in order to maintain the water volume and counterbalance the external perturbation. 

In an FE solver, the pressure that enforces the volume constraint is a Lagrangian multiplier and can be estimated through an augmented Lagrangian method (see Equation 1):  

\begin{equation}\tag{1}{p_{k + 1}} = {p_k} + \varepsilon *\left( {V - {V_0}} \right){\rm }\end{equation}

where Display Formula\(p\) is the fluid pressure, k is the augmentation iteration, Display Formula\(\varepsilon \) is a user-defined penalty factor, V and V₀ are the current (i.e., final) and initial volume of the vitreous body, respectively. This pressure is applied to the entire enclosed surface of the constrained volume and it is updated until the current (i.e., final) volume matches the original volume. Hence, this pressure represents the pressure change (i.e., the OPA) required to counterbalance the external perturbation (i.e., choroidal expansion) to maintain the volume of the vitreous body.  

FE Sensitivity Studies: Parameters Affecting the OPA

We aimed to understand whether scleral stiffness, ophthalmic artery pressure, and IOP had an impact on the OPA. For simplicity, one parameter was varied at a time and ranged from 80% to 120% of the baseline parameter values (scleral stiffness and ophthalmic artery pressure). A total of seven FE models were run (one baseline model at IOP = 15 mm Hg, two for low and high scleral stiffness, two for low and high ophthalmic artery pressure, and two for elevated IOPs of 30 and 45 mm Hg). 

FE Processing to Predict the OPA, Choroidal Pulse Volume, and ONH Deformations

All FE models were solved using FEBio v2.6.3 (Musculoskeletal Research Laboratories, University of Utah, Salt Lake City, UT, USA). To estimate the OPA in all models, a “volume constraint” was imposed on the inner limiting membrane and the augmented Lagrangian was turned on. The tolerance and penalty factor were set as 0.00001 and −0.001, respectively. Those values were chosen because they maintained the volume of the vitreous within 0.000001% during the cardiac cycle. All pressure estimates (to maintain the volume of the vitreous body constant) were extracted from the FEBio output log file (.log) under the section “Volume Constraint” at each converged step. 

Each model was run in two steps: the initial step and the volume constraint step during which the OPA was generated (Fig. 4). During the initial step, we applied boundary conditions and initial loading conditions, that is, the baseline IOP, CSFP, ophthalmic artery pressure, and vortex vein pressure. During the second step, we imposed a volume constraint on the inner limiting membrane to ensure that the vitreous volume remained constant. In addition, the CSFP and vortex vein pressures were kept constant; the ophthalmic artery pressure was allowed to vary from diastole to systole periodically, and we reported the resulting change in IOP (i.e., the OPA). Note that elevating the baseline IOP to 30 or 45 mm Hg would affect the conditions at the end of step 1 by making the sclera stiffer and the eye volume slightly higher, which is consistent with reality. Under such a scenario, we should expect a different OPA for a different given baseline IOP. 

Figure 4

 

Profiles of different pressure loads.

The results were analyzed using MATLAB (v2015b; MathWorks, Natick, MA, USA). For each model, we reported the resulting OPA, the pulse volume, the diastole-to-systole displacements of the LC and of the prelamina (radial component perpendicular to the corneoscleral shell surface; positive: outward the eye; negative: inward the eye), the change in LC depth, and the average first and third principal strains in the LC and prelamina. The pulse volume was calculated as the change in choroidal volume during one cardiac cycle. LC depth was defined as the distance from the center of the anterior LC surface to the plane passing through the anterior LC boundary (Fig. 2e). Finally, we also estimated the ocular rigidity for each model based on the Friedenwald equation (Equation 2)²⁰:  

\begin{equation}\tag{2}Ocular \ Rigidity = {{log(IO{P_1}) - {\rm{log}}(IO{P_2})} \over {\Delta V}}{\rm }\end{equation}

where IOP₁ and IOP₂ are the systolic and diastolic IOP, respectively, and Display Formula\(\Delta V\) is the choroidal volume change (or pulse volume) during one cardiac cycle.  

Contributions of OPA and Choroidal Expansion to ONH Deformations

During the cardiac cycle, there are potentially two contributors that could induce LC deformations: choroidal expansion and the OPA. To better understand the contribution of each to LC deformations, we ran an eighth FE model that simulated choroidal expansion but not the OPA. This was achieved by removing the volume constraint imposed on the inner limiting membrane in the baseline model. 

Results

Prediction of OPA and Pulse Volume During the Cardiac Cycle

In the baseline FE model, a change in arterial pressure (from 70.8 to 93 mm Hg) during the cardiac cycle resulted in choroidal expansion (diastolic thickness: 135 μm; systolic thickness: 136 − 142 μm [anterior to posterior]; resulting pulse volume: 3.81 μL) that in turn induced a change in IOP (OPA: 2.27 mm Hg; diastolic IOP: 13.90 mm Hg; systolic IOP: 16.16 mm Hg). IOP and pulse volume are shown in Figure 5 as a function of time for two cardiac cycles. 

Figure 5

 

IOP and choroidal volume over two cardiac cycles in the baseline model.

Diastole-to-Systole LC Displacements, Strains, and Depth Changes

During the cardiac cycle, we found that the ONH “pulsed” and moved posteriorly (diastole to systole) and anteriorly (systole to diastole; Figs. 6a, 6b). The diastole-to-systole radial displacement of the central anterior LC point was 7.81 μm (posterior) and that of the central anterior prelamina point was 7.23 μm (Fig. 6c). We also found that the averaged strains generated within the LC and prelamina at systole were relatively small and less than 1% (first principal strains in the LC and prelamina: 0.12% and 0.17%, respectively; third principal strains in the LC and prelamina: −0.18% and −0.16%, respectively). Note that the strains were computed with diastole as the reference state. In addition, LC depth increased by 2.71 μm from diastole to systole. 

Figure 6

 

Cross section of the deformed eye at (a) diastole and (b) systole. Displacements were exaggerated 15 times to ease visualization. The highlighted points (red: center of prelamina; black: center of LC) were tracked and (c) their radial displacements are shown over two cardiac cycles.

Diastole-to-Systole Shearing of the Neuroretinal Rim

From diastole to systole, we observed shearing of neural tissues at the neuroretinal rim. Specifically, in the baseline model, the prelamina and the peripapillary retina moved in opposite directions (Figs. 7a, 7b); the prelamina moved posteriorly (7.23 μm), and the peripapillary retina moved anteriorly (−2 μm). The difference between these two displacements at systole was used to define the amount of neuroretinal shear (here: 9.23 μm). This resulted in strain ring patterns (of relatively large magnitude) within the neuroretinal rim (max first principal strain: 1.37%; max third principal strain: −1.46%; Figs. 7c, 7d). 

Figure 7

 

(a) Cross section of the ONH. (b) The radial displacement of the highlighted points (red: center of prelamina; black: a point on the retinal surface located 0.65 mm away from Bruch's membrane opening) was plotted over two cardiac cycles. (c, d) First and third principal strain color maps showing a circular pattern of large strains in the neuroretinal rim.

Contributors to LC Deformations During the Cardiac Cycle

We found that the OPA deformed the ONH during the cardiac cycle but was not the only contributor; choroidal expansion also deformed the ONH, but its contribution was lower. Specifically, in the FE model with choroidal expansion but no OPA, LC strains were 0.03% (first principal) and −0.04% (third principal); in the FE model with choroidal expansion and OPA, LC strains were 0.12% (first principal) and −0.18% (third principal; Fig. 8). Hence, on average, choroidal expansion contributed to 23% of LC strains, while the OPA contributed to the rest. We also found that the change in LC depth (diastole to systole) was 2.71 μm in the baseline model and 0.66 μm after removing the OPA. 

Figure 8

 

LC principal strains (left: first; right: third) resulting from choroidal expansion and the OPA (red), or from choroidal expansion alone (black).

Effect of Scleral Stiffening on OPA and Diastole-to-Systole ONH Deformations

We found that stiff scleras increased the OPA but reduced the pulse volume. Stiff scleras also increased diastole-to-systole LC strains, LC and prelamina displacements, LC depth changes, and the amount of neuroretinal shear (Table 2). We also found that the ocular rigidity was higher for a stiffer sclera. 

Table 2

 

Effect of Scleral Stiffening From a Low (−20% Baseline) to a High (+20% Baseline) Value

 

Effect of Ophthalmic Artery Pressure on OPA and Diastole-to-Systole ONH Deformations

We found that both the OPA and the pulse volume increased with increasing ophthalmic artery pressure. ONH displacements/strains and the amount of shearing in the neuroretinal rim were larger with an increase in ophthalmic artery pressure (Table 3). 

Table 3

 

Effect of Increasing the Ophthalmic Artery Pressure From a Low (−20% Baseline) to a High (+20% Baseline) Value

 

Effect of IOP Elevation on OPA and Diastole-to-Systole ONH Deformations

We found that an elevated baseline IOP increased the OPA, but decreased the pulse volume. Both the amount of shearing in the neuroretinal rim and diastole-to-systole LC strains increased with a higher baseline IOP (Table 4). 

Table 4

 

Effect of Increasing the Baseline IOP From 15 mm Hg to 30 and 45 mm Hg

 

Discussion

In this study, we used FE modeling to better understand the origin of the ocular pulse and its biomechanical impact on the ONH. Our models predicted that a change in arterial pressure (diastole to systole) resulted in choroidal expansion,^33,34 which in turn induced a change in IOP. It was also found that both choroidal expansion and the OPA contributed to deform the LC and the prelamina during the cardiac cycle with a characteristic shearing of neural tissues in the neuroretinal rim. Changes in scleral stiffness, ophthalmic arterial pressure, and IOP affected the OPA and pulse volume, as has been hypothesized/observed clinically. 

A Change in Arterial Pressure Resulted in a Change in IOP

In our baseline FE model, we found that a change in arterial pressure during the cardiac cycle resulted in choroidal expansion, which in turn induced a change in IOP. In other words, we were able to model the origin of the ocular pulse with FE, and our predicted OPA (2.27 mm Hg) and pulse volume (3.81 μL) were consistent with those measured experimentally (OPA: 0.9–7.2 mm Hg among 148 subjects,¹ 0.71–3.09 mm Hg at IOP = 15 mm Hg,⁹ pulse volume: 2.61–8.74 μL at IOP = 15 mm Hg⁹). Clinically, DCT and pneumotonometry are the most commonly used tools to assess the OPA. Pneumotonometry is able to measure the IOP continuously by applying a column of flowing gas on the corneal surface. It defines OPA as the difference between the lowest and highest point of the pulse wave and typical OPA measurement ranges from 2 to 3 mm Hg,^35,36 which are in close agreement with our predicted result. Similarly, DCT measures the ocular pulse wave by recording the IOP continuously and defines the OPA as the difference between the minimum and maximum of the pulse wave contour.¹ DCT is a relatively new technology for noninvasive IOP measurement using the principle of “contour matching,” which is less dependent on the effect of individual corneal properties.³⁷ Similarly, our choroidal thickness changes (7.4 μm) were consistent with those observed experimentally (16.7 ± 10.9 μm).³⁴ To the best of our knowledge, while several models of the ONH have been proposed,^{12,14,15,18,38,39} none have yet reproduced both choroidal expansion and the OPA. 

Diastole-to-Systole ONH Displacements and Strains Were Small

During the cardiac cycle, we found that the ONH “pulsed” and moved posteriorly (diastole to systole) and anteriorly (systole to diastole). Radial LC and prelamina displacements in our baseline model were 7.81 and 7.23 μm, respectively. Such displacements are consistent with experimental measurements in human ONH tissues that ranged between 2 and 9 μm.⁶ We also found that the averaged strains generated within the LC and prelamina were relatively small and less than 1%. To the best of our knowledge, no studies have yet reported diastole-to-systole LC strains. Several studies have measured IOP-induced LC strains ex vivo^40–42 (between 2.8% and 8% for IOP changes between 5 and 40 mm Hg), and in vivo^43,44 (between 6.41% and 8.6% for IOP changes between 12 and 21 mm Hg). Our reported strain values were found to be considerably smaller, and this was expected, given the small OPA value (2.27 mm Hg). 

Both Choroidal Expansion and the OPA Contributed to Diastole-to-Systole ONH Deformations

Interestingly, we found that both choroidal expansion and the OPA contributed to diastole-to-systole ONH strains. We were able to tease out such a contribution by running a model in which choroidal expansion could occur but not a change in IOP. On average, we found that choroidal expansion contributed to 23% of the total diastole-to-systole LC strains, whereas the remaining contribution was attributed to the OPA alone. We believe that, as the choroid thickens, it can push and bend the border tissues of Elschnig and Jacoby, which may in turn compress the LC and deform it posteriorly. It is therefore possible that the morphology and biomechanical properties of the border tissues (currently unknown) may play important roles in deforming the LC during the cardiac cycle. 

Our Models Predicted Shearing of Neural Tissues in the Neuroretinal Rim

From diastole to systole, we found that choroidal expansion made the peripapillary retina move anteriorly, but the OPA made the prelamina and LC move posteriorly. The net result was shearing of neural tissues in the neuroretinal rim—the region where we observed the maximum values for the first and third principal strains (1.37% and −1.46%, respectively). This shearing amount was positively correlated with the maximum shear strain in the prelamina (see Supplementary Material). While not explicitly reported, this phenomenon (anterior retina movement together with posterior prelamina movement) has already been observed in vivo in humans using low-coherence tissue interferometry.⁶ Note that in the neuroretinal rim region, several groups have measured the Bruch's membrane opening (BMO) minimum rim width (MRW),^45–49 characterized as the maximum aperture at the level of BMO through which retinal ganglion cell axons can pass.⁴⁵ Compared to other conventional rim parameters, BMO-MRW was found to exhibit a higher diagnostic accuracy for glaucoma and a stronger association with visual field parameters.^46–48 BMO-MRW was also found to significantly decrease with age,⁴⁵ and was found to be smaller in glaucoma subjects.^46,47 It is not yet known whether shearing of neural tissues in the neuroretinal rim (as observed herein) could possibly induce axonal damage and be responsible for a decrease in BMO-MRW with age and glaucoma. While it is generally accepted that the main site of axonal damage is the LC,^50,51 nothing excludes the possibility that axonal damage could also occur in the neuroretinal rim as it is a sensitive region where retinal ganglion cell axons perform a sharp turn to enter the disc. This also fits with clinical observations on angle-closure suspects who showed a widening and deepening of the optic cup, decrease in neuroretinal rim width in particular in the temporal region, and thinning of the LC after a darkness-induced IOP increase of >15 mm Hg.⁵² 

Stiffening the Sclera Increased the OPA, Diastole-to-Systole LC Strains, and Neural Tissue Shear, but Reduced the Pulse Volume

We found that increasing the stiffness of the sclera resulted in an increase in OPA. Since it is well known that the sclera gets stiffer with age,^21,53,54 our results are consistent with the clinical observation that OPA is higher in older subjects from a glaucoma population.¹⁰ 

Interestingly, we found that a stiffer sclera resulted in larger diastole-to-systole LC strains and a larger amount of neural tissue shear in the neuroretinal rim. This result may appear counterintuitive at first. Using computational modeling, Sigal et al.^14,38,55 have previously reported that scleral stiffness was the major determinant of IOP-induced LC strains. In other words, a stiffer sclera should reduce IOP-induced LC strains. This has been confirmed experimentally in porcine eyes by two separate studies in which the peripapillary sclera was stiffened through cross-linking agents.^56,57 It is important to emphasize that our data do not contradict such findings, as we reported diastole-to-systole LC strains originating from the OPA and choroidal expansion, but not LC strains derived from a baseline IOP as was performed in those studies. Since the OPA increased with a stiffer sclera, it would seem logical to also observe an increase in diastole-to-systole LC strains. Overall, our work, and that of others, suggests that a stiff sclera is beneficial to protect the LC from deformations arising from an elevated baseline IOP, but a stiff sclera may be detrimental to protecting the LC and neural tissues from an increased OPA. We also found that a stiffer sclera reduced the pulse volume, which may result in abnormal choroidal blood supply in stiffer eyes; for instance, decreased pulsatile ocular blood flow has been reported in open-angle glaucoma and ocular hypertension.^58,59 It is yet unclear if any of these scenarios would facilitate the development and progression of glaucoma. 

In our models, an increase in scleral stiffness also resulted in an increase in ocular rigidity, and our predicted ocular rigidity values (range, 0.0125–0.0222 1/μL) were comparable to experimental measurements (range, 0.0112–0.0149 1/μL).⁶⁰ Clinically, an increased ocular rigidity has been reported in older subjects^60,61 and in subjects with primary open-angle glaucoma,^11,60 suggesting that ocular rigidity is representative of scleral stiffness. However, we also found that ocular rigidity decreased with increasing baseline IOP, suggesting that ocular rigidity does not represent only scleral stiffness. 

A Decrease in Ophthalmic Artery Pressure Reduced the OPA, the Pulse Volume, and Diastole-to-Systole ONH Deformations

Our models predicted that a smaller ophthalmic artery pressure resulted in a smaller pulse volume due to the reduced arteriovenous driving force. This latter can be assessed through the ocular perfusion pressure, typically defined as the difference between the ophthalmic artery pressure and IOP.⁶² It should be noted that a vascular mechanism has been proposed in the pathogenesis of glaucoma.^5,63–65 Vascular deficiencies including inadequate or unstable ocular blood supply can lead to ischemic damage or reperfusion injury to the optic nerve tissues and axons. Interestingly, glaucoma and high-risk ocular hypertensive patients were found to have reduced pulse volume when compared with normal subjects,^7,8 but it is not yet known if this is due to a decrease in ophthalmic artery pressure or an increase in scleral stiffness. According to our models, a decrease in ophthalmic artery pressure could be detrimental as it would limit blood supply to the choroid, but it could be beneficial if the goal is to reduce ONH deformations. A decrease in ophthalmic artery pressure reduced the OPA, diastole-to-systole LC and prelamina strains and displacements, LC depth changes, and the amount of tissue shear in the neuroretinal rim. It should also be noted that our models would predict the same trends (including a reduced OPA) if one were to increase the vortex vein pressure. This is because both a decrease in artery pressure and an increase in venous pressure would reduce the arteriovenous driving force (and thus the pulse volume). 

Increasing Baseline IOP Increased the OPA, Diastole-to-Systole LC Strains, and Neural Tissue Shear, but Reduced the Pulse Volume

Our models predicted that a higher baseline IOP (30 or 45 mm Hg) resulted in a larger OPA. This is not surprising because the sclera is nonlinear and stiffens with IOP. Since our model considered a nonlinear sclera, we were able to reproduce IOP-induced scleral stiffening, and thus an increase in OPA. This phenomenon has been observed in both population-based measurements¹ (0.12 mm Hg of OPA / 1 mm Hg of IOP) and in experiments where saline was injected into the eye to artificially increase IOP (0.05–0.073 mm Hg of OPA / 1 mm Hg of IOP).^9,66 

We also observed a decreased pulse volume for a higher baseline IOP, which could be due to a decreased ocular perfusion pressure. This has also been reported in experimental measurements (29% pulse volume decrease as IOP increased from 20 to 40 mm Hg).⁹ 

Finally, we found that increasing the baseline IOP increased neural tissue shear in the neuroretinal rim. This result is intuitive because an increase in IOP in our models will automatically stiffen the posterior sclera due to its stretch-induced (or nonlinear) biomechanical properties. 

A Framework to Assess ONH Biomechanics In Vivo Without Artificially Manipulating IOP

To assess the biomechanics of the ONH in vivo, one is typically required to artificially manipulate IOP (e.g., through ophthalmodynamometry) in order to assess the resulting displacements, strains, or changes in shape.^67–70 While such approaches may have value from a research point of view, they may not be easily translated clinically due to eye discomfort during testing (IOP elevation, sometimes large, needs to be maintained for several minutes). If one were to measure the OPA and diastole-to-systole choroidal expansion in vivo, such data could be combined with our FE models to derive the biomechanical properties of the ONH tissues. 

Clinical Implications: Effect of High-Frequency IOP Fluctuations on Glaucoma Pathogenesis

Very little is known about the implications of high-frequency IOP fluctuations in the development and progression of glaucoma. To date, the consensus is that a constant IOP level could be beneficial or detrimental for glaucoma patients. For instance, maintaining IOP to a constant low level is the mainstay of treatment in glaucoma patients. In addition, increasing IOP to a constant high level is the primary method to induce experimental glaucoma in most animal studies. From a purely mechanical point of view, IOP fluctuations at any timescale (seconds, days, or years) have the potential to injure the retinal ganglion cell axons in the ONH. Using a custom IOP telemetry system in monkeys, Downs et al.⁷¹ found that IOP could fluctuate significantly on short timescales because of the ocular pulse (fluctuations of 0.6–1.8 mm Hg) or because of eye blinks and saccades (changes in IOP up to 12 mm Hg). By examining a population of 183 glaucoma suspects, McMonnies⁷² suggested that such short-term IOP elevation episodes could have prognostic significance for glaucoma. Interestingly, pulsatile mechanical loading (on short timescales) may have a different effect on cell physiology when compared to constant steady loads. It has been shown that the response of cells to deformations (through a mechanism known as mechanotransduction) is dependent on both the magnitude and the rate of mechanical strain.^73,74 Some studies have shown that acute strain (short timescale) in neurons and neuron-like cells can lead to cellular injury.⁷⁴ There is also evidence (both in vivo and in vitro) showing that cyclic mechanical stress is more harmful to neurons than constant stress.⁷⁵ Overall, it is highly plausible that IOP fluctuations on short timescales could harm the retinal ganglion cell axons in the ONH, but further research is still needed to fully establish this link if it exists. 

Limitations

Eight limitations warrant further discussion. First, viscoelastic effects were not included. We would expect that ONH displacements and strains should change for different cardiac cycles (e.g., before and after exercising). We aim to take these effects into account in future models, as we believe it will become critical to understand ONH viscoelasticity in order to extract biomechanical properties from pulse data. 

Second, the choroid was simplified and modeled as a biphasic material, and consisted of a continuous porous solid matrix (connective tissues) mixed with a fluid (blood). Unfortunately, this material was unable to fully capture the complex microvasculature architecture of the choroid such as the segmental distribution of PCAs²⁹ and its autoregulatory capacity.⁷⁶ Furthermore, several simplifying assumptions were made. For instance, the artery pressure applied at the PCAs was taken from experimental measurements of the ophthalmic artery pressure in monkeys.³¹ For simplicity, we also assumed that the vortex vein blood pressure was constant and equal to IOP in our baseline model. However, it was found to be slightly higher than IOP in monkeys,⁷⁷ and to range between 10 and 15 mm Hg in rabbits.⁷⁸ In addition, blood was considered nonviscous as it is a limitation of the FEbio solver. However, blood viscosity is less important in our models since we did not model the microcapillary network of the choroid. Despite those limitations, it should be emphasized that our model was able to reproduce the physics of the ocular pulse. Future work should consider more complex permeability models to better describe the behavior of the choroid. 

Third, the permeability of the choroid may be a parameter with a large influence on the OPA. However, in this study, we only grossly approximated the permeability of the choroid using information about blood viscosity, the vascular resistance of the choroid, and the average vessel diameter of the choroid. To better understand the effect of permeability on the OPA, we performed additional simulations in which we varied the permeability from −80% to +80% of its baseline value (45,037 mm²/MPa.s). We found that increasing the permeability increased the pulse volume and the OPA (see Supplementary Material). This phenomenon can be simply understood through Darcy's law.⁷⁹ For a same pressure gradient driving blood flow (e.g., artery pressure minus vein pressure), a larger permeability will result in a larger blood flow, thus a larger pulse volume that will in turn generate a larger OPA. From our modeling estimates, a change in 20% permeability appears to be similar to a 2% to 12% change in ophthalmic artery pressure. This result suggests that the permeability of the choroid may be an important parameter affecting the OPA, and it would need to be better assessed experimentally. 

Fourth, our study did not account for regional variations in scleral thickness and elastic stiffness that are known to exist in humans and monkeys.^57,80 Such regional variations might possibly affect the resulting OPA and the measurements of ocular rigidity. Modeling studies that link regional variations in scleral stiffness/thickness to changes in OPA and ocular rigidity may be of interest to pursue. 

Fifth, our models excluded the circulation of blood within the LC and other pulsations including those from the central retinal artery, the retina, and the CSFP. It is highly plausible that the pulsation of the central retinal artery would also deform the LC. In addition, the location of the central retinal artery might affect such deformations, as it was recently shown that the position of the central retinal vessel trunk affected LC depth.^81,82 CSFP is also pulsatile in nature, and its pulsations are closely related to the retinal venous pulsations. Interestingly, CSFP pulsations are out of sync with the ocular pulse: CSFP experiences a peak slightly before IOP.^83,84 This difference in timing may lead to dynamic fluctuations of the translaminar pressure gradient, and this could potentially play a role in glaucoma.⁸⁵ The contribution of these pulsatile components may be important and should be considered in future models. 

Sixth, in our models we included a choroidal layer in the anterior part of the eye. While this is not physiologically accurate, it should be emphasized that most choroidal deformations during the cardiac cycle occur in the posterior part of the eye (since choroidal blood comes from the PCAs near the optic disc). Therefore, the effect of having a choroidal layer in the anterior part of the eye should be minimal. To verify this assumption, we performed an additional FE simulation (baseline model) in which the choroid and the retina were removed from the anterior chamber of the eye. We found that the OPA and the change in choroidal volume varied by less than 3% when compared to our original model. Hence, we believe that keeping choroidal and retinal layers in the anterior portion of the eye globe, albeit not realistic, should not have a significant effect on our results. Future studies may consider more realistic anterior chamber geometries while considering the flow of aqueous humor, as this may also have an impact on the OPA. 

Seventh, we reconstructed the eye and orbit geometry from a single eye only, and the ONH from average measurements in the literature. However, ONH and optic nerve geometries also vary across individuals⁵⁵ and, thus, might have significant influences on diastole-to-systole ONH strains. In addition, our model simplified the geometry of the retina, choroid, and BM with uniform thicknesses throughout the eye. Regional variations in tissue thickness may be taken into consideration in future studies. In the literature, we could not find volume measurements for the entire choroid to compare with our data. Volume measurements have typically been made in the macula region using optical coherence tomography. One could attempt to estimate the volume of the choroid using several anatomic landmarks. For instance, the mean subfoveal choroidal thickness was found to be 272 to 302 μm and is known to decrease nasally and temporally.^86,87 The mean peripapillary choroidal thickness was found to be 134 μm.⁸⁸ By using this information, we can estimate the volume of the choroid to be 219 μL, assuming that the total area covered by the choroid is similar to that of the retina (1094 mm²).⁸⁹ Birngruber et al.⁹⁰ estimated the volume of the human choroid to be 110 μL with several assumptions: a constant mean thickness of 200 μm that covers parts of a spheroid fundus with a diameter of 17.5 mm. In our study, the choroid was a uniform layer over the globe with an initial volume of 212.92 μL. After the initial step, the choroidal volume increased to 222.63 μL (for a mean ophthalmic artery pressure) and varied between 220.73 μL (diastole) and 224.53 μL (systole; see Fig. 5). The choroid geometry was adapted from our previous model, and the thickness of the choroid was chosen as the mean global parapapillary choroidal thickness measured by spectral-domain optical coherence tomography.⁸⁸ Future experimental studies should be considered to better estimate the total volume of the choroid as such measurements could be used to improve our models. 

Eighth, for simplicity, each artery or vein was represented by a single node due to its small size with respect to each element of the FE mesh. Technically, it is possible to assign more nodes for each vessel. We found that increasing the number of artery nodes increased the OPA, but increasing the number of vein nodes decreased the OPA (see Supplementary Material). This is consistent with the fact that increasing (or decreasing) the number of nodes can affect the blood flow. Using our approach (one node per vessel), we found that the predicted blood inflow velocity at the PCAs (systole: 7.3 cm/s; diastole: 5.0 cm/s) and the vortex veins (systole: 9.6 cm/s; diastole: 7.0 cm/s) was close to experimental results obtained using color Doppler imaging (PCAs: 9–19 cm/s at systole, 4–8 cm/s at diastole; vortex veins: 5–8 cm/s at systole, 2–5 cm/s at diastole).⁹¹ Hence, our approach may be a good start to model the dynamic behavior of choroidal blood flow. More complex loading conditions and/or choroidal blood flow models that also include regulation may be required to better match these experimental data. 

Conclusions

In this study, we modeled the origin of the ocular pulse and studied its impact on the ONH. We aimed to understand the links between ocular blood pressure and ONH biomechanics. We found that our models were able to reproduce the physics of the ocular pulse as observed clinically. Our models indicate that the OPA and choroidal expansion can deform the ONH with a net shearing of neural tissues within the neuroretinal rim. Future studies are needed to explore potential links with axonal loss in glaucoma. 

Acknowledgments

Supported by the Singapore Ministry of Education, Academic Research Funds, Tier 2 (MJAG; R-397-000-280-112). 

Disclosure: Y. Jin, None; X. Wang, None; L. Zhang, None; J.B. Jonas, None; T. Aung, None; L. Schmetterer, None; M.J.A. Girard, None