Characterizing the Collagen Network Structure and Pressure-Induced Strains of the Human Lamina Cribrosa

Abstract

Purpose: The purpose of this study was to measure the 2D collagen network structure of the human lamina cribrosa (LC), analyze for the correlations with age, region, and LC size, as well as the correlations with pressure-induced strains.

Methods: The posterior scleral cups of 10 enucleated human eyes with no known ocular disease were subjected to ex vivo inflation testing from 5 to 45 mm Hg. The optic nerve head was imaged by using second harmonic generation imaging (SHG) to identify the LC collagen structure at both pressures. Displacements and strains were calculated by using digital volume correlation of the SHG volumes. Nine structural features were measured by using a custom Matlab image analysis program, including the pore area fraction, node density, and beam connectivity, tortuosity, and anisotropy.

Results: All strain measures increased significantly with higher pore area fraction, and all but the radial-circumferential shear strain (E_rθ) decreased with higher node density. The maximum principal strain (E_max) and maximum shear strain (Γ_max) also increased with larger beam aspect ratio and tortuosity, respectively, and decreased with higher connectivity. The peripheral regions had lower node density and connectivity, and higher pore area fraction, tortuosity, and strains (except for E_rθ) than the central regions. The peripheral nasal region had the lowest E_max, Γ_max, radial strain, and pore area fraction.

Conclusions: Features of LC beam network microstructure that are indicative of greater collagen density and connectivity are associated with lower pressure-induced LC strain, potentially contributing to resistance to glaucomatous damage.

The lamina cribrosa (LC) is a connective tissue–containing structure that spans the optic nerve head (ONH) through which pass the unmyelinated retinal ganglion cell (RGC) axons as they go from cell bodies in the retina to the optic nerve. The LC network structure mechanically supports the RGC axons to withstand the translaminar pressure difference between the intraocular and extraocular space and the hoop stresses of the adjacent sclera imposed by the intraocular pressure (IOP).¹ In both human glaucoma and experimental models with increased IOP, RGC axonal injury typically occurs at the LC.² Thus, the structure and biomechanical behavior of the LC are clearly important potential pathogenic features of the disease. 

In humans and monkeys, the LC is composed of connective tissue beams whose core contains elastin and collagen types I and III fibrils.³ The beams are organized into a series of stacked cribriform plates, forming a beam and pore network structure. LC beams are covered by astrocytes, while fibroblast-like cells (lamina cribrosacytes) and microglia reside within the beams.^4–6 An intricate network of capillaries consisting of endothelial cells and pericytes also reside within the LC beams.^7,8 The structure of the LC is deformed by the IOP and the translaminar pressure gradient. ONH astrocytes and lamina cribrosacytes have been shown to become reactive to changes in the mechanical environment in animal models and in vitro experiments.^5,9–12 Thus, there are both cellular and noncellular response elements that are thought to be involved in mechanism of glaucomatous axonal damage.^6,13–15 The LC structure also undergoes significant connective tissue remodeling in glaucoma.² This may alter the physiological support for the resident cells and LC capillaries,¹⁶ which are likely contributing factors to glaucoma injury. 

Studies of LC structure using scanning electron microscopy (SEM) suggest that the connective tissue density is lower in the inferior and superior regions of the ONH.¹⁷ The authors speculate that the lower connective tissue density may subject the axons in these regions to higher pressure-induced strains leading to their earlier dysfunction and death. Subsequent studies of the LC pore and beam structure by histology,^18–22 scanning electron microscopy,²³ polarized light microscopy,²⁴ optical coherence tomography (OCT),^25–28 and confocal adaptive optics scanning laser ophthalmoscopy²⁹ have found similar regional variations, as well as variations between the anterior and posterior regions of the LC. Recent experimental and clinical studies applying image correlation methods to map the deformation response of the human ONH to pressure changes^30–32 have shown that mechanical strains within the LC exhibit local regions of high shear, tensile, and compressive strains. Midgett et al.³² have used a laser scanning microscope and second harmonic generation (SHG) imaging to study the collagen beam structure of the LC at a reference pressure and higher pressure. Digital volume correlation (DVC) was used to correlate the image volume at high pressure to the image volume at the reference pressure to calculate the displacement field and strain field of the LC. The maximum shear strain is larger in the peripheral versus central LC, while the maximum principal strain centrally is larger in the inferior and temporal quadrants than in the nasal and superior quadrants. Comparing the quadrants showed that the maximum principal strain was smallest in the nasal quadrant. These findings support the observations that lower levels of connective tissue are associated with greater strain and are consistent with clinical observations that structural and functional glaucoma injury occurs preferentially in RGC axons traversing the inferotemporal ONH.^17,33–36 However, past investigations have not measured the beam and pore network structure of the LC and compared it to the IOP-induced strain in the ONH of the same specimen. 

The aim of this work was to characterize the beam network microstructure of the LC of human donor eyes, including six eyes tested by inflation in Midgett et al.³² along with four additional eyes, and to correlate measured LC strains with the network structure. We developed methods to analyze the image volume at the reference pressure to measure features of the LC beam and pore network structure, including the area fraction of the pores; node number density, where a node is defined as the junction of collagen beams; connectivity, defined as the number of beams joining at a node; beam width, length, aspect ratio, dominant orientation, and degree of alignment. The structural features and strain measures were compared also for the correlations with age, region, and LC size. 

Methods

The specimen preparation, imaging, DVC, and strain calculation methods were described in detail in Midgett et al.³² The following section briefly summarizes these methods and presents in detail methods for image preprocessing, segmentation, and network analysis used to identify and quantify the features of the LC network structure in the imaged volumes, as well as the statistical methods for analysis of age and regional variations, and correlations between measured outcomes. 

Specimen Preparation

Ten human eyes from eight donors with no history of glaucoma were obtained from the National Disease Research Interchange and Eversight 12 to 24 hours post mortem and tested by inflation within 48 hours post mortem (Table 1). These included the six specimens presented previously in Midgett et al.³² and four newly tested specimens. The optic nerves of the enucleated eyes were cut flush to the sclera to remove the post-LC myelinated optic nerve and to reveal the trabecular structure of the LC. The eyes were glued to a custom plastic ring at a location 2 to 5 mm posterior to the equator. The ONH was aligned with the center of the plastic ring. The anterior chamber and intraocular components were removed by trimming along the coronal plane of the eye with a scalpel. 

Table 1

 

