Morphometric, Hemodynamic, and Biomechanical Factors Influencing Blood Flow and Oxygen Concentration in the Human Lamina Cribrosa

Abstract

Purpose: We developed a combined biomechanical and hemodynamic model of the human eye to estimate blood flow and oxygen concentration within the lamina cribrosa (LC) and rank the factors that influence LC oxygen concentration.

Methods: We generated 5000 finite-element eye models with detailed microcapillary networks of the LC and computed the oxygen concentration of the lamina retinal ganglion cell axons. For each model, we varied the intraocular pressure (IOP) from 10 mm Hg to 55 mm Hg in 5-mm Hg increments, the cerebrospinal fluid pressure (13 ± 2 mm Hg), cup depth (0.2 ± 0.1 mm), scleral stiffness (±20% of the mean values), LC stiffness (0.41 ± 0.2 MPa), LC radius (1.2 ± 0.12 mm), average LC pore size (5400 ± 2400 µm²), the microcapillary arrangement (radial, isotropic, or circumferential), and perfusion pressure (50 ± 9 mm Hg). Blood flow was assumed to originate from the LC periphery and drain via the central retinal vein. Finally, we performed linear regressions to rank the influence of each factor on the LC tissue oxygen concentration.

Results: LC radius and perfusion pressure were the most important factors in influencing the oxygen concentration of the LC. IOP was another important parameter, and eyes with higher IOP had higher compressive strain and slightly lower oxygen concentration. In general, superior–inferior regions of the LC had significantly lower oxygen concentration than the nasal–temporal regions, resulting in an hourglass pattern of oxygen deficiency.

Conclusions: To the best of our knowledge, this study is the first to implement a comprehensive hemodynamical model of the eye that accounts for the biomechanical forces and morphological parameters of the LC. The results provide further insight into the possible relationship of biomechanical and vascular pathways leading to ischemia-induced optic neuropathy.

Open angle glaucoma (OAG) is an optic neuropathy characterized by a progressive loss of retinal ganglion cells (RGCs) that could be caused by excessive mechanical stress or other less well understood mechanisms, such as vascular deficiency. The primary risk factor for OAG is elevated intraocular pressure (IOP); however, the causal link between elevated IOP and RGC loss has yet to be clearly established. Indeed, up to a third of OAG patients do not have elevated IOP, and not all individuals with significantly elevated IOP (greater than 21 mm Hg) develop glaucoma.¹ Numerous sources of mechanical load, including IOP,² cerebrospinal fluid pressure (CSFP),³ and optic nerve traction,^4,5 have been speculated to cause RGC axonal death by inducing biomechanical insult to the ONH tissue. It is envisaged that these insults may cause direct mechanical damage to the RGC or/and disruption to the axonal transport of essential trophic factors.^6–8 Another possible pathway for RGC damage is vascular deficiency, and it is speculated that alterations in the structures of the optic nerve head (ONH), due to a modified biomechanical environment or poor ocular hemodynamics, could lead to ischemic RGC death.^9–12 

To date, both mechanical and vascular factors have been assumed to be involved in glaucoma pathogenesis, but the relationship between these factors remains unclear. In our previous work, we investigated the hemodynamical aspects of OAG by developing computational fluid dynamics (CFD) models of the lamina cribrosa (LC), a sieve-like network of collagenous connective tissue and microcapillaries that mechanically supports and nourishes the RGCs.¹³ We found that morphological parameters of the LC, particularly the increase in LC radius, had an adverse impact on LC blood flow and could lead to ischemia.¹² Given the technical limitations involved in experimentally characterizing blood flow in the LC microcapillary network, our computational model represented a powerful alternative to study the involvement of blood flow in glaucoma pathogenesis. However, in our previous study we did not include the elasticity of tissues, IOP and CSFP loads, or their potential impacts on blood flow. 

A few computational studies have emerged to study the effects of IOP on LC hemodynamics; for example, a CFD study by Carichino et al.¹⁴ suggested that an increase in IOP reduces the blood flow velocity in the central retinal artery, and a poroelastic model by Causin et al.¹⁵ suggested that an IOP increase could disrupt ONH blood flow. However, no studies have modeled the detailed morphometric features of the LC microcapillary network. With recent experimental data by Brazile et al.,¹⁶ suggesting that IOP could affect the morphology of the LC microcapillary and with the availability of improved imaging techniques (e.g., photoacoustic imaging) capable of visualizing the LC microcapillary in vivo,¹⁷ we believe that it is of great importance to study the influence of biomechanical forces on the LC microcapillary hemodynamics to complement the experimental data. 

The goal of this study was to apply finite-element (FE) analysis to estimate the compressive strain experienced by the LC microcapillary network under IOP and CSFP loads, incorporating a CFD approach to estimate the LC blood flow and oxygen concentration. Specifically, we aimed to study the effects of biomechanical parameters (IOP, CSFP, scleral stiffness, and LC stiffness), morphometric parameters of the LC (pore size, LC radius, cup depth, and microcapillary arrangement), and a hemodynamics parameter (perfusion pressure) on the blood flow and oxygen distribution within the LC. 

Methods

