Factors Influencing Lamina Cribrosa Microcapillary Hemodynamics and Oxygen Concentrations

Abstract

Purpose: To identify and rank the lamina cribrosa (LC) morphologic factors that influence LC microcapillary hemodynamics and oxygen concentrations using computational fluid dynamics (CFD).

Methods: We generated 12,000 ‘artificial' LC microcapillary networks and predicted blood flow velocities and oxygen concentrations within the microcapillaries using CFD. Across models, we varied the average pore size of the LC (5500 ± 2400 μm²), the microcapillary arrangement (radial, isotropic, or circumferential), the LC diameter (1.9 ± 0.3 mm), the inferior–superior curvature (340 ± 116 m⁻¹), and the nasal–temporal curvature (−78 ± 130 m⁻¹). We assumed that blood flow originated from the Circle of Zinn-Haller, fed the LC uniformly at its periphery, and was drained into the central retinal vein. Arterial (50 ± 6 mm Hg) and venous (17.7 ± 6 mm Hg) pressures were applied as boundary conditions and were also varied within our simulations. Finally, we performed linear regression analysis to rank the influence of each factor on LC hemodynamics and oxygen concentrations.

Results: The factors influencing LC hemodynamics and oxygen concentrations the most were: LC diameter, arterial pressure, and venous pressure, and to a lesser extent: the microcapillary arrangement (anisotropy) and nasal–temporal curvature. Lamina cribrosa pore size and superior–inferior curvature had almost no impact. Specifically, we found that LCs with a smaller diameter, a radial arrangement of the microcapillaries, an elevated arterial pressure and a decreased venous pressure had higher oxygen concentrations across their networks.

Conclusion: This study described LC hemodynamics using a computational modeling approach. Our study may provide clinically relevant information for the understanding of ischemia-induced neuronal cell death in optic neuropathies.

Open angle glaucoma (OAG) is an optic neuropathy disease, characterized by a progressive loss of retinal ganglion cells (RGC), and permanent vision loss. A well-known risk factor for OAG is an elevated IOP and it is currently considered the only treatable factor. However, the link between the pathogenesis of OAG and elevated IOP is not firmly established.¹ Up to one-third of glaucoma patients do not have elevated IOP² and a significant number of people with elevated IOP never develop glaucoma.³ 