Human Donor Posterior Sclerae Subjected to Inflation Testing

 

Image Acquisition

Posterior scleral specimens were mounted in a custom inflation chamber, such that the posterior surface of the LC was aligned with the objective of a Zeiss 710 laser-scanning microscope (LSM 710 NLO; Zeiss, Inc., Oberkochen, Germany). Specimens were kept hydrated with 1M phosphate-buffered saline throughout the inflation test. Pressure was applied by using a water column. Each specimen was pressurized first to a baseline of 5 mm Hg then to 45 mm Hg. At both pressures, the specimen was allowed to equilibrate for at least 30 minutes before two sets of 2 × 2 tiled z-stacks were acquired, starting from a z-depth of 300 μm below the posterior surface and moved toward the posterior surface with a z-increment of 3 μm, using a Chameleon Ultra II laser tuned to 780 nm (Coherent Inc., Santa Clara, CA, USA), a 390- to 410-nm band pass filter to isolate the SHG signal, and a 10× 0.45 NA Apochromat objective. The zoom factor was varied between 0.6 to 0.8 depending on the size of the LC, which corresponded to an in-plane resolution of 2.08 to 2.77 μm/pixel. The tiles were imaged at 512 × 512 pixels and stitched with 15% overlap. The resulting 947 × 947-pixel images were exported as TIFF files (Fig. 1A). 

Figure 1

 

A series of image-processing steps were applied to the SHG image of a representative specimen (specimen 6) at 5 mm Hg to measure the features of the beam network structure, showing (A) the 30th z-slice of the stitched SHG image acquired by the laser scanning microscope after deconvolution, (B) after contrast-limited adaptive histogram equalization, (C) and after application of the Frangi filter⁴⁰ to enhance the contrast of the collagen beams. (D) A 198.9 × 219.3-μm area selected from (C) after binarization using the Otsu thresholding method.⁴¹ (E) Smoothing, dilation, and erosion were performed to obtain the final binary mask used to segment the beam network for structural characterization. The accuracy of the binary mask was confirmed by (F) overlaying the outline of the segmented structure on the preprocessed images. (G) Skeletonization of collagen beams was obtained by using Matlab function ‘bwmorph' option ‘thin'. (H) Nodes were identified as the junctions where two or more skeletonized beams meet. (I) The width of each collagen beam was measured at the midpoint of the beam by searching for black pixels along the direction perpendicular to the measured beam orientation. Scale bar in (A–C): 200 μm.

Image Preprocessing

The SHG volumes (Fig. 1A) were preprocessed by three-dimensional (3D) deconvolution using a theoretical point spread function in Huygens Essentials (SVI, Hilversum, The Netherlands) to reduce noise and blur. Contrast enhancement was then performed by two-dimensional (2D) contrast-limited adaptive histogram equalization³⁷ in FIJI³⁸ (Fig. 1B). 

The SHG image volumes were used for DVC analysis of the 3D displacement field caused by inflation from 5 to 45 mm Hg, using the Fast Iterative DVC algorithm by Bar-Kochba et al.^32,39 The maximum z-projection of the SHG volume at 5 mm Hg was used to generate a 2D LC mask by manually tracing the boundary of the peripapillary sclera on the z-projection image, using the Matlab function ‘roipoly' (Matlab R2017a; MathWorks, Inc., Natick, MA, USA). The sclera is denser in collagen than the LC and was clearly identifiable as an oversaturated region in the SHG imaged volumes. The number of pixels within the traced LC boundary of each specimen was recorded as the LC area. The preprocessed images acquired at 5 mm Hg were further analyzed to measure 2D features of the LC network structure as described below. Unless stated otherwise, all following image processing and structural characterizations were performed on 2D z-slices. 

Segmentation and Skeletonization of the LC Network Structure

A modified 3D Frangi vessel filter,⁴⁰ with settings α = 0.05, β = 1, and γ = 7.5, was applied to the preprocessed image stacks at 5 mm Hg to further enhance the contrast of the beam structures (Fig. 1C).^41,42 The Frangi filter vesselness threshold was set to 1/6 of the brightness of a representative beam,⁴³ and the scale range was set to 0.1 to 7.1 with a step size of 1 pixel to account for a wide distribution of beam widths within the LC. The 2D LC mask was applied uniformly through the Frangi-filtered z-stack to segment the LC. A few specimens exhibited locally dark areas with no visible beams in select z-slices. These were further masked and removed from the image volumes to avoid errors in calculating the pore area fraction. 

A 3 × 3 median filter was applied, replacing the intensity of each pixel with the median intensity of the surrounding 9 pixels, to locally smooth the image volume. The resulting volumes were binarized by using the Otsu thresholding method⁴¹ to segment the collagen beams from their background (Fig. 1D). The Otsu thresholding method separates gray-scale images into foreground (white) and background (black) by calculating an optimal thresholding value that minimizes the between-class variance. Here, the white pixels designated the collagen beams and black pixels represented the pores. 