The general overview of our methodology is as follows: First, we generated 5000 eye models, each with varied morphometric parameters (LC radius, cup depth, LC pore size, and LC microcapillary arrangement, referred to herein as anisotropy), biomechanical parameters (LC and scleral stiffness), loading conditions (IOP and CSFP), and a hemodynamics characteristic (perfusion pressure), using custom-made software. The methods used to generate the network and descriptions of the parameters were described in our previous work.¹⁸ So, for each eye model we varied eight unique parameters (excluding the IOP) simultaneously by randomly assigning a value to each parameter based on a predefined normal distribution (described in Table 1). The IOP was varied from 10 mm Hg to 55 mm Hg in increments of 5 mm Hg, and 500 models were built for each IOP value. Second, we performed FE analysis using the FEBio Software Suite¹⁹ to estimate the compressive strain in the LC under IOP and CSFP loads. Third, we calculated the change in dimension (compression) of the microcapillary lumen diameter based on the compressive strain derived from FE analysis. Fourth, we employed CFD analysis to approximate blood flow characteristics in the LC microcapillary and oxygen concentration of the RGC axons situated within the LC. 

Table 1.

 

Study Parameters (Independent Variables), Normal Distribution, and Range of Values

 

General Assumptions of the LC Microcapillary Network Morphology and Blood Flow

We made the following key assumptions in all our models: (1) The LC microcapillary network is a curved surface that follows the shape of the anterior LC surface. (2) Each LC collagenous beam contains a capillary in its center supplying nutrients and oxygen to the surrounding axons.²⁰ (3) The LC beams, and thus the microcapillaries, are parallel to the LC surface. (4) The microcapillary lumen have a uniform diameter of 8 µm prior to being loaded with biomechanical forces.²¹ (5) Oxygenated blood, originating from the Circle of Zinn–Haller or the short posterior arteries, feeds the LC uniformly at its periphery.²² (6) Venous drainage of the LC microcapillaries occurs through the central retinal vein in the LC center.²³ (7) The central retinal vein pressure is equal to the value of IOP,²⁴ and the arterial pressure is equal to the central retinal vein pressure plus the perfusion pressure. (8) RGC axons are homogeneous and are distributed uniformly across the LC. (9) The LC microcapillary network and the oxygen-consuming neural tissues are situated in the middle of the LC thickness (in this study, we assumed that the LC tissue thickness was 300 µm and the LC microcapillary network was situated at a depth of 150 µm with respect to the anterior LC surface). (10) Finally, the microcapillaries lumen retain a circular cross section under compressive load. 

FE Model of the Eye

We adapted a half-eye model (Fig. 1a) with a plane of symmetry along a transverse direction, based on previous works by Wang et al.²⁵ and Sigal et al.,²⁶ incorporating a scleral shell inner radius of 12 mm and an optic nerve length of 2 mm. The thicknesses of three layers of the globe—sclera, choroid and retinal—were 0.9 mm, 0.2 mm, and 0.249, mm respectively.²⁷ The ONH region was comprised of the LC, prelamina tissue, choroid, and peripapillary sclera. The optic nerve region consisted of dura, pia, and post-LC tissue. The thickness of the LC was 300 µm. Constitutive models and their biomechanical properties were assigned for the eye tissues according to Supplementary Material A. 

Figure 1.

 

(a) Finite-element model of the eye. (b) Optic nerve head LC microcapillary network embedded in the LC tissue. Blood flow originates from the LC periphery and drains in the central retinal vein. (c) LC posterior displacement under 50 mm Hg IOP load. (d) Illustration of a single-element LC with a single microcapillary passing through the middle of the element. As IOP and CSFP loads are applied, the element will undergo compression along the axial axis (y-axis). Microcapillaries embedded in the element will proportionally change their shape.

Constitutive models and their biomechanical properties were assigned for the eye tissues according to Table 2. Specifically, the sclera was modeled as a Mooney–Rivlin material with reinforced collagen fibers (following a von Mises distribution for orientation).²⁸ The collagen fibers in the peripapillary sclera that surrounded the peripapillary border tissues (PBTs) were organized into a ring.²⁹ The PBTs, the pia, and the dura mater were assumed to be isotropic elastic (neo-Hookean) with the same stiffness. For these tissues, the material properties were obtained from Wang et al.⁴ The LC was also modeled as an isotropic elastic material (neo-Hookean) for simplicity. Incompressibility and hyperelasticity were assumed for all soft tissues. 

Table 2.

 

Constitutive Model and Biomechanical Properties for Eye Tissues

 

A convergence test, using the LC effective strain as reference, was conducted to derive the appropriate mesh size of the model. Specifically, we increased the number of elements from 8640 in increments of 1000 until the error of the LC average effective langrage strain (compared with the previous number of elements used) was less than 1e-5. Based on the convergence test, an appropriate number of elements for all FE models was determined to be 13,640. All tissues were meshed using eight-node hexagonal elements. Mesh generation and solving of the FE models were accomplished by using FEBio, a nonlinear FE solver designed for biomechanical studies.¹⁹ 

For boundary conditions, the posterior end of the optic nerve was fixed in place, and the nodes that belonged to the plane of symmetry were not allowed to displace out of this plane. There were no contact definitions, as all of the nodes along the tissue interfaces were shared. We applied two loading conditions, the IOP and the CSFP, at the inner limiting membrane and within the arachnoid space, respectively. Across FE models, we varied the following characteristics simultaneously to study their influence on the LC blood flow and oxygen concentration—LC radius, cup depth, LC stiffness, sclera stiffness, and CSFP—based on the normal distributions of the variables (mean ± SD) obtained from other literature (Table 1), whereas IOP was varied in increments. We defined cup depth as the distance from the insertion point of the anterior surface of the LC to the lowest point of the anterior surface (Fig. 1a). Scleral stiffness was varied by multiplying a ratio value (1 ± 0.2) by baseline values of Mooney–Rivlin Von Mises distributed fiber coefficients. 