Intraocular pressure is known to deform the posterior structures of the eyes, particularly the structurally weak lamina cribrosa (LC)^4,5—a major site of RGC damage in glaucoma. The LC is a sieve-like network of collagenous connective tissue beams located in the optic nerve head (ONH) containing a comprehensive network of capillaries that provides nutritional and oxygen support to the RGCs.⁶ The LC undergoes morphologic changes (e.g., ‘cupping') and significant remodeling during OAG development^7–10 and these changes have become the focus of many studies that have attempted to link changes in IOP to RGC death. A prominent hypothesis suggests that IOP-induced ONH strains damage the RGCs either directly¹¹ or indirectly through a mechanical disruption of axonal transport in RGCs, thus depriving RGCs from essential trophic factors.^12–14 An alternative hypothesis suggests that ONH strains can alter hemodynamics within the LC, which will in turn reduce the diffusion of nutrients to astrocytes and/or induce ischemia, thus resulting in RGC death. 

To date, there is evidence to support the latter vascular hypothesis. Several studies have reported a correlation between poor ocular hemodynamics and glaucoma^15–17; for instance, Kanakamedal et al.¹⁸ have suggested that people of African descent had a much higher prevalence of OAG compared with those of European descent due to their relatively weaker vascular components. The study also revealed a significant correlation between changes in blood flow and ONH morphology. Unfortunately, accurate measurement of blood circulation in vivo, especially at the level of the LC, has remained extremely challenging. This represents a barrier to establish impaired ocular hemodynamics as a contributing cause of OAG. 

In view of the technical limitations, computational fluid dynamics (CFD) is a powerful alternative to understand how deep ocular blood flow could be involved in the development and progression of OAG. A CFD study by Carichino et al.¹⁹ revealed that IOP-induced LC deformation could reduce the lumen size of the central retinal artery (CRA), and thus the overall ocular blood flow volume. A computational study by Causin et al.²⁰ also suggested that an increase in IOP could disrupt LC hemodynamics through biomechanical interactions. However, no CFD studies have yet modeled the LC as a detailed capillary network and studied how changes in ONH morphology could affect hemodynamics and oxygen concentrations within such a capillary network. We believe it is of critical importance to understand how ONH morphologic changes impact the LC blood flow at the microscale, because this could be a major reason for RGC death in OAG pathophysiology. 

The goal of this study was to provide a modeling study of LC hemodynamics and oxygen concentration at the microscale. Specifically, we aimed to identify the main LC morphologic parameters that could affect LC blood flow and oxygen distribution. Because a majority of these morphologic parameters can now be measured in vivo with optical coherence tomography (OCT), our study may be relevant for the clinical management of glaucoma. 

Methods

Our General Approach to Study LC Hemodynamics

In this study, our aim was to identify and rank the morphologic factors influencing the hemodynamics and oxygen concentrations in the microcapillaries of the LC. To this end, we generated 12,000 artificially-reconstructed LC microcapillary networks with varied morphologic characteristics. A CFD approach was then used to predict blood flow velocities and oxygen concentrations in all LC networks. 

Assumptions on LC Morphology and Microcapillary Blood Flow

For all our models, several general assumptions were made: (1) The LC is a curved surface with principal curvatures aligned along the inferior–superior and nasal–temporal directions,²¹ (2) each LC collagenous beam contains a microcapillary in its center that supplies nutrients and oxygen to adjacent axons,²² (3) the LC beams, and thus the microcapillaries, are parallel to the LC surface, (4) the LC microcapillaries have a uniform diameter of 8 μm,²³ (5) oxygenated blood, originating from the Circle of Zinn-Haller or the short posterior arteries, feeds the LC uniformly at its periphery,²⁴ and (6) venous drainage of the LC microcapillaries occurs through the central retinal vein that passes through the center of the LC.²⁵ 

Generation of Artificial LC Microcapillary Networks

Under our general assumptions, we used a custom-written algorithm (Matlab 2014b; The Mathworks, Inc., Natick, MA, USA) to generate artificial LC microcapillary networks that were representative of a healthy (i.e., nonglaucomatous) human LC. Across networks, we varied the pore size of the LC, the microcapillary arrangement, the LC diameter, and the main LC curvatures (Table 1). 

Table 1

 

Input Parameters Used for the CFD Simulations

 

Table 2

 

Table of Constants Used in the Governing Equations

 

LC Network Generator Algorithm.

Given a set of control parameters (defined in subsequent sections), we numerically generated multiple ellipses in an enclosed two-dimensional (2D) space (Fig. 1a). Each new ellipse was discarded if it overlapped with its neighbors, but kept if it did not. We iterated the process until the space was densely packed and additional ellipses could not be fitted into the space. 

Figure 1

 

(a) An initial image with large circles and radially aligned ellipses. (b) A LC capillary network generated from image (a). (c) An initial image with small circles and radially aligned ellipses (in the peripheral region). Diagonal dashed-lines define the N-T and S-I regions. The dashed-circles section the central and peripheral regions of the LC. The outermost solid circle represents the LC diameter and is used to crop the LC network. (d) A LC capillary network generated from (c). The locations of arterial inlets (red circles) and venous outlets (blue) are highlighted. The Letter A corresponds to the central retinal artery and the letter V to the central retinal vein.

The negative ‘ellipse' image (white region in Fig. 1a) was then ‘skeletonized', a process whereby a thin line (1-pixel thick) was generated within the white region, equidistant to the region's boundaries. A ‘pruning' algorithm was then used to remove isolated branches to ensure that all LC capillaries were connected except at the periphery where blood flowed in. The central retinal artery and vein (off-centered by 0.05 mm toward the nasal direction) were added to complete each LC capillary network.²⁶ The size of the artery (0.02 mm²) and that of the vein (0.03 mm²) were fixed for all networks²⁶ independently of the total LC area. 

For each LC network, our algorithm received four input parameters (LC pore size, anisotropy, diameter, superior–inferior curvature, and nasal–temporal curvature), each of which was a normally distributed random number. 

The final LC capillary image was then converted to a vector file (DXF format) using ImageJ²⁷ and Matlab in order to generate nodal and connectivity information (i.e., mesh) of each LC capillary network, as required for CFD processing. 

LC Pore Size.

Generally, the LC exhibits larger pores and less connective tissue in the superior and inferior regions than its temporal and nasal regions^28–30 (Table 1). Furthermore, the pores are smaller and more uniform in the central area.⁷ To mimic such heterogeneity, all our networks were subdivided into nine regions (Fig. 1c): the central region (inner 20% of the LC diameter), four regions (superior [S], inferior [I], nasal [N], and temporal [T],) in the middle area (defined between 20% and 60% of the LC diameter), and four additional regions (S, I, N and T) in the peripheral area (defined between 60% and 100% of the LC diameter). 

For each LC network, a pore size ‘seed' value δ was selected randomly from a normal distribution with a mean of 0 and a SD of 5.4 μm. From such a value, the minor axis of each ellipse in the central region was randomly chosen between δ + 30 μm and δ + 32 μm; in the N-T peripheral regions between δ + 45 μm and δ + 50 μm; and in the S-I peripheral regions between δ + 58 μm and δ + 81 μm. This method allowed us to generate heterogeneously distributed LC pores with a single control parameter. For instance, Figure 1b was generated with a δ value of 8 μm and Figure 1d with a value of 0 mm. 

LC Anisotropy.

In its periphery, the LC was either considered: (1) anisotropic with a radial arrangement of the LC beams (and thus capillaries) as observed in humans and primates,³¹ (2) isotropic (also observed in humans),³¹ or (3) anisotropic with a circumferential arrangement (hypothetical scenario). For our LC network generator algorithm, the anisotropy input parameter (σ) was defined as the ellipse major-to-minor axis ratio minus 1. For each network, such a parameter was chosen randomly from a normal distribution (mean: 0; SD: 0.05). Here, σ less than 0 indicates anisotropy with a circumferential arrangement, σ equal to zero indicates isotropy, and σ greater than 0 indicates anisotropy with a radial arrangement. Note that a higher absolute value corresponds to more elongated ellipses, and thus a higher anisotropy degree. The central and middle regions of the LC were considered isotropic for all LC networks. 

LC Diameter.

For each network, the LC diameter was used as an input to define the height and the width of the original 2D space (before ellipses packing). This method ensured that LC pore size was independent of the LC diameter. Once ellipses were generated, the LC diameter was further used to crop the LC to a circular shape. We used the following normal distribution of healthy LCs' diameters: 1.9 ± 0.3 mm (Table 1). 

LC Curvature.

For each LC network, we varied its N-T and S-I curvatures K1 and K2 by ‘morphing' the initially flat geometry onto a hyperbolic three-dimensional (3D) plane (of principal curvatures K1 and K2). K1 and K2 followed normal distributions, K1: 340 ± 116 m⁻¹ and K2: −78 ± 130 m⁻¹ (Table 1), as recently measured by our group in a healthy population.²¹ From K1 and K2, and for each LC network, we computed the corresponding global shape index (GSI) – a quantity describing the global shape of the LC and that varied between −1 (cup shape) and +1 (cap shape) (Fig. 3c).²¹ 

Hemodynamics Boundary Conditions

We applied arterial inlet pressure (50 ± 6 mm Hg) and venous (17.7 ± 6 mm Hg)³² outlet pressure as boundary conditions to each LC network and each value was randomly selected from normal distributions (Table 1). Because there was no reported data on the arterial inlet pressure of the LC microcapillary network, we assumed the mean arterial inlet pressure to be 50 mm Hg, which would give an average blood flow velocity of approximately 1.0 mm/s³³—a reasonable assumption for microcapillaries. We also assumed arterial pressure SD to be the same as that of the venous pressure. 

Governing Equations for LC Blood Flow Simulation

Blood Flow Equations.

As a first approximation, we assumed that blood was an incompressible Newtonian fluid,³⁴ whose flow behavior can be described with the Hagen-Poiseuille equation:  

 

Here, η (Pa·s) the blood viscosity, Q the volumetric flow rate (m³/s), Δp is the pressure drop along the vessel (Pa) and L is the vessel length (m) and R_c is the capillary radius.  

Oxygen mass transport inside the microcapillary was solved with the following equations:  

 

where f is the rate of convective oxygen transport within the microcapillary segments (mol/s), C is 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 H and C_HB to be constant for all capillaries segment. SO₂(C) represents the oxygen-hemoglobin saturation function approximated by the Adair equation³⁵:  

 

where α is the oxygen solubility coefficient in blood (mol/m³/mm Hg), α₁ = 0.01524, α₂ = 2.7 × 10⁻⁶, α₃ = 0, and α₄ = 2.7 × 10⁻⁶.³⁵  

Oxygen Diffusion and Consumption Equations.

To predict oxygen diffusion and consumption inside the nervous tissues, we followed the standard Krogh cylinder-type model³⁶ (Fig. 2), in which a tissue cylinder is assumed to receive oxygen supply from a single capillary that lies in its center. The following governing equation (Fick's law of diffusion³⁶) was imposed to the nervous tissues:  

 

where D is the diffusion coefficient of oxygen of the nervous tissues (m²/s), C is the oxygen concentration inside the nervous tissues (mol/m³), M(C) is the oxygen consumption rate of the nervous tissues (mol/m³/s). We assumed that the nervous tissues were homogeneous with the same oxygen diffusion coefficient and consumption regardless of location. M(C) can be approximated by Michaelis-Menten kinetics—the most common model used to model tissue oxygen consumption rate³⁶:  

 

where M₀ is 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₀. We assumed M₀ to be uniform throughout the nervous tissues.  

Figure 2

 

Geometry of the Krogh cylinder-type model. The inner cylinder represents a capillary and the outer cylinder represents the nervous tissue. R_t, nervous tissue radius; R_c, capillary radius; Z, distance along axial direction. A no-flux boundary condition was imposed on the outer nervous tissue boundary (R = R_t). The oxygen flux from the capillary to the nervous tissue occurs at the capillary wall (R = R_c).

Figure 3

 

Quantification of pore size, anisotropy, and curvature. (a) A histogram of pore size distribution in different regions of a standard LC. (b) A vector plot (based on a Fourier Transform method) shows the major alignment of the network in each region of interest. (c) Global Shape Index distribution of the generated LC geometries, showing the morphologic shapes (exaggerated) that correspond to each GSI value.

We further assumed that the transfer of oxygen from the LC microcapillaries to the nervous tissues occurred at the capillary wall (R_c) as governed by the following equation³⁷:  

 

where r is the radial distance (m) from the vessel centerline, R_c is the capillary radius (m), and q is the oxygen flux across the vessel wall (mol/m/s). The term   represents the oxygen flux at the capillary wall, which is a function of the average oxygen concentration inside the capillary.  

In our study, we simplified the 3D Krogh model to a one-dimensional (1D) line model to reduce the computational time to process a large number of complex vascular networks. Equation 1 was first used to solve for blood flow rate inside all capillary segments and Equation 2 was used to approximate plasma oxygen concentration in the segments. The oxygen consumption–diffusion equation (Eq. 4) was applied to the nervous tissue, and by conservation of mass (Eq. 6), we determined the flux of oxygen from the capillaries to the nervous tissues and quantified the oxygen concentration drop along the axial length of each capillary (refer to Table 2 for list of contants used in the equations). Only the average oxygen concentration and average flow velocity were computed in the capillary. These values were representative of the oxygen concentration in the surrounding nervous tissue (Supplementary Material). 

Data Processing and Statistical Analysis

A total of 12,000 simulations were performed in this study. The simulation outputs consisted of three parameters—the oxygen concentration, the blood pressure, and the blood flow velocity inside the microcapillaries. We also identified capillary segments of each LC that had oxygen concentrations of less than 10 mm Hg O₂ as hypoxic regions.³⁸ Because no study has been conducted on the metabolism of LC's nervous tissue, this value was assumed from the study of oxygen-intensive human tumors and it should serve as an upper threshold for normal cells. We performed multiple linear regression on a normalized mean value of oxygen concentration and flow velocity for each LC network, and used the regression coefficient (β) to rank the importance of each LC characteristic (Table 1) in the LC hemodynamics. 

Verification of Generated LC Networks

Each generated LC network was re-evaluated using a digital image processing algorithm to obtain each individual pore size in μm² and its regional distribution. A Fourier transform algorithm³⁹ was used to compute the principal alignment and the degree of anisotropy of the microcapillaries in each local region. 

Results

Verification of Generated LC Networks

We analyzed the various study parameters of our LC networks and found the following. First, the pore size distribution (Fig. 3a) of our models closely matched that of normal human eyes, with S-I regions having significantly larger pores than N-T regions (pore size by region: Nasal: 2200 ± 1061 μm², Temporal: 2351 ± 950 μm², Superior: 5093 ± 1750 μm², Inferior: 5100 ± 1765 μm²). Next, for anisotropy, our LC networks had normally distributed LC anisotropy values that described a range of anisotropic scenarios from radial to tangential capillary arrangements (Fig. 3b). Finally, the distribution of the Global Shape Index (GSI) values (a measure of LC curvature) of our models had an average value of −0.38, corresponding to a saddle-rut shape and similar to that reported in human eyes (Fig. 3c). 

Regional Variations in Oxygen Concentration and Flow Velocity

The average flow velocity (m/s) and average oxygen concentration (mm Hg O₂) across the LC network are shown in Figure 4d. Notably, the average oxygen concentration decreased significantly (P < 0.001) from the peripheral region (64.8 mm Hg O₂) toward the mid (35.8 mm Hg O₂) and central regions (25.5 mm Hg O₂). Variations in oxygen concentration between N-T-S-I regions were not significant. The overall oxygen concentration across all LC networks was 50.4 ± 9.5 mm Hg O₂. 

Figure 4

 

Input geometry and output parameters from the CFD simulations. (a) An example of a LC geometry that is used for CFD simulations. (b) Color-map (top-view) of flow velocity (m/s) within the microcapillary network. (c) Color map (top-view) oxygen concentration (mm Hg O₂) within the microcapillary network. (d) Mean values of oxygen concentration and flow velocity (from all 12,000 simulations) for each region.

Flow velocity increased significantly (P < 0.001) from the peripheral region (0.7 × 10⁻³ m×/s) toward the mid (1.2 × 10⁻³ m/s) and central region (2.9 × 10⁻³ m/s). No significant variations were observed between N-T-S-I regions and the LC networks had an average flow velocity of 1.2 × 10⁻³ ± 0.3 × 10⁻³ m/s. 

Influence of the LC Parameters on Oxygen Concentration, Flow Velocity, and Hypoxic Area

The association between LC parameters and the resulting LC blood oxygen concentration and flow velocity are plotted in Figures 5 and 6, respectively. In addition, we also quantified a third characteristic—hypoxic regions (<10 mm Hg in oxygen) as a proportion of total LC microcapillary area—and looked at its association with all LC parameters ( Supplementary Figure S1). 

Figure 5

 

Scatter density plot of oxygen concentration (mm Hg O₂) with respect to each of the study parameter. Each point represents the average oxygen concentration of a given network (out of 12,000).

Figure 6

 

Scatter density plot of flow velocity (m/s) with respect to each of the study parameter. Each point represents the average flow velocity of a given network (out of 12,000).

We ranked the influence of each LC parameter on the three characteristics based on the magnitude of the coefficient (β) obtained from multiple linear regression analysis (a positive coefficient suggests a positive correlation between the LC parameter and the hemodynamics characteristics and vice versa). First, the factors influencing average oxygen concentration the most were, in order of importance, LC diameter (β = +0.72), arterial pressure (β = +0.58), venous pressure (β = −0.37), anisotropy (β = +0.14), K1 (β = −0.07), pore size (β = −0.01), and K2 (β = +0.002). Second, for average flow velocity, the factors were ranked as follows, arterial pressure (β = +0.69), diameter (β = −0.57), venous pressure (β = −0.44), pore size (β = +0.82), K1 (β = −0.06), anisotropy (β = +0.047), and K2 (β = +0.006). Finally, the ranking of factors associated with the proportion of hypoxic areas was similar to that of average oxygen distribution. The rankings based on regression results are summarized in Figure 7. 

Figure 7

 

Multiple linear regressions coefficient plot of the study parameters with respect to (a) average oxygen concentration (b) average blood flow velocity, and (c) percentage of hypoxic area (<10 mm Hg O₂). Green bar represents positive correlation and red bar represents negative correlation. The error bar is not presented here due to a very small error (<0.1%) for each β value.

Discussion

In this study, we developed a novel method to generate LC microcapillary networks that were incorporated into a CFD model for blood flow analysis. We were able to finely control the important LC morphologic parameters (pore size, anisotropy, diameter, and curvature) and study their effects on LC hemodynamics and oxygen concentration. Our models predicted that LC hemodynamics was highly influenced by the LC diameter, the arterial pressure, and the venous pressure. Notably, LCs with large diameters (>2.0 mm) had significantly poorer oxygen concentrations. Other parameters, such as pore size, anisotropy, and curvature had much less influence on LC hemodynamics. 

We found that the regional variations in both the flow velocity and oxygen concentration were significant along the radial direction, but not in the circumferential N-T-S-I direction (Fig. 4d). On average, oxygen concentrations decreased significantly from 64 mm Hg O₂ in the peripheral region to 25 mm Hg O₂ in the central region. While there is no literature to date that measured human ONH oxygen concentrations in vivo, our range of values is consistent with that observed in mice and porcine models.^40,41 The absolute velocity values obtained from our simulations also expand our knowledge of ocular blood flow in the human LC, because conventional techniques (e.g., laser speckle flowgraphy) yield only relative velocities and cannot provide deep flow measurements within the LC.⁴² 

Lamina cribosa diameter was found to be the most influential parameter affecting both oxygen concentration and hemodynamics. Lamina cribosas with large diameters were found to exhibit a lower average oxygen concentration, a lower average flow velocity and a higher percentage of hypoxic area (up to 80% for LCs larger than 2.0 mm). Previous studies have reported an increase in optic disc diameter in glaucoma patients.^43,44 Combined with our findings, this suggests a plausible mechanism for RGC axonal death in glaucoma: elevated IOP could result in an increase in LC diameter (acutely and chronically),^45–47 which may in turn result in lower blood flow velocities and oxygen concentrations within the LC, starving RGC axons of oxygen, and possibly other nutrients. Further research needs to be done on the role of optic disc enlargement as a triggering cause for glaucoma. Our finding may also explain why certain races have a greater predisposition to OAG.¹⁸ For instance, people of African descent have been shown to exhibit significantly larger optic discs⁴⁸ and poorer vascular components compared with those of European descent.¹⁸ 

Not surprisingly, arterial and venous blood pressures were the other two most significant parameters that affected oxygen concentration and flow velocity. A decrease in arterial pressure and an increase in venous pressure both resulted in lower average oxygen concentration and flow velocity. This result is reasonable, as the aforementioned changes in pressure will reduce the overall perfusion pressure of the LC. It has also been shown that both an increase in central venous pressure and a decrease in arterial pressure are risk factors for the development of OAG,^49–52 and our result suggests that changes in blood pressure can significantly impact LC hemodynamics and potentially contribute to OAG progression. 

We found LC pore size to have a significantly smaller influence on LC hemodynamics when compared with other factors. Changes in pore size had almost no effect on the oxygen concentration, while larger pores slightly increased the average blood flow velocity. Several monkey model studies have shown that there are no significant differences in pore size between normal and glaucomatous eyes.^53,54 A histologic study by Roberts et al.⁸ found an increase in connective tissue beams in experimental glaucoma eyes, but the connective tissue volume fraction (analogous to pore area) remained relatively constant. Our findings further suggest that changes in collagenous beams and the LC microcapillary pore size may not significantly impact LC hemodynamics. 

Anisotropy of the microcapillary network was shown to have greater impact on average oxygen concentration than pore size. Lamina cribosas with radially oriented capillaries at the peripheral region had higher average oxygen concentration than those with isotropic or circumferential configurations. Note that human LCs naturally exhibit a radial arrangement of the microcapillaries at the periphery³¹ and any deviations from such an arrangement (possible through focal defects⁵⁵) may contribute to the development of OAG. 

We found that an increase in curvature along the N-T direction (K1) resulted in both lower average oxygen concentration and blood flow velocity, whereas changes in S-I curvature (K2) had almost no impact on the two output parameters. The increase in N-T curvature corresponds to an increase in LC ‘cupping' and LC depth,²¹ which are well-known to be associated with glaucoma development.^56,57 It was notable in our study that K1 was significantly more influential than K2. This difference may arise from the heterogeneity in LC capillary density. The N-T region has a greater density of microcapillaries than the S-I region (due to smaller pore area). This may mean that changes to curvature along the N-T direction have a greater influence on microcapillary morphology and hemodynamics. 

It should be mentioned that a computational study on LC hemodynamics was previously conducted by Causin et al.²⁰ The authors used a poroelastic model, whereby blood capillaries were modeled as homogeneous pores in an elastic matrix. The work accounted for the biomechanical action of IOP, retrolaminar tissue pressure, and scleral tension on hemodynamics, but it did not take into account complex capillary networks and oxygen transport. On the other hand, our work took those latter into account, but ignored the effects of mechanical stress. We believe both works should be regarded as complementary and contribute to the current understanding of LC hemodynamics. 

The followings are the main limitations of our study. First, we modelled the LC as a surface. This is not biologically accurate as the LC has a thickness of approximately 270 μm.⁵⁸ Considering that the LC may exhibit highly heterogeneous stresses (induced by IOP, the retrolaminar pressure, or during eye movements⁵⁹), it would be logical to observe variations in LC hemodynamics and oxygen concentrations through the thickness of the LC, as predicted by Causin et al.²⁰ LC thickness has also been shown to change with glaucoma progression, which may in turn compromise blood flow further in glaucoma subjects.⁸ While using ‘surface-like' LCs (as performed herein) should provide a first degree of understanding of LC hemodynamics, future improvements to our models should prioritize on incorporating LC thickness and interactions with LC stresses. 

Another limitation is that we simplified the study by reducing our 3D models into 1D ‘line' models in which only the average blood flow velocity and oxygen concentration were solved for each cross section of the capillary. According to our 3D preliminary work (Supplementary Material), the average oxygen concentration inside the capillary segment was found representative of the surrounding neural tissue's oxygenation state, because the oxygen concentration dropped at most by 10% from the capillary to the adjacent LC nervous tissue. However, this assumption undermined the effect of very large pores in our findings as there could be a significant drop in oxygen concentration in the middle of the pore that was not taken into account in our models. Further models, which include nervous tissue between the LC pores, should be built, to better understand the effect of pore size. 

We also imposed several constraints on the CFD simulation. The constraints associated with blood flow equations and oxygen diffusion in tissue are mostly justified based on similar computational studies conducted on capillary blood flow.^35,36 However, the blood viscosity, especially in microcirculation, may not be constant as we assumed. Blood becomes more particulate in small capillaries and its bulk viscosity can be a function of hematocrit, flow velocity, and capillary diameter.⁶⁰ One could try to account for this particulate nature by increasing the bulk viscosity of blood by 30% as proposed by several studies.^61,62 Interestingly, we found that increasing the viscosity did not considerably change the linear regression coefficients of all parameters (coefficients values changed by 5% at most when the viscosity was increased by 30% to 5.2 × 10⁻³ Pa/s) and the ranking of factors remained the same. Changes in average flow velocity and oxygen concentration were observed but remained relatively small (20% and 12%, respectively). However, increasing viscosity may not necessarily make our models more biofidelic as blood viscosity in the microcapillary is determined by many factors, including the hematocrit. Furthermore, variations in viscosity could also happen locally in each capillary segment. 

In addition, our assumption of constant capillary diameter may not be biologically accurate. A computational study by Causin et al.²⁰ also suggested that an increases in IOP could induce changes in the capillaries' diameter, which could compromise the lamina's perfusion. Autoregulation process can dynamically change the LC capillary diameter and autoregulation impairment at the ONH has been considered as a possible cause for glaucoma.^63,64 Unfortunately, no studies have yet provided a comprehensive map of capillary size within the LC. It is plausible that capillary radius could change regionally, and this may be further influenced by autoregulation, and IOP-induced LC stresses (these latter being highly heterogeneous). Note that all our other parameters were derived from experimental and clinical observations. Due to a lack of information, we have considered not to vary capillary radius in our models. Further studies are necessary to quantify the size of capillaries in the human LC. Such information could considerably enhance our models. 

Our findings suggest that the central region of the LC may be more susceptible to ischemia than other regions. This is in contrast with other studies, which indicate that glaucomatous axonal loss starts from the peripheral region.⁶⁵ Also, an hour-glass pattern of axonal loss, in which S-I nerve fibers are selectively more damaged than the N-T nerve fibers,⁶⁵ was not evident from this study. These discrepancies may arise from the inhomogeneity of both the microcapillary characteristics (diameter, autoregulation capabilities, breakage, etc.) and nerve fiber density in each region of the LC. It is also possible that mechanical insult plays a greater role than hemodynamics in creating these damage patterns. Further studies that include the regional variations of both nerve fibers and microcapillaries may reconcile the observed axonal loss pattern and our hemodynamics result. 

Lastly, accurate in vivo measurement of oxygen concentration and blood flow velocity are not yet possible with current techniques. This makes the validation of our model difficult. Nevertheless, our output lies within normal biological values in capillary networks of comparable size. 

In summary, our study revealed, for the first time, the relationship between clinically relevant LC parameters and the hemodynamics and oxygen distribution of the LC. The clear ranking of each parameter's importance presented in the study has the potential to be clinically useful in glaucoma management, diagnosis, and risk profiling. Finally, our approach could serve as a numerical testbed for hypotheses related to LC hemodynamics. 

Acknowledgements

Supported by grants from the Ministry of Education, Academic Research Funds, Tier 1 (MJAG; R-397-000-181-112; Singapore) and from an NUS Young Investigator Award (MJAG; NUSYIA_FY13_P03, R-397-000-174-133; Singapore). 

Disclosure: T. Chuangsuwanich, None; K.E. Birgersson, None; A. Thiery, None; S.G. Thakku, None; H.L. Leo, None; M.J.A. Girard, None