Following the method described in D'Amore et al.,⁴⁴ a series of morphologic operations was performed on the thresholded images in Matlab to remove artifacts generated by binarization, which included removing isolated pixels. The thresholded stacks were first treated with the Matlab function ‘bwmorph' using the options ‘thin', ‘majority', ‘clean' to remove pixel-level noise. Morphologic erosion (Matlab function ‘imerode'), dilation (Matlab function ‘imdilate'), and a final erosion step were applied to smooth the beam boundaries (Fig. 1E). Disk element sizes of 1/6 and 1/3 of the representative beam width were used for the erosion and dilation operations, respectively. A representative beam width was first estimated by counting the number of pixels spanning the width of 10 random beams on the 30th z-slice of specimen 6. The representative beam width of 8 pixels was then determined by iteratively increasing or decreasing the estimation by a step size of 1 pixel, and verifying the results by overlaying the boundaries of the segmentation on the preprocessed images (Fig. 1F). To quantitatively verify the accuracy of the segmentation, the binarization from autosegmentation was compared to that from manual tracing by calculating the percentage differences in the resulting pore area fractions. The calculation of pore area fraction is described in detail in the Structural Measurements section. Two local regions of size 200 μm × 200 μm on the 30th z-slice of all specimens were selected for manual segmentation. Resulting pore area fraction from representative beam width of 7 and 9 pixels were also calculated for validation. The iterative process ensured that the width of the collagen beams was minimally affected by the image-processing steps. The segmented beam network in each z-slice (Fig. 1E) was skeletonized by using the Matlab function ‘bwmorph' with options ‘thin' and ‘inf' (Fig. 1G). Skeletonization shrinks the beams to 1 pixel in width. Two-dimensional skeletonization was used, since available 3D skeletonization methods for the plate-like meshes produce spurious branches that are not representative of the LC structure. 

Regional Division

The LC was divided into eight regions centered about the central retinal artery and vein (CRAV) (Fig. 2). A CRAV center point was manually picked by marking a point on the maximum intensity projection of the SHG images. The CRAV region was defined as a cylindrical volume of radius of 200 μm about the CRAV center. The eight regions were defined by dividing the LC into central and peripheral regions at the radial midpoint between boundary of the CRAV and the LC mask. The central and peripheral regions were further divided into four quadrants—superior (S), inferior (I), temporal (T), and nasal (N)—using 45° and 135° bisectors. 

Figure 2

 

Regional division of the LC for statistical analysis. The center of the CRAV was manually identified on the maximum intensity projection of the SHG volumes acquired at 5 mm Hg. The CRAV was defined as the cylindrical region with radius of 200 μm from the CRAV center. The LC was divided into eight regions centered about the CRAV by first dividing the LC into central (Cen-) and peripheral (Peri-) regions at the location that was halfway between boundary of the CRAV and the LC mask. The central and peripheral regions were further divided into four quadrants—superior (S), inferior (I), temporal (T), and nasal (N)—by using 45° and 135° bisectors.

Structural Measurements

Nine structural features were identified and calculated from the segmented beam and pore structures (Fig. 1E) and the skeletonized network structure (Fig. 1G). The collagen structures appeared sparsely in the top-most z-slices because of the roughness of the specimen cut surface. The bottom-most slices were also sparse in features because of light attenuation. We only used the z-slices in which at least 2/3 of the LC was visible to measure the structural features. For the 100-slice (300 μm) deep image volume, this included the central 16 to 26 slices depending on the specimen. The structural features were measured on each z-slice, then averaged either within the eight regions (Fig. 2) or within the LC masked volume for statistical analysis. 

The structural features were defined as follows. 

Pore Area Fraction

The pore area fraction was calculated as the number of black pixels within the area of interest (either the LC masked area or region) by the total number of pixels in that area (Fig. 1E). A pore area fraction of 1 represents an LC with no collagen beam structures. 

Node Density

The pixels where two or more skeletonized beams (representing collagen beams) intersect in the skeletonized network were identified initially as a node (Matlab function ‘bwmorph', option ‘branchpoints') (Fig. 3B). Each node was dilated to a disk size of the representative beam width (Figs. 1H, 3C) to merge nodes that were closer to each other than a representative beam width. This step was necessary to remove an artifact of skeletonization, where wide beam junctions were separated into two nodes (Figs. 3A, 3B). The node density was defined in the area of interest as the number of nodes divided by the area. 

Figure 3

 

Illustration of the method used for node detection and connectivity measurement, showing (A) the skeletonized collagen beams within a 48.5 × 35.7-μm area from Figure 1g. (B) Nodes were initially identified at every junction of the skeletonized beams. (C) Each node was dilated to a disk to merge nodes that were closer to each other than the representative beam width. (D) The connectivity of a node was calculated as the number of skeletonized beams extruding from the node. The gray region represented the dilated node and the green region represented pixels around the boundary of the node that were used to search for a connecting beam represented by a white pixel.

Connectivity

The connectivity was defined as the number of network beams connected to a node. This was calculated by searching for pixels within a 1-pixel radial distance immediately surrounding the boundary of each dilated node and recording the number of white pixels (highlighted in dark green in Fig. 3D). 

Beam Length

The pixels representing each dilated node (highlighted in gray in Fig. 3D) were set to 0 (black) to separate the connected skeletonized beams. The straight line distance of the beam was calculated as the euclidean distance between the starting and ending pixels of each disconnected skeletonized beam (Matlab function ‘pdist2’). 

Beam Tortuosity

The beam tortuosity was defined as the ratio of the contour length divided by the straight line distance of the skeletonized beam. The contour length of a beam was evaluated as the sum of the euclidean distance between each pixel in the skeletonized beam (Matlab function ‘bwdistgeodesic' option ‘quasi-euclidean'). 

Beam Orientation

An ellipse was fitted to each skeletonized beam, where the major axis was identified. The angle between the major axis and the horizontal, nasal-temporal, axis of the LC was defined as the beam orientation, which ranged from −π/2 to π/2 in the X-Y plan (Matlab function ‘regionprops' option ‘orientation'). A preferred beam orientation for an area of interest was calculated as the circular average of beam orientation weighted by the beam length, using the Matlab circular statistics toolbox developed by Berens.⁴⁵ 

Anisotropy

The anisotropy of the beams was defined for an area of interest by fitting the semicircular von Mises probability density function to the distribution of beam orientation angles relative to the areal average beam orientation. The resulting dispersion parameter, which signifies the degree of alignment along the average orientation, was defined as the anisotropy. 

Beam Width

The width of each collagen beam was measured at the midpoint of the skeletonized beam. The beam angle at the midpoint was recorded as the angle between the line connecting points on the skeletonized beam that were 2 pixels away from the midpoint and the nasal-temporal axis. A line perpendicular to the beam angle was drawn at the midpoint, and the beam width was defined as the number of white pixels on the perpendicular line extending in both directions from the midpoint until it reached a black boundary (Fig. 1I). 

Beam Aspect Ratio

The aspect ratio of a beam was determined as the ratio of the beam length over the width. 

Strain Calculations

The methods for DVC analysis, DVC error analysis, and displacement post processing for strain calculation have been described in detail in previous studies.^32,46 In the following, we briefly outline the method for strain calculation. The Fast Iterative DVC algorithm, originally developed by Bar-Kochba et al.,³⁹ was applied to the preprocessed SHG volumes at 5 mm Hg and 45 mm Hg to determine the components of the 3D displacement field, U_X (X, Y, Z), U_Y (X, Y, Z), and U_Z (X, Y, Z). The (X, Y, Z) directions are aligned with the nasal-temporal, inferior-superior, and anterior-posterior axes, respectively (Fig. 2). The components of the Green-Lagrange strain tensor E in the Cartesian (X, Y, Z) coordinate system were calculated from the displacement gradients as  

\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\begin{equation}\tag{1}{E_{XX}} = {{\partial {U_X}} \over {\partial X}} + {1 \over 2}\left( {{{\left( {{{\partial {U_X}} \over {\partial X}}} \right)}^2} + {{\left( {{{\partial {U_Y}} \over {\partial X}}} \right)}^2} + {{\left( {{{\partial {U_Z}} \over {\partial X}}} \right)}^2}} \right)\end{equation}

 

\begin{equation}\tag{2}{E_{YY}} = {{\partial {U_Y}} \over {\partial Y}} + {1 \over 2}\left( {{{\left( {{{\partial {U_X}} \over {\partial Y}}} \right)}^2} + {{\left( {{{\partial {U_Y}} \over {\partial Y}}} \right)}^2} + {{\left( {{{\partial {U_Z}} \over {\partial Y}}} \right)}^2}} \right)\end{equation}

 

\begin{equation}\tag{3}{E_{XY}} = {1 \over 2}\left( {{{\partial {U_X}} \over {\partial Y}} + {{\partial {U_Y}} \over {\partial X}} + {{\partial {U_X}} \over {\partial X}}{{\partial {U_X}} \over {\partial Y}} + {{\partial {U_Y}} \over {\partial X}}{{\partial {U_Y}} \over {\partial Y}} + {{\partial {U_Z}} \over {\partial X}}{{\partial {U_Z}} \over {\partial Y}}} \right).\end{equation}

 

The strain components E_XX and E_YY signify the elongation or compression in the nasal-temporal direction and inferior-superior direction, respectively, and E_XY signifies the angle distortion between the X-Y directions. We calculated the out-of-plane strain components, E_ZZ, but do not report them here because preliminary analysis showed large errors for E_ZZ (Supplementary Section S3). The in-plane strain components were used to calculate the maximum principal strain, E_max, which signifies the tensile strain along the direction in the X-Y plane at which it is maximum; and the maximum shear strain, Γ_max, which signifies the shear strain in the direction of the X-Y plane at which it is maximum:  

\begin{equation}\tag{4}{E_{max}} = {{{E_{XX}} + {E_{YY}}} \over 2} + \sqrt {{{\left( {{{{E_{XX}} - {E_{YY}}} \over 2}} \right)}^2} + E_{XY}^2} \end{equation}

 

\begin{equation}\tag{5}{{\rm{\Gamma }}_{max}} = \sqrt {{{\left( {{{{E_{XX}} - {E_{YY}}} \over 2}} \right)}^2} + E_{XY}^2} .\end{equation}

 

We also calculated the cylindrical components of strain, that is, the radial strain, E_rr, the circumferential strain, E_θθ, and the radial-circumferential shear strain, E_rθ, using coordinate transformation because of the cylindrical shape of the ONH. To calculate the cylindrical coordinates of each point in the ONH, we assumed that the origin of the cylindrical coordinate system was the center of the ellipse fit to the boundary of LC, using the Matlab function ‘fit_ellipse' (Ohad Gal, 2003). An orientation angle θ was determined for each point in the LC area by defining a line from the pixel to the center of the ellipse and calculating the angle between the line and the horizontal X axis. The angle θ was used to transform the strain components as follows:  

\begin{equation}\tag{6}{E_{rr}} = {E_{XX}}{\cos^2}\theta + 2{E_{XY}}\cos \theta \sin \theta + {E_{YY}}{\sin^2}\theta \end{equation}

 

\begin{equation}\tag{7}{E_{\theta \theta }} = {E_{XX}}{\sin^2}\theta - 2{E_{XY}}\cos \theta \sin \theta + {E_{YY}}{\cos^2}\theta \end{equation}

 

\begin{equation}\tag{8}{E_{r\theta }} = \left( {{E_{YY}} - {E_{XX}}} \right)\cos \theta \sin \theta - {E_{XY}}\left( {{\sin^2}\theta - {\cos^2}\theta } \right).\end{equation}

 

We estimated the DVC displacement correlation error and strain error for each specimen (Supplementary Table S3) by applying numerically a rigid body displacement and triaxial strain to one of the duplicate z-stacks acquired at 5 mm Hg, correlating the deformed z-stack to the undeformed z-stack using DVC, and calculating the difference between the DVC correlated and the numerically applied displacements and strains (see Midgett et al.³² for more details). The absolute displacement Display Formula\({e_{{U_{I_j}}}}\)and strain Display Formula\({e_{{E_{I{J_j}}}}}\) error between the correlated displacement U_I and strain E_IJ components and the applied displacement Display Formula\(U_I^{app}\)and Display Formula\(E_{IJ}^{app}\)strain components are defined at each point j in the LC as  

\begin{equation}\tag{9}{e_{{U_{I_j}}}} = \sqrt {{{\left( {U_{I_j}^{app} - {U_{I_j}}} \right)}^2}} \end{equation}

 

\begin{equation}\tag{10}{e_{{E_{I{J_j}}}}} = \sqrt {{{\left( {E_{I{J_j}}^{app} - {E_{I{J_j}}}} \right)}^2}} ,\end{equation}

where I and J are X , Y, and Z. This error includes the effects of tissue creep during imaging, characteristics of the collagen structure used as the natural speckles for DVC and speckle distortion by deformation. The average absolute displacement error was under 2.43 ± 0.95 μm (0.88 ± 0.34 pixels) for U_X and U_Y and under 5.256 ± 3.232 μm (1.90 ± 1.17 pixels) for U_Z for all specimens. The absolute strain error for each specimen was below 0.003 for E_XX, E_YY, and E_XY, and below 0.045 for E_ZZ (Supplementary Table S3).  

Statistical Analysis

The nine structural features; strain components E_rr, E_θθ, and E_rθ; maximum principal strain E_max; and maximum shear strain Γ_max were averaged for each specimen. Circular mean was used for average beam orientation. Linear regression was used to estimate the correlations with age and with LC area for each averaged outcome, and to test correlations between the nine structural features, as well as correlations between the structural features and five strain measures. The analysis was performed in Matlab using the function ‘fitlm'. The correlation between left and right eyes was not taken into account in the linear regression analysis because of the small sample size (n = 10) and because only two donors out of eight contributed data from both eyes. One-way analysis of variance (ANOVA) was performed by using the Matlab function ‘annova1’ and followed by ‘multcompare' for multiple comparison test of the specimen-averaged strain outcomes. 

The strains and structural outcomes were also averaged in the eight regions of each specimen, resulting in 80 regional measurements for each outcome. Six regional measurements were found to be extreme outliers and were eliminated for the regional analyses. The extreme outliers were greater than the third quartile for the measured outcome by more than three times the interquartile range. Linear mixed models were used to take into account the clustering of two eyes from a single donor and the repeated measurements within a single eye for the regional analyses. Measurements from each of the eight regions of the LC were assumed to have a compound symmetry correlation structure, in which the measurements from any two regions have the same correlation. Least square means were used to estimate the mean outcome for the central region, peripheral region, full quadrants, central quadrants, and peripheral quadrants. The Tukey-Kramer method was used to adjust the significance levels upwards to reduce the probability of a type I error to less than 0.05 for the multiple pairwise comparisons.⁴⁷ The P values that were adjusted were labeled as “adjusted P values” in all following tables. The mixed models were used to analyze for (1) the differences in means of each of the five strains and nine structural features between central and peripheral regions, (2) the differences in means of strains and structural features between the quadrants, and (3) the correlations between the nine structural features and five strain measures. The results of (3) were further analyzed to determine which structure has the greatest effect on E_max and Γ_max. Structural features significantly associated with the average E_max and Γ_max strains were first identified. For each of these structures, the absolute change in strain for a change of 1 standard deviation in structure was calculated from the estimated mixed model parameters and the standard deviation of the regional measurements of structure. For each strain, the mixed models for the selected structures were fit by using the same set of observations, and the Akaike Information Criteria (AIC) from the mixed models for a strain were then used to identify the structure that best predicts the observed differences in strain. Regional analyses were performed with SAS 9.2 (SAS Institute, Cary, NC, USA). For all analyses, a P value or an adjusted P value lower than 0.05 was considered as a significant observation. 

Results

Strain Outcomes

The strains showed little variation through the thickness of the imaged volume (through Z) but were highly variable within the nasal-temporal (X-Y) plane (Fig. 4, Supplementary Figs. S1–S10). The normal strain components E_rr and E_θθ were predominantly positive (tensile) and concentrated at the boundary with the peripapillary sclera. The shear strains Γ_max and E_rθ also exhibited localized regions of large positive and negative strains at the boundary. For the 10 specimens, the average radial strain, 0.027 ± 0.014, was similar in magnitude to the circumferential strain, 0.024 ± 0.014 (P = 0.99). The radial-circumferential shear strain E_rθ, 0.0026 ± 0.004, was one order of magnitude lower than the radial and circumferential strain components (P < 0.005), though it had small, local regions of high negative and positive strains. The specimen-averaged maximum shear strain Γ_max, (0.018 ± 0.007), was significantly lower than the maximum principal strain E_max (0.044 ± 0.018, P = 0.0005), which agrees with prior studies,^32,46 though all specimens exhibited localized regions of high shear. 

Figure 4

 

The thickness-averaged strain response for inflation from 5 to 45 mm Hg for specimen 6 showing (A) the maximum principal strain E_max, (B) maximum shear strain Γ_max, (C) radial strain E_rr, (D) circumferential strain E_θθ, and (E) radial-circumferential shear strain E_rθ. (F) A box plot of the specimen-averaged strains.

Structural Features

The average pore area fraction for the 10 specimens was 0.43 ± 0.05 μm²/μm², while the specimen-averaged connectivity ranged from 3.03 to 3.18 beams per node. The beams were relatively straight, with an average tortuosity of 1.13 ± 0.02 μm/μm and had a low aspect ratio (length/width) of 1.59 to 2.13 μm/μm. The average beam orientation angle was −2.22° ± 2.77° with respect to the nasal-temporal axis. However, the beam orientations were highly dispersed with an average anisotropy of 0.54, where a value of 0 indicates a random distribution and infinity indicates that all the beams are aligned along the average orientation. 

Linear regression models were used to examine for pairwise correlations between the specimen-averaged structural features (Table 2, Supplementary Table S6). A higher pore area fraction was associated with fewer beams joining at a node (P = 0.0002), a higher beam aspect ratio (P = 0.04), and a more negative beam angle with respect to the nasal-temporal axis (P = 0.02). A longer beam was associated with a lower node density (P = 0.001), a more aligned beam (P = 0.04), and a higher beam aspect ratio (P = 0.02). 

Table 2

 

Results of the Linear Regression Analysis for the Correlation Between the Specimen-Averaged Structural Features (n = 10), Showing Estimated Regression Parameter, the P Value for the Significance of the Structural Correlation and the R² Correlation Coefficient

 

Variations With Age

The structural features did not vary significantly with age (Supplementary Figs. S12, S13). The radial strain E_rr, circumferential strain E_θθ, maximum tensile strain E_max, and maximum shear strain Γ_max all decreased significantly with age, while E_rθ, which was an order of magnitude smaller than E_rr and E_θθ, increased with age (Fig. 5; Supplementary Table S8). 

Figure 5

 

Linear regression analysis of the specimen-averaged strain outcomes for variation with age, showing (A) E_rr radial strain, (B) E_θθ circumferential strain, and (C) E_rθ radial-circumferential shear strain.

Variations With LC Area

The LC area for the 10 specimens ranged from 2.61 to 4.14 mm². The pore area fraction (P = 0.001) and beam aspect ratio (P = 0.001) increased with larger LC area, while connectivity decreased with larger LC area (Fig. 6, Supplementary Fig. S15). Interestingly, higher normal (tensile and compressive) strains, E_rr, E_θθ, and E_max, were also associated with higher LC area, while the shear strains, Γ_max and E_rθ, did not vary significantly with the LC area (Fig. 7, Supplementary Fig. S14). 

Figure 6

 

Linear regression analysis of the specimen-averaged structural features for variation with LC area, showing (A) pore area fraction, (B) connectivity, and (C) beam aspect ratio.

Figure 7

 

Linear regression analysis of the specimen-averaged strain outcomes for variation with LC area, showing (A) E_max maximum principal strain and (B) E_rr radial strain.

Regional Variations

Radial Variations

Compared to the central region, the peripheral region had a 15% higher pore area fraction (P < 0.0001), 2% higher tortuosity (P = 0.01), 21% lower node density (P < 0.0001), and 3% lower connectivity (P < 0.0001) (Fig. 8, Supplementary Fig. S16). The maximum principal strain E_max was 37% higher (P < 0.0001), the maximum shear strain Γ_max was 46% higher (P < 0.0001), and the radial strain was 25% higher (P = 0.04) in the peripheral than in the central LC (Fig. 9). The peripheral region also exhibited higher E_θθ (P = 0.02). The radial-circumferential shear strain E_rθ was more negative in the peripheral region (P = 0.02). 

Figure 8

 

Effects of radial positions on structural features, showing (A) pore area fraction, (B) node density, (C) connectivity, and (D) tortuosity.

Figure 9

 

Effect of radial position on strain outcomes, showing (A) E_max maximum principal strain, (B) Γ_max maximum shear strain, (C) E_rr radial strain, (D) E_θθ circumferential strain.

Quadrant Variations

The strain and structural variations were analyzed also for variations between the four quadrants within the central, peripheral, and full LC (Supplementary Tables S13–S18). 

In the full LC, the beam orientation (P = 0.02), anisotropy (P = 0.01), and length (P = 0.02) varied significantly among the quadrants (Supplementary Table S13). Multiple post hoc pairwise comparisons showed that the beam orientation angle was smaller (i.e., more negative with respect to the nasal-temporal direction) in the superior quadrant than the inferior quadrant (P = 0.02). Beam anisotropy in the superior (P = 0.01) quadrants was smaller (less aligned) than that in the inferior quadrant, and beam length was shorter in the temporal quadrant than in the inferior quadrant (P = 0.01). However, the differences in anisotropy were small, ranging from 0.56 to 0.87. The maximum principal strain E_max, maximum shear strain Γ_max, and radial-circumferential shear strain E_rθ also varied significantly among the quadrants (Supplementary Table S16). For all strain outcomes except E_rθ, the strain was highest in the temporal quadrant and lowest in the nasal quadrant (P = 0.03 for E_max and P = 0.01 for Γ_max from multiple pairwise comparisons). 

Within the central LC, only the beam anisotropy varied significantly among the quadrants (P = 0.004) and was highest in the inferior quadrant and lowest in the superior quadrant (P = 0.01, Supplementary Table S14). The maximum principal strain E_max (P = 0.02), maximum shear strain Γ_max (P = 0.03), and radial-circumferential shear strain E_rθ (P = 0.03) also varied among the central quadrants. The E_max (P = 0.02) and Γ_max (P = 0.04) were higher in the inferior region than nasal region. E_max was also higher in inferior than superior region (P = 0.05). 

In the peripheral LC, the pore area fraction varied significantly between quadrants (P = 0.05), and was smaller in the nasal than superior quadrant (P = 0.04, Fig. 10A). The nasal region also had lower E_max (P = 0.001), Γ_max (P < 0.001), and E_rr (P = 0.02) than the temporal region (Fig. 10B; Supplementary Table S18). 

Figure 10

 

Effect of peripheral quadrant position in (A) pore area fraction and (B) maximum principal strain. The P values were calculated from linear mixed models. The superscript asterisk (*) indicates significant adjusted P value using Tukey-Kramer method for multiple pairwise comparison.

Structure-Strain Correlations

Specimen Averages

Linear regression models were used to investigate the pairwise correlations between specimen-averaged strain and structural outcomes (Table 3, Supplementary Tables S19–S23). The maximum principal strain E_max increased with higher pore area fraction (P = 0.01). The circumferential strain also increased with larger pore area fraction (P = 0.006) and beam aspect ratio (P = 0.004). The radial strain increased with larger pore area fraction (P = 0.007) and tortuosity (P = 0.03), and decreased with higher connectivity (P = 0.03). The strongest correlations (highest R² values) were between the radial strain E_rr and the pore area fraction and between the circumferential strain E_θθ and the beam aspect ratio (Fig. 11). 

Table 3

 

Results of the Linear Regression Models for Correlations Between Structures Measured at 5 mm Hg and Strains Measured for Inflation From 5 to 45 mm Hg Using Specimen-Averaged Data, Showing the Estimated Regression Parameter, the P Value for the Significance and the R² Correlation Coefficient

 

Figure 11

 

The specimen-averaged radial strain E_rr increased with (A) larger pore area fraction and (B) larger tortuosity. The specimen-averaged circumferential strain E_θθ increased with (C) larger pore area fraction and (D) higher beam aspect ratio (n = 10).

Regional Averages

We next applied linear mixed models to investigate the pairwise correlations between regionally averaged strain and structural outcomes, because significant regional variations were found for both structural and strain outcomes (Table 4, Supplementary Tables S24–S28). All strain components increased significantly with increased pore area fraction and all strains (except for E_rθ) decreased with higher node density. Moreover, E_max and Γ_max decreased with higher connectivity; E_max and E_θθ increased with higher beam aspect ratio; and Γ_max increased with higher tortuosity (Table 4). The shear strain E_rθ decreased with increasing beam orientation angle and increased with longer beams. 

Table 4

 

Results of the Linear Mixed Models for Correlations Between Regionally Averaged Structural Features Measured at 5 mm Hg and Strains Measured for Inflation from 5 to 45 mm Hg, Showing the Number of Measurements After Elimination of Outliers, the Estimated Regression Parameter and the P Value for the Significance

 

A total of four structural measures were significantly correlated with E_max and Γ_max. These were compared to identify the structure with the greatest impact on the average strain. AIC statistics were used to measure the information lost when the model is used to predict the observed E_max and Γ_max strain. This allowed for calculation of the relative likelihood for each model, which is the statistical probability that the model will best predict the strain outcome, relative to the best fitting model. For E_max and Γ_max, the structures that had the greatest estimated absolute change in strain corresponding to a change of 1 standard deviation in structure were pore area fraction and connectivity, respectively (Fig. 12). The predicted strains using these structures also best fit the observed strains. 

Figure 12

 

(A) The maximum principal strain E_max increased with larger pore area fraction (n = 80); (B) the maximum shear strain Γ_max increased with lower connectivity (n = 79); (C) the radial strain E_rr (n = 80) and (D) the circumferential strain E_θθ (n = 79) decreased with larger node density.

Discussion

Our method for measuring the LC microstructure from SHG image volumes of the human LC produced quantitative findings similar to those previously reported by other methods. The average beam width was 36.23 ± 2.65 μm, which is comparable to previous beam width measurements by OCT (38.1 ± 1.4 μm by Nadler et al.⁴⁸ and 46.7 ± 3.2 μm by Wang et al.²⁸). The average pore area fraction, 0.43 ± 0.05, also was similar to the 0.48 value for proportionate pore area for the posterior LC in histologic sections,¹⁸ but somewhat higher than the 34.5% ± 7.6% measured by SEM.²³ The average tortuosity of the collagen beams (1.13 ± 0.02) was also larger than the range of tortuosities (1–1.12) measured by Brazile et al.⁴⁹ for the collagen fibers within the LC beams of sheep eyes. While our specimens had been inflation tested, they were not chemically fixed and remained in a fully hydrated state, compared to histologic sections or dried scanning microscopy specimens. Differences in the structural measures between methods could also be caused by variations between species, differences in the imaging and segmentation approaches, and differences in the specimen preparation methods. We confirmed that the segmentation method was appropriate by comparing the difference in resulting pore area fraction from manual segmentation and that from the Matlab algorithm. Two local regions of size 200 μm × 200 μm on the 30th z-slice of all specimens were selected for manual segmentation. For a representative beam width of 7, 8 (used in this study), and 9 pixels, the average absolute percentage difference in resulting pore area fraction when compared to manual segmentation was 16.91%, 3.87%, and 37.97%, respectively; while the average difference between two manual tracings by the same operator on 10 of the regions was 4.37%. This showed that despite an overestimation or underestimation that may be uniformly applied across images during autosegmentation, the approach eliminates human error in manual identification especially in regions with lower beam contrast. In addition, we compared the resulting pore area fraction from the three representative beam widths (7, 8, and 9 pixels) with the maximum principal strain E_max in our linear regression model. The change in parameter did not affect the conclusion that specimen-averaged pore area fraction increased with increasing pore area fraction. Overlaying the segmented structures onto the original SHG image volume also served to visually verify that the segmentation method was sound. 

There were significant correlations between LC structural features and strains measured in the same specimens, both globally and regionally. Strains were larger with higher pore area, lower beam connectivity, and more tortuous, thinner beams. These findings show that an LC network structure with a lower connective tissue content and fewer beam connections is associated with greater IOP-induced strains. All strain measures increased significantly with larger pore area, and greater pore area and greater connectivity were the most predictive factors determining regional maximum principal and maximum shear strains, respectively. Our measurements reflected the tissue structure and mechanical behavior at a single point in time for each LC. It is reasonable to speculate that areas with larger tensile and shear strains would induce stretch activation of astrocytes and lamina cribrosacytes that reside on and within the LC beams. Such stimulation would potentially cause regionally greater connective tissue remodeling of the beam structure, with either protective or adverse effects on RGC axons passing through those zones.^50–53 

We found correlations among LC structural features that have not been detected previously. Lower connectivity and longer beam length were associated with larger pore area fraction, while a longer beam length was associated with a higher degree of beam alignment. Regions with higher pore area fraction likely have longer, but fewer beams to form junctions and connections, and regions with longer but fewer beams would increase anisotropy. We also observed that longer beams were associated with a lower node density and a higher beam aspect ratio. 

In both the analysis of the specimen-averaged and regionally averaged outcomes, variations in structural features that would lead to a more compliant network structure were associated with greater strains. Different structural features affected the normal strains, E_rr and E_θθ, and shear strain E_rθ. We found significant differences between the central and peripheral LC regions for four structural features and for all five strain measures. The pore area fraction was 15% larger in the peripheral region, which agrees with the findings of lower connective tissue volume fraction in the periphery in previous studies using SEM and histology.^17,54 We also found that the peripheral LC had more curved beams and lower connectivity and node density than the central LC, which may contribute to its more compliant network structure and higher strains. Strains and structural features also covaried among LC quadrants. In the peripheral LC, the nasal quadrant had the lowest E_max, Γ_max, and E_rr strains; the lowest pore area fraction; and highest node density. Lower pore area fraction and higher node density both suggest a stiffer network structure, which can lead to lower strains in the peripheral nasal region of the LC and may promote the resilience of the nasal neural-retinal rim in advanced glaucoma.^34,35 While the strains were more uniform in the central LC, its inferior quadrant had the highest E_max and Γ_max. The beam alignment (anisotropy) was highest in the central inferior region. Computational modeling is needed to evaluate the extent that the small, but significant variation in the anisotropy contributed to the more compliant strain response of the central inferior quadrant. 

The correlations between LC structure and strain suggest that these parameters may be useful predictors of the risk and progression of glaucoma damage. As described above, the zones of greater RGC axon damage correspond to the inferior and peripheral LC in glaucoma. The regional findings of our analysis are consistent with this greater susceptibility. In vivo measurements of the pore area fraction and beam width are becoming more accessible with the developments of adaptive optics scanning laser ophthalmoscopy,²⁹ swept-source OCT,²⁸ and adaptive optics spectral-domain OCT.^55–57 These advanced imaging methods permit high-quality visualization of the LC structure, including beam thickness and pore geometry and depth variations of these features.⁵⁸ 

However, variations in structural features could not explain all of the observed strain variations. All strain components varied significantly with age, but none of the structural features varied with age. While the differences in the structural features between the central and peripheral LC were modest (3%–20%), the strains differed by 25% to 50%. Larger differences were also measured for strains than for structural features among the four LC quadrants. These suggest that factors other than LC microstructure as measured by the present methods contributed to variations in strains. The increased stiffness of the LC may be caused by an increase in the modulus of the extracellular matrix material of the LC, for example, by the loss of glycosaminoglycans,⁴⁶ accumulation of elastin and collagen with age,⁵⁹ or by age-related cross-linking of the collagen structure. Previous ex vivo inflation tests have also shown that the inflation response (pressure-strain) of the sclera became stiffer with age.⁶⁰ The anisotropy of the collagen structure of the posterior sclera also decreased with age, and inverse finite element analysis using specimen specific collagen structure and scleral shape showed that the elastic modulus of the sclera, specifically of the components not associated with the aligned collagen fibers, increased with age.^61,62 While a similar inverse analysis is needed to determine the material properties of the LC, the findings here suggest that the age-related stiffening of the pressure-strain response of the LC is not caused by changes in the pore and beam microstructure of the LC. In addition, the LC strain is influenced by the behavior of the parapapillary sclera. The parapapillary sclera is thicker than the LC, has a higher connective tissue density, and has a highly aligned circumferential collagen structure. The mismatch between parapapillary sclera and LC in response to IOP change likely induces a strain concentration at the scleral junction that produces a significant difference between the central and peripheral LC strain independent of LC internal structure.^63,64 

Another parameter that affected the relationship between LC structural elements and strains was the LC area, which varies dramatically among human eyes and has been shown to be an epidemiologic risk factor for glaucoma prevalence, with larger discs at greater risk.^65,66 In larger LCs, there were larger values for pore area fraction, tortuosity, and beam aspect ratio, and lower connectivity. Each of these tendencies was associated with higher strains. Our previous study has not found a statistically significant variation between the strain and LC area,³² perhaps because of the smaller sample size. 

There were several limitations to this study. The strain and structural features were measured from SHG images of the LC collagen structure, which suffers from significant blurring in the Z direction. We applied a series of image processing techniques to reduce the blur and noise and to enhance contrast, but the features remained elongated in the anterior-posterior direction of the LC. This resulted in higher errors for the out-of-plane compressive strain E_ZZ, and E_ZZ was excluded from the analysis for age and region variations and for strain and structure correlations. We estimated the DVC displacement and strain error for every specimen to ensure that the strain variations measured here were all larger than the specimen-averaged absolute error (Supplementary Section S3). The blurring in Z should not have affected the measurement of the beam structures, because we assumed that the beams mainly lie in the plane. 

The LC mask was segmented by using the maximum z-projection image of each LC volume, which could result in a smaller LC area and hide the beams in the peripheral LC from the microstructural analysis. We then applied a 2D skeletonization method to each z-slice to generate the skeletonized beam network that was used to calculate the node density, beam length, aspect ratio, tortuosity, connectivity, orientation, and anisotropy. Histology, polarized light microscopy, and SHG images of LC longitudinal and cross sections of normal human eyes showed that the beams mainly lie in the plane.^67–69 However, a small out-of-plane orientation of the beam may lead to a shorter beam length and other errors in the connectivity and pore density. The LC undergoes significant remodeling in glaucoma and can appear deeply cupped in donors with advanced glaucoma. The present study examined eyes that had no history of glaucoma. The 2D skeletonization method may underestimate the beam length, orientation, anisotropy, and tortuosity in the LC of glaucoma eyes. 

The specimens contained one eye from a Hispanic donor. There may be racioethnic differences in the structural features and biomechanical responses of the LC between the Hispanic donor and Caucasian donors. We repeated the statistical analysis for correlations between the specimen-averaged strain outcomes and pore area fraction and beam aspect ratio. Removing the Hispanic donor eye changed the P values but did not alter the significant findings of this analysis. 

While structural features were not associated with age for the 10 specimens in the age range of 26 to 90 years, we found a borderline significant (P = 0.06) decrease in the specimen-averaged tortuosity with increasing age. The comparison may become significant with more specimens from younger donors. 

Finally, the biomechanical response could be affected by the thickness of the LC after enucleation. If a sizable stump was left on the inflation-tested specimen, the dura may increase the stiffness of the sclera canal opening and the RGC axons may constrain the posterior deflection of the LC under inflation. We removed the optic nerve up to the LC posterior surface to image the LC structure for strain and structural measurements. To ensure that we did not remove a significant portion of the LC, the specimen was imaged after each thin cut under a dissecting microscope during the specimen preparation process. The SHG method was able to image through 300 μm of the posterior LC surface, but the more anterior and posterior z-slices had large DVC error and unreliable structural measurements, either because of surface roughness or light attenuation, and were excluded from the average strain and structure calculations. For all specimens, the most anterior slice used for structural characterization was between 225 to 270 μm from the posterior surface, and the most posterior slice used was between 150 to 210 μm from the posterior surface. Previous histologic studies have shown that the number of pores increases and the pore size decreases anterior-posteriorly.¹⁸ There may be differences in structural parameters between more anterior and more posterior sections that were not captured within the slices selected for structural analysis. Further, the thickness of LC was not measured, since the tissues were preserved for wide-angle X-ray scattering analysis⁷⁰ after inflation testing. Variation in regional LC thickness may contribute to the regional variation in inflation strain response. Improved depth penetration of the imaging method is needed to investigate the differences across z-depth. 

Conclusions

A total of nine structural features and five strain outcomes were measured in 10 inflation-tested human LC specimens. The structural features measured in 2D sections were found to vary with LC area, region, and to be correlated with strains. Specifically, the main findings include the following: 

1.  
   Of the nine structural features, pore area fraction was the most predictive of maximum principal strain E_max, while connectivity was most predictive of maximum shear strain Γ_max.
2.  
   The peripheral regions had higher pore area fraction, tortuosity, maximum principal strain E_max, maximum shear strain Γ_max, radial strain E_rr, and circumferential strain E_θθ. The peripheral region also had a lower node density and connectivity than the central region.
3.  
   The maximum principal strain E_max, radial strain E_rr, and pore area fraction were larger in the inferior, superior, and temporal quadrants than in the peripheral nasal quadrant.
4.  
   Maximum principal strain E_max, radial strain E_rr, circumferential strain E_θθ, pore area fraction, and beam aspect ratio were larger with larger LC area; connectivity and beam orientation were smaller as with larger LC area.
5.  
   All strain outcomes except for radial circumferential strain E_rθ decreased with age but none of the structural measures decreased significantly with age.
6.  
   The pore area increased as the connectivity and beam orientation decreased and beam aspect ratio increased. A longer beam length was associated with lower node density and a larger beam anisotropy and aspect ratio.

Overall, variations in the LC network structure were associated with variations in strains, which may affect the biomechanical and physiological support for RGC axons. Thus, variations in the LC beam structure may be predictive of the susceptibility and progression of glaucomatous axonal damage. The difference between central and peripheral regions may explain the early vision loss in mid-peripheral regions associated with glaucoma.⁷¹ 

Acknowledgments

The authors thank the eye donors and their families who graciously provided ocular tissue to our project through Eversight. 

Supported by National Institutes of Health Grants EY01765 and EY02120, National Science Foundation CMMI Award 1727104, Brightfocus Foundation G2015132, Public Health Service Research Grants EY021500 and EY001765, Microscopy and Imaging Core Module, Wilmer Core Grant for Vision Research and the Croucher Foundation (YTTL). 

Disclosure: Y.T.T. Ling, None; R. Shi, None; D.E. Midgett, None; J.L. Jefferys, None; H.A. Quigley, None; T.D. Nguyen, None