FE Analysis to Predict Compressive Strain and Microcapillary Lumen Deformation

In order to translate the strain experienced in the LC to the microcapillary network, we assumed that the microcapillary network was situated in the middle of the LC thickness (Fig. 1b) and that any strain and displacement experienced in the LC element (eight-node hexahedron) were translated to the segment of the microcapillary network embedded within the element (Fig. 1c). Because of the IOP and CFSP loads (Fig. 1d), it was observed that the LC most often exhibited large compressive strains along the axial direction. Assuming homogeneity (same biomechanical properties) between the microcapillary and the LC tissue element, we applied LC compressive strain directly to the microcapillary in the middle of the element, resulting in a compression of the microcapillary from a circular cross section to an oval cross-section with a minimal radius along the axial direction (Fig. 1d).³⁰ 

Calculation of the average compressive strain in each LC element was simplified, and it was similar to the calculation of strain based on uniaxial compression.³¹ Given the hexahedral mesh, we took an average strain of the four edges along the axial direction (y-direction in Fig. 1d) to be an approximation of the average compressive strain within the LC element according to Equation 1:  

\begin{equation}{\varepsilon _{avg}} = \frac{{\left| {\frac{{\Delta {L_a}}}{{{L_a}}} + \frac{{\Delta {L_b}}}{{{L_b}}} + \frac{{\Delta {L_c}}}{{{L_c}}} + \frac{{\Delta {L_d}}}{{{L_d}}}} \right|}}{4}\end{equation}

(1)

Here, ∆L is the change in length along each edge (a, b, c, and d) as defined in Figure 1d, and L is the original length of each edge. Assuming that the microcapillary lies in the middle of the element, its diameter after compression can be approximated with the following equation:  

\begin{equation}{D_{compressed}} = {D_{original}} - {{\rm{\varepsilon }}_{avg}}*{D_{original}}\end{equation}

(2)

Here, D_compressed is the diameter of the microcapillary after applying the compressive strain, and D_original is the original diameter of the microcapillary (8 µm). We assumed in this study that the lumen of the blood microcapillary always retained its circular cross-section with a decrease in diameter from a baseline of 8 µm to a circular cross section with an equivalent minimum radius of an oval under compressive strain. 

Hemodynamics Model of Blood Flow and Oxygen Consumption

We modeled blood flow as a non-Newtonian fluid whose behavior can be described using the incompressible Navier–Stokes equations (Equations 3 and 4):³² 

\begin{equation}{\rm{\rho }}\frac{{\partial {{\bf u}}}}{{\partial t}} - {\rm{\ \eta }}{\nabla ^2}{{\bf u}} + {\rm{\rho }}\left( {{{\bf u}} \cdot \ \nabla } \right){{\bf u}} = F - \ \nabla p\end{equation}

(3)

 

\begin{equation}\nabla \ \cdot \ {{\bf u}} = 0\end{equation}

(4)

Here, ρ is the blood density (kg/m³), η (Pa·s) is the blood viscosity, u is the blood velocity field (m/s), p is the blood pressure (Pa), and F is the volume force (N/m³). We assumed a fully developed blood flow profile, a zero volume force, and a steady-state condition. The apparent blood viscosity η in microvessels was modeled as a function of vessel diameter and discharge hematocrit according to previous works by Pries et al.^33,34 In addition, the hematocrit partitioning at vessel bifurcations modeled as a function of flow rates, diameters of the vessels, and the hematocrit of the parent vessel, according to the equations (not listed herein) described by Pries et al.³⁵ 

We describe the diffusion of oxygen in tissues and microcapillaries with Equation 5 (Fick's law of diffusion) and consumption of oxygen in tissue with Equation 6 (Michaelis–Menten relationship):  

\begin{equation}{\rm{D\alpha }}{\nabla ^2}{\rm{C}} = {\rm{M(C)}}\end{equation}

(5)

Here, D is the diffusion coefficient of oxygen of the nervous tissues (m²/s), α is the solubility coefficient of oxygen in nervous tissue, C is the oxygen concentration inside the nervous tissues (mol/m³), and M(C) is the oxygen consumption rate of the nervous tissues (mol/m³/s) which can be approximated by the following equation:  

\begin{equation}{\rm{M(C)}} = {{\rm{M}}_0}{\rm{C}}/\left( {{{\rm{C}}_{{\rm{MM}}}} + {\rm{C}}} \right)\end{equation}

(6)

Here, M₀ is the tissue-dependent maximum oxygen consumption rate, and C_MM is the Michaelis–Menten constant corresponding to the oxygen concentration where the consumption rate is 50% of M₀. The value of M₀ for the RGC tissue has not been measured, and we approximated it to be 0.000225 mol/m³/s; this value was close to the M₀ of a tumor tissue³⁶ and resulted in our average eye model (IOP = 15 mm Hg, and all other parameter values set at mean) being barely above the hypoxic level of 5 mm Hg O₂.³⁷ We assumed M₀ to be uniform throughout the nervous tissues.  

\begin{equation}f\left( {{{\rm{C}}_P}} \right) = Q\left( {4{{\rm{C}}_{Hb}} \times H \times {\rm{S}}{{\rm{O}}_2}\left( {{{\rm{C}}_P}} \right) + {{\rm{C}}_P}} \right)\end{equation}

(7)

 

Oxygen transport along a microcapillary segment can be described by Equation 7, where f is the rate of convective oxygen transport within the microcapillary segments (mol/s), C_P is the oxygen concentration in plasma (mol/m³), Q is the volume flow rate (m³/s), C_Hb is the concentration of hemoglobin-bound oxygen in a fully saturated red blood cell, and H is the hematocrit. We assumed C_Hb to be constant for all capillary segments. SO₂(C) represents the oxygen–hemoglobin saturation function approximated by the Adair equation.³⁸ The blood flow and oxygen transport along the microcapillary segment are described in Figure 2a. Conservation of oxygen mass was governed by the following equation:  

\begin{equation}\frac{{{\rm{df}}\left( {{{\rm{C}}_P}} \right)}}{{{\rm{ds}}}} = - q\,(s)\end{equation}

(8)

Here, s is the distance along the microcapillary length, and q is the rate of oxygen diffusion out of the microcapillary per unit length. The oxygen concentration for C_P at the inlet node (s = 0) is 80 mm Hg O₂.³⁹ 

Figure 2.

 

(a) Blood flow parameters (velocity and pressure) were calculated for each microcapillary branch. Oxygen concentration was also calculated along the length of the capillaries as oxygen diffused out of the vessel to be consumed by the RGCs. Hemoglobin travels with blood flow and acts as the main transporter of oxygen. The hematocrit (Hc) was calculated at each branching point of the microcapillaries, such as when the Hc of the parent branch (Hc_Parent) is redistributed to two daughter branches. (b) Oxygen-consuming RGCs are modeled as point sources (blue dots) embedded throughout the microcapillary network (red lines).

At the microcapillary wall–tissue interface, the following boundary condition for diffusion of oxygen out of the microcapillary was applied:  

\begin{equation}q\,(s) = - D{\rm{\alpha R}}\int_{0}^{{2\pi }}{{\frac{{\partial C}}{{\partial r}}d{\rm{\theta }}}}\end{equation}

(9)

Here, R is the vessel radius and r is the radial distance from the centerline. At the vessel–tissue interface (radius R), we assumed that C_P is equal to the tissue oxygen concentration (C). 

We employed a Green's function method, according to computational work by Secomb et al.,⁴⁰ whereby the oxygen-consuming neural tissue nodes were modeled as point sources distributed throughout the microcapillary network (Fig. 2b). The density of the neural tissue node was fixed at 300 nodes/mm² with respect to the LC surface area. On average, each LC contained 1580 neural tissue nodes representing the RGCs. The boundary condition for the neural tissue region was considered to be an infinite boundary with no tissue nodes found outside the LC microcapillary region. A list of constants used in hemodynamics and the Green's function method of computation are provided in Supplementary Material B. 

Study Parameters

In this study, each model was built with a unique set of 9 parameters (listed in Table 1) that defined the morphology, biomechanics, and hemodynamics of the ONH. IOP was varied stepwise, and all other parameters were varied based on a predefined normal distribution. 

Main Outputs

A total of 5000 models (FE model and LC hemodynamics) were solved and the main outputs were the oxygen concentration in the LC neural tissue nodes at steady state, the compressive strain field in the LC, and the LC displacement field. The main parameter we were investigating in this study was the average oxygen concentration, which we defined as the average value of oxygen concentration across all LC tissue nodes in a model. 

Data Processing and Statistical Analysis

We performed a multivariable linear regression for each parameter in Table 1 with respect to the average value of oxygen concentration of all tissue nodes in each LC network, the compressive strain, and the displacement of the LC. We then ranked the influence of each parameter based on the magnitude of the coefficient (β) obtained from multivariable linear regression analysis. 

Specifically, the multivariable linear regression analysis that we used is described by the following equation:  

\begin{equation}{{\rm{y}}_i} = {\rm{\beta }}{X_{i\ }} + {{\rm{e}}_{i\ }}\end{equation}

(10)

Here, y_i is a vector of values for dependent variables (i.e., average oxygen concentration), X_i is the matrix of our study parameters (e.g., LC radius, pore size, cup depth), β is a vector of regression coefficients for each independent variable, and e_i is the vector of error terms with respect to the regression coefficient. The multi-variable regression was computed in MATLAB with a built-in function mvregress, which employs an ordinary multivariate normal maximum likelihood estimation to deduce β. In addition, prior to regression, the value of each study parameter was normalized using standard z-distribution (to obtain a mean of 0 and standard deviation of 1) to account for a difference in each variable's distribution. 

We also identified the hypoxic neural tissue nodes in each model that had oxygen concentrations of less than 5 mm Hg O₂. The level of oxygen concentration that could lead to RGC hypoxia is unknown, and our definition of 5 mm Hg O₂ is in between the hypoxic levels of brain tissue (10 mm Hg O₂)³⁷ and skeletal muscle tissue (1 mm Hg O₂).⁴¹ We further defined hypoxic LC as having at least one hypoxic tissue node (out of 1580 nodes). 

To compare the values of the dependent variables (e.g., average oxygen concentration) between two groups with equal sample sizes (eyes with 15 mm Hg vs. eyes with 55 mm Hg), we used a standard t-test with 95% confidence intervals (MATLAB ttest2 function; MathWorks, Natick, MA, USA). 

Case Study: The Effects of Study Parameters on the Change in Average Oxygen Concentration when IOP was Raised

We performed a comparative study across 2000 additional models, where we increased the IOP from 18 mm Hg to 55 mm Hg and kept the CSFP constant at 13 mm Hg. The rest of the parameter values followed the distributions given in Table 1. We then performed a multivariable linear regression to rank the influence of each parameter in Table 1, except the IOP and CSFP, with respect to the magnitude of the change in oxygen concentration after the IOP increase. 

Results

LC Displacement and Change in Microcapillary Lumen Diameter

We observed that, for every model, the LC displaced posteriorly (cupping) under IOP and CSFP loads, with an average posterior displacement of 50 µm, ranging from 8 to 133 µm. The LC microcapillary network, which was embedded in the middle of the LC tissue, also displaced proportionally. The average magnitude of compressive strain in the LC for all models was 0.7 ± 0.3%, resulting in a 0.06 ± 0.02 µm average decrease in the microcapillary lumen diameter. According to multivariable linear regression analysis (Table 3c), an increase in IOP was found to be the most important parameter that increased the LC posterior displacement and decreased the microcapillary lumen diameter. 

Table 3.

 

Multivariable Regression Results

 

Differences in Average Oxygen Concentration Between Each LC Region

The average oxygen concentration (mm Hg O₂), with respect to each region of the LC, across all 5000 models is shown in Figure 3a. Specifically, the average oxygen concentration in the peripheral region (40 mm Hg O₂) was significantly higher (t-test P < 0.05) than it was in the central region (30 mm Hg O₂). In addition, the average oxygen concentration at the nasal–temporal region (38 mm Hg O₂) was significantly greater (t-test P < 0.05) than in the superior–inferior region (32 mm Hg O₂), resulting in a distinct hourglass pattern of oxygen distribution as shown in Figure 3b. The overall average oxygen concentration across all LC networks was 38.2 ± 7.1 mm Hg O₂. We also found that 80% of the hypoxic neural tissue nodes across all models were located in the superior–inferior region. 

Figure 3.

 

(a) Average values of oxygen concentration with respect to each region of the LC across all 5000 models. (b) Heat map of the density of hypoxic tissue (less than 5 mm Hg O₂) nodes across all models. This was created by selecting hypoxic tissue nodes for every model and plotting its relative location on a single map. Red areas have a high density of hypoxic tissue, and blue areas have a low density of hypoxic tissue. This distinct hourglass pattern reveals a superior–inferior region containing more hypoxic tissue.

Influence of Study Parameters on the Average Oxygen Concentration

The factors influencing average oxygen concentration the most were, in order of importance, perfusion pressure (β = +4.79), LC radius (β = –3.42), pore size (β = –1.65), anisotropy (β = +0.68), IOP (β = –0.49), LC stiffness (β = +0.30), scleral stiffness (β = +0.21), CSFP (β = –0.11), and cup depth (β = –0.05). The rankings based on regression results are summarized in Figure 4a. 

Figure 4.

 

(a) Bar chart showing the result of multivariable linear regression to rank the importance of all independent variables with respect to average LC oxygen concentration. (b) Scatterplot showing the effects of LC radius on average LC compressive strain and average LC oxygen concentration. (c) Scatterplot showing the effects of scleral stiffness on compressive strain and average oxygen concentration.

Effects of LC Radius on the Average Oxygen Concentration

LC radius was one of the most influential factors that affected both the compressive strain and average oxygen concentration in the LC. In general, an increase in LC radius of 1 SD (+0.1 mm) resulted in a 0.07% increase in compressive strain and a 3.43-mm Hg O₂ decrease in average oxygen concentration. Scatterplots showing the effects of LC radius on compressive strain and average oxygen concentration are provided in Figure 4b. 

Effects of Sclera Stiffness on the Average Oxygen Concentration

Sclera stiffness was ranked second to last among all the parameters according to its influence on oxygen concentration; however, scleral stiffness still had notable effects on both compressive strain and average oxygen concentration. An increase in scleral stiffness of 1 SD (+20% in stiffness coefficients) resulted in a 0.04% decrease in compressive strain and a 0.17-mm Hg O₂ increase in average oxygen concentration. Scatterplots showing the effects of scleral stiffness on compressive strain and average oxygen concentration are given in Figure 4c. 

Effects of IOP on the Average Oxygen Concentration

IOP was the most important parameter in affecting the compressive strain according to the ranking Table 3b. The relationship between compressive strain and IOP is provided in Figure 5a. An increase in IOP from 15 mm Hg to 30 mm Hg resulted in a 0.3% increase in compressive strain. With respect to IOP influence on oxygen concentration, IOP was ranked fifth among the nine parameters (Fig. 4a). 

Figure 5.

 

(a) Box plot showing the relationship between average compressive strain and IOP for 5000 models. Each bar represents 500 models with the same value of IOP. (b) Box plot showing relationship between average oxygen concentration and IOP for 5000 models. (c) A single model comparison of strain in the optic nerve head under IOP of 15 mm Hg and IOP of 55 mm Hg. The values for other parameters are kept at mean. (d) A single model comparison of average oxygen concentration in the LC under IOP of 15 mm Hg and IOP of 55 mm Hg.

Between eyes with IOPs of 15 mm Hg (n = 500) and 55 mm Hg (n = 500), the eyes with IOP of 15 mm Hg had higher (+1.8 mm Hg O₂) average oxygen concentrations than the eyes with IOP of 55 mm Hg (t-test P < 0.05). Eyes with IOP of 55 mm Hg were also more susceptible to having hypoxic tissue, with 8.1% of the high (55 mm Hg) IOP eyes containing at least one hypoxic node in comparison to 3.2% in the low (15 mm Hg) IOP eyes. The adverse effects of IOP increase were less pronounced when we considered the eyes with a more conventional pathological IOP value of 30 mm Hg. These eyes had a slightly lower (–0.7 mm Hg O₂), albeit statistically significant (t-test P < 0.05), average oxygen concentration than eyes with IOP of 15 mm Hg. Figures 5c–5d depict the influence of elevated IOP on compressive strain and average oxygen concentration in a single model. 

Effects of Study Parameters on the Change in Average Oxygen Concentration when IOP was Elevated

When the IOP was raised from 18 mm Hg to 55 mm Hg (2000 cases), the average change in oxygen concentration was –1.75 ± 2.3 mm Hg O₂. The factors influencing the change in average oxygen concentration the most were, in order of importance, cup depth (β = +1.28), LC stiffness (β = –1.05), radius (β = +0.45), perfusion pressure (β = –0.25), anisotropy (β = –0.038), scleral stiffness (β = –0.033), and pore size (β = –0.008). The rankings based on regression results are summarized in Figure 6. 

Figure 6.

 

Bar chart showing the result of multivariable linear regression to rank the importance of independent variables (except IOP and CSFP) with respect to the magnitude of change in oxygen concentration when IOP is increased from 15 mm Hg to 55 mm Hg (other parameters remained constant) across 2000 models.

Discussion

In this study, we developed a novel combined CFD and FE eye model to predict blood flow and oxygen concentration of the LC under biomechanical loads. We were able to use the model to finely control and investigate the most influential morphometric, hemodynamic, and biomechanical parameters. Our models predicted that LC oxygen concentration was highly influenced by the LC radius and perfusion pressure. Biomechanical factors, such as the IOP, CSFP, and scleral stiffness, were ranked lower in their influence as compared with LC morphometric and hemodynamic factors, but they still had an important impact on the LC oxygen concentration. Specifically, higher IOP increased the compressive strain in the LC and resulted in a lower average oxygen concentration. 

Across models, our FE analysis predicted an increase in ONH posterior deformation at larger IOPs, with an average posterior displacement of 50 µm (consistent with reported experimental values).⁴² We also found that compressive strain increased gradually from the peripheral region of the LC toward its center, resulting in a greater reduction of microcapillary radius in the central region. This implied that the microcapillary was more prone to collapse in the central region of the LC under compressive strain; however, the average magnitude of the decrease in blood capillary lumen diameter was relatively small (0.06 µm) compared with the baseline diameter of 8 µm, and we did not observe any total collapse of the capillaries within our model's assumptions. A decrease in blood capillary diameter increased the flow resistance and decreased the blood flow rate (according to our preliminary study, an average collapse of 60 nm resulted in a 5% reduction in flow rate), resulting in less efficient oxygen transport and lower oxygen concentrations in the LC tissue. 

With respect to average oxygen concentrations in each region of the LC (Fig. 3a), the superior and inferior quadrants had notably lower average oxygen concentrations than the nasal and temporal quadrants. This, in effect, created an hourglass pattern of oxygen distribution that could be related to the pattern of neural axon loss reported in previous literature,^43,44 where the neural tissue in the superior and inferior regions of the optic nerve was more susceptible to damage in glaucoma. We also plotted the distribution of hypoxic tissue nodes across all models (Fig. 3b) and found that approximately 80% of hypoxia instances were located in superior and inferior regions. In addition, we found that the inner region of the LC (inner 20% of the diameter) had markedly lower oxygen concentrations than the outer region, consistent with our previous work.¹² We found no significant differences in oxygen concentration between superior versus inferior regions and nasal versus temporal regions. It should be noted that the difference in average oxygen concentration between regions was an expected result given that our morphological definitions created an hourglass pattern of blood vessel density within the LC. Also, it could be inferred without the use of simulations that the inner region of the LC should have lower average oxygen concentrations. However, we have shown here, for the first time, that the differences in oxygen concentration between each region were statistically significant, using a large number of computational models based on the detailed morphological characteristics of the LC. 

Similar to our previous work,¹² we found that the most influential parameters with regard to LC average oxygen concentrations were perfusion pressure and LC radius. One notable difference was that LC pore size is now a more important parameter in our predictions, due to the inclusion of neural tissues in the LC pores in this study, with an increase in average pore size resulting in a considerable decrease in average oxygen concentration. This may account for the observed hourglass pattern of oxygen distribution, as superior and inferior regions of the lamina cribrosa had larger pores than the nasal and temporal regions.⁴⁵ We also found that the LC oxygen concentration was primarily affected by the LC morphometric features (pore size and LC radius) and hemodynamics parameter (perfusion pressure), rather than the biomechanical parameters (IOP, CSFP, and tissue stiffness), according to our ranking. Thus, the LC morphometric features, in particular the radius, could be the main predictors of neural tissue insult due to ischemia. 

We found that the LC radius had a twofold effect on the average oxygen concentration. First, an increase in LC radius significantly increased the LC compressive strain (Fig. 4b), resulting in greater compression of capillaries and worse oxygen transfer. This effect can be explained biomechanically as the increase in lamina cribrosa diameter will increase the distance between a soft LC tissue and a supportive peripapillary sclera ring, resulting in impeded transfer of IOP-induced stresses and strains away from the LC.⁴⁶ Second, the increase in LC radius would increase the total length of the microcapillaries in the network, resulting in an increase in flow resistance (according to Poiseuille's equation) and a decrease in blood flow rate. Therefore, it can be said that the LC radius was important both biomechanically and hemodynamically. A slight increase in LC radius (+0.1 mm from a mean of 1.2 mm) would decrease the average oxygen concentration by 3.5 mm Hg O₂. This may explain why people of African descent, who have large optic discs,⁴⁷ are more susceptible to glaucoma.⁴⁸ In addition, it was observed that the normal tension glaucoma population tends to have larger optic discs,⁴⁹ which could suggest that one of the underlying causes for normal tension glaucoma could be due to poor LC hemodynamics, but this remains to be demonstrated. 

We also found that eyes with stiffer sclera (both peripapillary and peripheral sclera) had lower compressive strain, which resulted in better blood flow and oxygen distribution (Fig. 4c). Despite having a small impact as compared with other factors in this study, scleral stiffness was still an important parameter to investigate given the fact that the sclera has been shown to become stiffer with aging.⁵⁰ It has also been suggested that a stiff sclera provides protective mechanisms for the ONH from IOP-induced deformation and strains.⁵¹ Our result support the idea and further suggest that a stiff sclera could also be beneficial to the LC blood flow. 

The biomechanical effects of IOP on the eye have been extensively studied; however, the impact on ONH blood flow is still not well understood. In this study, we found that an increase in IOP is predicted by our models to significantly increase compressive strain and decrease the average oxygen concentration (Fig. 5a). In particular, looking at the increase of IOP from 15 mm Hg to 55 mm Hg in one eye (Fig. 5d), the average oxygen concentration decreased from 40.3 mm Hg O₂ to 36.8 mm Hg O₂, and the hourglass pattern became more pronounced. Also, we found that eyes with higher IOP (55 mm Hg) were more likely to contain hypoxic tissue. Our study provides a first glimpse of the direct impact of IOP on LC blood flow and oxygen distribution, and the results suggest that an elevated IOP may contribute directly to neural tissue ischemia. 

When the IOP was raised from 18 mm Hg to 55 mm Hg in 2000 eyes, we observed that the cup depth and the LC stiffness had become the most important parameters that influence the magnitude of change in oxygen concentration (Fig. 6). We found that the decline of oxygen concentration (when IOP was raised) was greatly magnified with a greater cup depth and a weaker LC; however, both of these parameters were found to be insignificant in affecting the average oxygen concentration in general (Fig. 4a). The effects of cup depth was particularly interesting because the effect of cup depth on stress and strain of the ONH was considered limited,⁵² and its influence only became prominent when we considered it together with IOP changes. In short, the LC cup depth and LC stiffness would aggravate the effects of the change in IOP. This example highlights the complex interactions that may arise among parameters and how we could utilize the model in a specific context to address clinically relevant questions. 

The IOP, in terms of ranking, had less of an impact on oxygen concentration as compared with other factors. Because the IOP-induced compressive strain in the LC was relatively low (up to 2.5%), it therefore makes sense that microcapillaries would not be affected that much by such deformations. It also means that an elevated pressure range normally observed in glaucoma patients (25–35 mm Hg⁵³) would have minimal effects on the average LC oxygen concentration. Despite the small impact that an elevated IOP had on the oxygen concentration, this effect may become significant over a long progression period of glaucoma disease. 

The following limitations apply to this study. First, we modeled the LC as a surface, similar to our previous study.¹² This is not biologically accurate, as the LC has a thickness of approximately 270 µm,⁵⁴ and recent studies have highlighted a degree of tortuosity of the LC microcapillaries through the LC depth.¹⁶ However, the LC thickness value was small compared with its radius, and our model should provide a good approximation of the general oxygenation state of the LC along its middle cross-section but not in its anterior and posterior regions. Future studies could incorporate a full three-dimensional model of the LC, providing detailed histological data can be obtained. 

Second, the artificial LC microcapillary network used in our work may not be an accurate representation of the physiological network. We constructed our LC microcapillary network based on the reported values of pore size (in terms of area, not volume), the arrangement of the collagen beam, the size of the LC, and curvature. Several important parameters that could impact our results have not been reported in earlier literature, such as the LC microcapillary lumen diameter, density of the microcapillary, and extent of connectivity among vessel branches. We assumed here that the LC microcapillary existed only inside the collagen beam and that the vessel branches were fully connected to each other, but these aspects remained to be investigated physiologically. In addition, as a result of the assumption that the LC had no thickness, the connective tissue volume fraction⁵⁵ was not represented in our study. In order to formulate a more accurate model for LC hemodynamics, we would need additional physiological data for the microcapillary network, particularly quantification of its morphology in three dimensions, and we would need to revisit our assumptions to recreate a model that accounts for the full thickness the LC tissue. 

Third, we modeled microcapillary compression under strain as a change in radius of a circular cross-section. This assumption simplified the solid–fluid interaction of our model and allowed us to perform an approximation of LC blood flow without using a three-dimensional mesh for the capillaries, which would be costly computationally. Hemodynamics in microcapillaries of this size (8 µm) is very complex because of the phase separation phenomenon between red blood cells and the plasma layer. We have accounted for this phase separation effect⁵⁶ to a certain degree by adopting the relationship between the viscosity and tube diameter (circular cross-section) as described by Pries et al.³⁴ However, as illustrated in Figure 1d, the microcapillaries would be more appropriately described by an ellipsoid cross-section³⁰ when they are compressed, rather than a circular one, and our approximation of the properties of blood flow (e.g., viscosity) may not be accurate. Unfortunately, there was no experimental information available for microcapillary hemodynamics in a complex lumen geometry. It was likely that we underestimated the effects of reduction in blood flow velocity when it was compressed, as irregularly shaped microcapillaries can further disrupt the lubricating plasma cell-free layer.⁵⁷ 

Fourth, we assumed steady-state approximations in our models, and the cyclic nature of IOP⁵⁸ and CSFP⁵⁹ were neglected. The IOP can vary at a rate of 2.7 mm Hg/s⁵⁸ due to ocular pulse. Also, the CSFP was shown to be pulsatile and out of sync with IOP fluctuations⁵⁹; thus, the translaminar pressure that is the main determinant of the compressive strain is dynamic. A transient model of both biomechanics and hemodynamics should be developed in the future to study the effects of dynamic loading on LC blood flow. 

Fifth, we modeled a one-way interaction of biomechanical forces and hemodynamics; thus, we did not account for the forces exerted by blood flow on the ocular tissues. Considering a macroscale model of the ONH (a homogeneous LC with uniform biomechanical properties), a change in microcapillary blood pressure in the LC tissue may alter its stiffness (e.g., the LC will become stiffer with an increase in microcapillary pressure) and the magnitude of microcapillary compression.⁶⁰ An eye model with solid–fluid interaction that includes a feedback loop will be the target for our next work, and it will be especially useful when we want to look at dynamic models accounting for variations in IOP and blood pressure. 

Sixth, autoregulation, a parameter that can directly affect the radius of the microcapillaries, was not included in our study. It has been suggested that autoregulation impairment at the ONH could be a contributing factor in glaucoma development.⁶¹ Our model could readily accommodate any parameters that affect the radius through the addition of further governing equations that fit the experimental data. However, autoregulation in the LC has not been well characterized, due to the limited penetration depth of conventional techniques to measure blood flow (e.g., laser speckle flowgraphy and laser Doppler flowmetry), and current measurements have been limited to the pre-laminar region only.^62,63 We believe that more information is required before we can come up with a biologically accurate model of LC autoregulation. Notwithstanding these challenges, inclusion of autoregulation in the model remains our primary goal for future improvements. 

Seventh, we did not model the geometry of blood vessels surrounding the LC, such as the Circle of Zinn-Haller (CZH) and the central retinal vein (CRV). The IOP elevation could compress the CRV,⁶⁴ reducing the perfusion pressure; thus, the impact of IOP on the LC hemodynamics could be more dramatic. The same reasoning applies for the CZH and relevant arteries⁶⁵ surrounding the ONH. The logical next step would be to include the CRV geometry in our model, which has been attempted in a computational study.⁶⁶ Because the CZH structure was more elusive, we would need more information about its morphology before we could include it in a computational model. 

Eighth, we assumed a linear relationship between each study parameter in Table 1 and oxygen concentration. We also used multivariable linear regressions to determine the importance of each parameter; however, the relationships may not be linear, and some relationships may also exhibit weak linear correlations. To accurately predict the relationships between each study parameter and the average oxygen concentration, each parameter would have to be isolated and studied in full detail with many more simulations (which is beyond the scope of the present work). Also, one should be cautious in interpreting the rankings of the less important variables—particularly scleral stiffness, LC stiffness, CSFP, and cup depth—as the magnitude of the beta coefficients were similar and the rankings may not be definitive. 

Ninth, the maximum oxygen consumption rate (M₀) of the RGC tissue and the oxygen concentration that can cause RGC tissue hypoxia have never been measured. We have to approximate the values for both of these parameters, which introduces uncertainties to our conclusion about hypoxia and oxygen concentration. We believe that it is very crucial to characterize the oxygen consumption rate of the RGC in order to fully understand the clinical implications of the variations in ONH blood flow. 

Finally, LC blood flow and oxygen concentration have still not been characterized accurately with current imaging techniques, which make it difficult to validate our models. Nevertheless, recently there have been promising efforts to image blood flow volume in the microcapillaries,⁶⁷ and we expect relevant experimental data to be available in the near future. 

In summary, we have implemented a comprehensive computational framework to study the effects of IOP and a wide range of other parameters on the hemodynamics of the LC. Our results complement clinical observations, such as the presence of an hourglass pattern for the oxygen concentration and the adverse effects of IOP, that could potentially be useful to further our understanding of glaucoma. In the future, our approach could readily be improved with the addition of experimental data and has the potential to ultimately become a numerical testbed for the study of glaucoma. 

Acknowledgments

Supported by the Ministry of Education, Academic Research Fund Tier 2 Grants (R-397-000-280-112, R-397-000-308-112[MJAG]). 

Disclosure: T. Chuangsuwanich, None; P.T. Hung, None; X. Wang, None; L.H. Liang, None; L. Schmetterer, None; C. Boote, None; M.J.A. Girard, None