June 2019
Volume 60, Issue 7
Open Access
Glaucoma  |   June 2019
Characterizing the Collagen Network Structure and Pressure-Induced Strains of the Human Lamina Cribrosa
Author Affiliations & Notes
  • Yik Tung Tracy Ling
    Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, Maryland, United States
  • Ran Shi
    Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, Maryland, United States
    Department of Materials Science and Engineering, Beijing Institute of Technology, Beijing, China
  • Dan E. Midgett
    Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, Maryland, United States
  • Joan L. Jefferys
    Wilmer Ophthalmological Institute, School of Medicine, The Johns Hopkins University, Baltimore, Maryland, United States
  • Harry A. Quigley
    Wilmer Ophthalmological Institute, School of Medicine, The Johns Hopkins University, Baltimore, Maryland, United States
  • Thao D. Nguyen
    Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, Maryland, United States
  • Correspondence: Thao D. Nguyen, Department of Mechanical Engineering, The Johns Hopkins University, Baltimore, MD 21218, USA; vicky.nguyen@jhu.edu
Investigative Ophthalmology & Visual Science June 2019, Vol.60, 2406-2422. doi:10.1167/iovs.18-25863
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Yik Tung Tracy Ling, Ran Shi, Dan E. Midgett, Joan L. Jefferys, Harry A. Quigley, Thao D. Nguyen; Characterizing the Collagen Network Structure and Pressure-Induced Strains of the Human Lamina Cribrosa. Invest. Ophthalmol. Vis. Sci. 2019;60(7):2406-2422. doi: 10.1167/iovs.18-25863.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
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) decreased with higher node density. The maximum principal strain (Emax) 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) than the central regions. The peripheral nasal region had the lowest Emax, Γ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).1 In both human glaucoma and experimental models with increased IOP, RGC axonal injury typically occurs at the LC.2 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.3 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.46 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,912 Thus, there are both cellular and noncellular response elements that are thought to be involved in mechanism of glaucomatous axonal damage.6,1315 The LC structure also undergoes significant connective tissue remodeling in glaucoma.2 This may alter the physiological support for the resident cells and LC capillaries,16 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.17 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,1822 scanning electron microscopy,23 polarized light microscopy,24 optical coherence tomography (OCT),2528 and confocal adaptive optics scanning laser ophthalmoscopy29 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 changes3032 have shown that mechanical strains within the LC exhibit local regions of high shear, tensile, and compressive strains. Midgett et al.32 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,3336 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.32 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.32 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.32 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
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 filter40 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.41 (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 (AC): 200 μm.
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 filter40 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.41 (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 (AC): 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 equalization37 in FIJI38 (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,40 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,43 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 method41 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.,44 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.
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.
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.45 
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.,39 was applied to the preprocessed SHG volumes at 5 mm Hg and 45 mm Hg to determine the components of the 3D displacement field, UX (X, Y, Z), UY (X, Y, Z), and UZ (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 EXX and EYY signify the elongation or compression in the nasal-temporal direction and inferior-superior direction, respectively, and EXY signifies the angle distortion between the X-Y directions. We calculated the out-of-plane strain components, EZZ, but do not report them here because preliminary analysis showed large errors for EZZ (Supplementary Section S3). The in-plane strain components were used to calculate the maximum principal strain, Emax, 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, Err, the circumferential strain, Eθθ, and the radial-circumferential shear strain, E, 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.32 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 UI and strain EIJ 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 UX and UY and under 5.256 ± 3.232 μm (1.90 ± 1.17 pixels) for UZ for all specimens. The absolute strain error for each specimen was below 0.003 for EXX, EYY, and EXY, and below 0.045 for EZZ (Supplementary Table S3).  
Statistical Analysis
The nine structural features; strain components Err, Eθθ, and E; maximum principal strain Emax; 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.47 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 Emax and Γmax. Structural features significantly associated with the average Emax 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. S1S10). The normal strain components Err and Eθθ were predominantly positive (tensile) and concentrated at the boundary with the peripapillary sclera. The shear strains Γmax and E 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, 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 Emax (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 Emax, (B) maximum shear strain Γmax, (C) radial strain Err, (D) circumferential strain Eθθ, and (E) radial-circumferential shear strain E. (F) A box plot of the specimen-averaged strains.
Figure 4
 
The thickness-averaged strain response for inflation from 5 to 45 mm Hg for specimen 6 showing (A) the maximum principal strain Emax, (B) maximum shear strain Γmax, (C) radial strain Err, (D) circumferential strain Eθθ, and (E) radial-circumferential shear strain E. (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 μm2/μm2, 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 R2 Correlation Coefficient
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 R2 Correlation Coefficient
Variations With Age
The structural features did not vary significantly with age (Supplementary Figs. S12, S13). The radial strain Err, circumferential strain Eθθ, maximum tensile strain Emax, and maximum shear strain Γmax all decreased significantly with age, while E, which was an order of magnitude smaller than Err 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) Err radial strain, (B) Eθθ circumferential strain, and (C) E radial-circumferential shear strain.
Figure 5
 
Linear regression analysis of the specimen-averaged strain outcomes for variation with age, showing (A) Err radial strain, (B) Eθθ circumferential strain, and (C) E radial-circumferential shear strain.
Variations With LC Area
The LC area for the 10 specimens ranged from 2.61 to 4.14 mm2. 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, Err, Eθθ, and Emax, were also associated with higher LC area, while the shear strains, Γmax and E, 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 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) Emax maximum principal strain and (B) Err radial strain.
Figure 7
 
Linear regression analysis of the specimen-averaged strain outcomes for variation with LC area, showing (A) Emax maximum principal strain and (B) Err 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 Emax 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 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 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) Emax maximum principal strain, (B) Γmax maximum shear strain, (C) Err radial strain, (D) Eθθ circumferential strain.
Figure 9
 
Effect of radial position on strain outcomes, showing (A) Emax maximum principal strain, (B) Γmax maximum shear strain, (C) Err 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 S13S18). 
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 Emax, maximum shear strain Γmax, and radial-circumferential shear strain E also varied significantly among the quadrants (Supplementary Table S16). For all strain outcomes except E, the strain was highest in the temporal quadrant and lowest in the nasal quadrant (P = 0.03 for Emax 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 Emax (P = 0.02), maximum shear strain Γmax (P = 0.03), and radial-circumferential shear strain E (P = 0.03) also varied among the central quadrants. The Emax (P = 0.02) and Γmax (P = 0.04) were higher in the inferior region than nasal region. Emax 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 Emax (P = 0.001), Γmax (P < 0.001), and Err (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.
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 S19S23). The maximum principal strain Emax 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 R2 values) were between the radial strain Err 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 R2 Correlation Coefficient
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 R2 Correlation Coefficient
Figure 11
 
The specimen-averaged radial strain Err 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).
Figure 11
 
The specimen-averaged radial strain Err 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 S24S28). All strain components increased significantly with increased pore area fraction and all strains (except for E) decreased with higher node density. Moreover, Emax and Γmax decreased with higher connectivity; Emax and Eθθ increased with higher beam aspect ratio; and Γmax increased with higher tortuosity (Table 4). The shear strain E 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
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 Emax 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 Emax 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 Emax 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 Emax increased with larger pore area fraction (n = 80); (B) the maximum shear strain Γmax increased with lower connectivity (n = 79); (C) the radial strain Err (n = 80) and (D) the circumferential strain Eθθ (n = 79) decreased with larger node density.
Figure 12
 
(A) The maximum principal strain Emax increased with larger pore area fraction (n = 80); (B) the maximum shear strain Γmax increased with lower connectivity (n = 79); (C) the radial strain Err (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.48 and 46.7 ± 3.2 μm by Wang et al.28). 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,18 but somewhat higher than the 34.5% ± 7.6% measured by SEM.23 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.49 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 Emax 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.5053 
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, Err and Eθθ, and shear strain E. 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 Emax, Γmax, and Err 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 Emax 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,29 swept-source OCT,28 and adaptive optics spectral-domain OCT.5557 These advanced imaging methods permit high-quality visualization of the LC structure, including beam thickness and pore geometry and depth variations of these features.58 
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,46 accumulation of elastin and collagen with age,59 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.60 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,32 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 EZZ, and EZZ 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.6769 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.18 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 analysis70 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 Emax, while connectivity was most predictive of maximum shear strain Γmax.
  2.  
    The peripheral regions had higher pore area fraction, tortuosity, maximum principal strain Emax, maximum shear strain Γmax, radial strain Err, and circumferential strain Eθθ. The peripheral region also had a lower node density and connectivity than the central region.
  3.  
    The maximum principal strain Emax, radial strain Err, and pore area fraction were larger in the inferior, superior, and temporal quadrants than in the peripheral nasal quadrant.
  4.  
    Maximum principal strain Emax, radial strain Err, 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 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.71 
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 
References
Burgoyne CF. The morphological difference between glaucoma and other optic neuropathies. J neuroophthalmol. 2015; 35: S8–S21.
Quigley HA, Green WR. The histology of human glaucoma cupping and optic nerve damage: clinicopathologic correlation in 21 eyes. Ophthalmology. 1979; 86: 1803–1827.
Albon J, Karwatowski WS, Avery N, Easty DL, Duance VC. Changes in the collagenous matrix of the aging human lamina cribrosa. Br J Ophthalmol. 1995; 79: 368–375.
Hernandez MR, Igoe F, Neufeld AH. Cell culture of the human lamina cribrosa. Invest Ophthalmol Vis Sci. 1988; 29: 78–89.
Hernandez MR. The optic nerve head in glaucoma: role of astrocytes in tissue remodeling. Prog Retin Eye Res. 2000; 19: 297–321.
Wallace DM, O'Brien CJ. The role of lamina cribrosa cells in optic nerve head fibrosis in glaucoma. Exp Eye Res. 2016; 142: 102–109.
Burgoyne CF, Downs JC, Bellezza AJ, Suh JK, Hart RT. The optic nerve head as a biomechanical structure: a new paradigm for understanding the role of IOP-related stress and strain in the pathophysiology of glaucomatous optic nerve head damage. Prog Retin Eye Res. 2005; 24: 39–73.
Tamm ER, Ethier CR, Dowling JE, et al. Biological aspects of axonal damage in glaucoma: a brief review. Exp Eye Res. 2017; 157: 5–12.
Hernandez MR, Pena JD, Selvidge JA, Salvador-Silva M, Yang P. Hydrostatic pressure stimulates synthesis of elastin in cultured optic nerve head astrocytes. Glia. 2000; 32: 122–136.
Morgan JE. Optic nerve head structure in glaucoma: astrocytes as mediators of axonal damage. Eye. 2000; 14: 437–444.
Sun D, Moore S, Jakobs TC. Optic nerve astrocyte reactivity protects function in experimental glaucoma and other nerve injuries. J Exp Med. 2017; 214: 1411–1430.
Clissmann D, Murphy R, Low R, et al. Biomimetic modelling of the glaucomatous lamina cribrosa region of the optic nerve head using tissue engineered scaffolds. 8th World Congr Biomech. 2017; 186: S479–S479.
Varela HJ, Hernandez MR. Astrocyte responses in human optic nerve head with primary open-angle glaucoma. J Glaucoma. 1997; 6: 303–313.
Hernandez MR, Agapova OA, Yang P, Salvador-Silva M, Ricard CS, Aoi S. Differential gene expression in astrocytes from human normal and glaucomatous optic nerve head analyzed by cDNA microarray. Glia. 2002; 38: 45–64.
Tovar-Vidales T, Wordinger RJ, Clark AF. Identification and localization of lamina cribrosa cells in the human optic nerve head. Exp Eye Res. 2016; 147: 94–97.
Anderson DR. Ultrastructure of human and monkey lamina cribrosa and optic nerve head. Arch Ophthalmol. 1969; 82: 800–814.
Quigley HA, Addicks EM. Regional differences in the structure of the lamina cribrosa and their relation to glaucomatous optic nerve damage. Arch Ophthalmol. 1981; 99: 137–143.
Ogden TE, Duggan J, Danley K, Wilcox M, Minckler DS. Morphometry of nerve fiber bundle pores in the optic nerve head of the human. Exp Eye Res. 1988; 46: 559–568.
Yan DB, Coloma FM, Metheetrairut A, Trope GE, Heathcote JG, Ethier CR. Deformation of the lamina cribrosa by elevated intraocular pressure. Br J Ophthalmol. 1994; 78: 643–648.
Downs JC, Yang H, Girkin C, et al. Three-dimensional histomorphometry of the normal and early glaucomatous monkey optic nerve head: neural canal and subarachnoid space architecture. Invest Ophthalmol Vis Sci. 2017; 48: 3195–3208.
Bellezza AJ, Rintalan CJ, Thompson HW, Downs JC, Hart RT, Burgoyne CF. Anterior scleral canal geometry in pressurised (IOP 10) and non-pressurised (IOP 0) normal monkey eyes. Br J Ophthalmol. 2003; 87: 1284–1290.
Sun D, Lye-Barthel M, Masland RH, Jakobs TC. The morphology and spatial arrangement of astrocytes in the optic nerve head of the mouse. J Comp Neurol. 2009; 516: 1–19.
Jonas JB, Mardin CY, Schlötzer-Schrehardt U, Naumann GO. Morphometry of the human lamina cribrosa surface. Invest Ophthalmol Vis Sci. 1991; 32: 401–405.
Jan NJ, Lathrop K, Sigal IA. Collagen architecture of the posterior pole: high-resolution wide field of view visualization and analysis using polarized light microscopy. Invest Ophthalmol Vis Sci. 2017; 58: 735–744.
Inoue R, Hangai M, Kotera Y, et al. Three-dimensional high-speed optical coherence tomography imaging of lamina cribrosa in glaucoma. Ophthalmology. 2009; 116: 214–222.
Strouthidis NG, Fortune B, Yang H, Sigal IA, Burgoyne CF. Effect of acute intraocular pressure elevation on the monkey optic nerve head as detected by spectral domain optical coherence tomography. Invest Ophthalmol Vis Sci. 2011; 52: 9431–9437.
Fatehee N, Paula KY, Morgan WH, Cringle SJ, Yu DY. The impact of acutely elevated intraocular pressure on the porcine optic nerve head. Invest Ophthalmol Vis Sci. 2011; 52: 6192–6198.
Wang B, Nevins JE, Nadler Z, et al. In vivo lamina cribrosa micro-architecture in healthy and glaucomatous eyes as assessed by optical coherence tomography. Invest Ophthalmol Vis Sci. 2013; 54: 8270–8274.
Ivers KM, Li C, Patel N, et al. Reproducibility of measuring lamina cribrosa pore geometry in human and nonhuman primates with in vivo adaptive optics imaging. Invest Ophthalmol Vis Sci. 2011; 52: 5473–5480.
Sigal IA, Grimm JL, Jan NJ, Reid K, Minckler DS, Brown DJ. Eye-specific IOP-induced displacements and deformations of human lamina cribrosa. Invest Ophthalmol Vis Sci. 2014; 55: 1–15.
Girard MJ, Beotra MR, Chin KS, et al. In vivo 3-dimensional strain mapping of the optic nerve head following intraocular pressure lowering by trabeculectomy. Ophthalmology. 2016; 123: 1190–1200.
Midgett DE, Pease ME, Jefferys JL, et al. The pressure-induced deformation response of the human lamina cribrosa: analysis of regional variations. Acta Biomater. 2017; 53: 123–139.
Drance SM. The glaucomatous visual field. Br J Ophthalmol. 1972; 56: 186.
Jonas JB, Fernández MC, Stürmer J. Pattern of glaucomatous neuroretinal rim loss. Ophthalmology. 1993; 100: 63–68.
Jonas JB, Budde WM, Panda-Jonas S. Ophthalmoscopic evaluation of the optic nerve head. Surv Ophthalmol. 1999; 43: 293–320.
Quigley HA. Glaucoma. Lancet. 2011; 377: 1367–1377.
Zuiderveld K. Contrast limited adaptive histogram equalization. In: Graphic Gems IV. San Diego, California: Academic Press Professional; 1994: 474–485.
Schindelin J, Arganda-Carreras I, Frise E, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012; 9: 676–682.
Bar-Kochba E, Toyjanova J, Andrews E, Kim KS, Franck C. A fast iterative digital volume correlation algorithm for large deformations. Exp Mech. 2015; 55: 261–274.
Campbell IC, Coudrillier B, Mensah J, Abel RL, Ethier CR. Automated segmentation of the lamina cribrosa using Frangi's filter: a novel approach for rapid identification of tissue volume fraction and beam orientation in a trabeculated structure in the eye. J R Soc Interface. 2015; 12: 20141009.
Otsu N. A threshold selection method from gray-level histograms. IEEE Trans Syst Man Cybern. 1979; 20: 62–66.
Kroon DJ. Fork of ‘Hessian-based Frangi Vesselness Filter'. 2013.
Antiga L. Generalizing vesselness with respect to dimensionality and shape. Insight J. 2007; 3: 1–14.
D'Amore A, Stella JA, Wagner WR, Sacks MS. Characterization of the complete fiber network topology of planar fibrous tissues and scaffolds. Biomaterials. 2010; 31: 5345–5354.
Berens P. CircStat: a MATLAB toolbox for circular statistics. J Stat Softw. 2009; 31.
Midgett DE, Jefferys JL, Quigley HA, Nguyen TD. The contribution of sulfated glycosaminoglycans to the inflation response of the human optic nerve head. Invest Ophthalmol Vis Sci. 2018; 59: 3144–3154.
Kramer CY. Extension of multiple range tests to group means with unequal numbers of replications. Biometrics. 1956; 12: 307–310.
Nadler Z, Wang B, Schuman JS, et al. In vivo three-dimensional characterization of the healthy human lamina cribrosa with adaptive optics spectral-domain optical coherence tomography. Invest Ophthalmol Vis Sci. 2014; 55: 6459–6466.
Brazile BL, Hua Y, Jan NJ, Wallace J, Gogola A, Sigal IA. Thin lamina cribrosa beams have different collagen microstructure than thick beams. Invest Ophthalmol Vis Sci. 2018; 59: 4653–4661.
Rogers RS, Dharsee M, Ackloo S, Sivak JM, Flanagan JG. Proteomics analysis of human optic nerve head astrocytes following biomechanical strain. Mol Cell Proteomics. 2011; 11: M111.
Balaratnasingam C, Kang MH, Yu P, et al. Comparative quantitative study of astrocytes and capillary distribution in optic nerve laminar regions. Exp Eye Res. 2014; 121: 11–22.
Wang R, Seifert P, Jakobs TC. Astrocytes in the optic nerve head of glaucomatous mice display a characteristic reactive phenotype. Invest Ophthalmol Vis Sci. 2017; 58: 924–932.
Tehrani S, Davis L, Cepurna WO, et al. Astrocyte structural and molecular response to elevated intraocular pressure occurs rapidly and precedes axonal tubulin rearrangement within the optic nerve head in a rat model. PLoS One. 2016; 11: 1–22.
Lockwood H, Reynaud J, Gardiner S, et al. Lamina cribrosa microarchitecture in normal monkey eyes part 1: methods and initial results. Invest Ophthalmol Vis Sci. 2015; 56: 1618–1637.
Hermann B, Fernández EJ, Unterhuber A, et al. Adaptive-optics ultrahigh-resolution optical coherence tomography. Opt Lett. 2004; 29: 2142–2144.
Nadler Z, Wang B, Wollstein G, et al. Automated lamina cribrosa microstructural segmentation in optical coherence tomography scans of healthy and glaucomatous eyes. Biomed Opt Express. 2013; 4: 2596–2608.
Nadler Z, Wang B, Wollstein G, et al. Repeatability of in vivo 3D lamina cribrosa microarchitecture using adaptive optics spectral domain optical coherence tomography. Biomed Opt Express. 2014; 5: 1114–1123.
Duan XJ, Jefferys JL, Quigley HA. Evaluation of automated segmentation algorithms for optic nerve head structures in optical coherence tomography images. Invest Ophthalmol Vis Sci. 2018; 59: 3816–3826.
Albon J, Karwatowski WS, Easty DL, Sims TJ, Duance VC. Age related changes in the non-collagenous components of the extracellular matrix of the human lamina cribrosa. Br J Ophthalmol. 2000; 84: 311–317.
Coudrillier B, Tian J, Alexander S, Myers KM, Quigley HA, Nguyen TD. Biomechanics of the human posterior sclera: age-and glaucoma-related changes measured using inflation testing. Invest Ophthalmol Vis Sci. 2012; 53: 1714–1728.
Fazio MA, Grytz R, Morris JS, Bruno L, Girkin CA, Downs JC. Human scleral structural stiffness increases more rapidly with age in donors of african descent compared to donors of european descent. Invest Ophthalmol Vis Sci. 2014; 55: 7189–7198.
Coudrillier B, Pijanka J, Jefferys JL, et al. Collagen structure and mechanical properties of the human sclera: analysis for the effects of age. J Biomech Eng. 2015; 137: 41006.
Jonas JB, Budde WM. Diagnosis and pathogenesis of glaucomatous optic neuropathy: morphological aspects. Prog Retin Eye Res. 2000; 19: 1–40.
Pijanka JK, Coudrillier B, Ziegler K, et al. Quantitative mapping of collagen fiber orientation in non-glaucoma and glaucoma posterior human sclerae. Invest Ophthalmol Vis Sci. 2012; 53: 5258–5270.
Quigley HA, Varma R, Tielsch JM, Katz J, Sommer A, Gilbert DL. The relationship between optic disc area and open-angle glaucoma: the Baltimore Eye Survey. J Glaucoma. 1999; 8: 347–352.
Healey PR, Mitchell P. The relationship between optic disc area and open-angle glaucoma. J Glaucoma. 2000; 9: 203–204.
Dandona L, Quigley HA, Brown AE, Enger C. Quantitative regional structure of the normal human lamina cribrosa: a racial comparison. Arch Ophthalmol. 1990; 108: 393–398.
Jones HJ, Girard MJ, White N, et al. Quantitative analysis of three-dimensional fibrillar collagen microstructure within the normal, aged and glaucomatous human optic nerve head. J R Soc Interface. 2015; 12: 20150066.
Gogola A, Jan NJ, Lathrop KL, Sigal IA. Radial and circumferential collagen fibers are a feature of the peripapillary sclera of human, monkey, pig, cow, goat, and sheep. Invest Ophthalmol Vis Sci. 2018; 59: 4763–4774.
Pijanka JK, Markov PP, Midgett DE, et al. Quantification of collagen fiber structure using second harmonic generation imaging and two-dimensional discrete fourier transform analysis: application to the human optic nerve head. J Biophotonics. 2018; 12: e201800376.
Brais P, Drance SM. The temporal field in chronic simple glaucoma. Arch Ophthalmol. 1972; 88: 518–522.
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 filter40 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.41 (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 (AC): 200 μm.
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 filter40 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.41 (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 (AC): 200 μm.
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.
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.
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.
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.
Figure 4
 
The thickness-averaged strain response for inflation from 5 to 45 mm Hg for specimen 6 showing (A) the maximum principal strain Emax, (B) maximum shear strain Γmax, (C) radial strain Err, (D) circumferential strain Eθθ, and (E) radial-circumferential shear strain E. (F) A box plot of the specimen-averaged strains.
Figure 4
 
The thickness-averaged strain response for inflation from 5 to 45 mm Hg for specimen 6 showing (A) the maximum principal strain Emax, (B) maximum shear strain Γmax, (C) radial strain Err, (D) circumferential strain Eθθ, and (E) radial-circumferential shear strain E. (F) A box plot of the specimen-averaged strains.
Figure 5
 
Linear regression analysis of the specimen-averaged strain outcomes for variation with age, showing (A) Err radial strain, (B) Eθθ circumferential strain, and (C) E radial-circumferential shear strain.
Figure 5
 
Linear regression analysis of the specimen-averaged strain outcomes for variation with age, showing (A) Err radial strain, (B) Eθθ circumferential strain, and (C) E radial-circumferential shear strain.
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 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) Emax maximum principal strain and (B) Err radial strain.
Figure 7
 
Linear regression analysis of the specimen-averaged strain outcomes for variation with LC area, showing (A) Emax maximum principal strain and (B) Err radial strain.
Figure 8
 
Effects of radial positions on structural features, showing (A) pore area fraction, (B) node density, (C) connectivity, and (D) tortuosity.
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) Emax maximum principal strain, (B) Γmax maximum shear strain, (C) Err radial strain, (D) Eθθ circumferential strain.
Figure 9
 
Effect of radial position on strain outcomes, showing (A) Emax maximum principal strain, (B) Γmax maximum shear strain, (C) Err radial strain, (D) Eθθ circumferential strain.
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.
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.
Figure 11
 
The specimen-averaged radial strain Err 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).
Figure 11
 
The specimen-averaged radial strain Err 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).
Figure 12
 
(A) The maximum principal strain Emax increased with larger pore area fraction (n = 80); (B) the maximum shear strain Γmax increased with lower connectivity (n = 79); (C) the radial strain Err (n = 80) and (D) the circumferential strain Eθθ (n = 79) decreased with larger node density.
Figure 12
 
(A) The maximum principal strain Emax increased with larger pore area fraction (n = 80); (B) the maximum shear strain Γmax increased with lower connectivity (n = 79); (C) the radial strain Err (n = 80) and (D) the circumferential strain Eθθ (n = 79) decreased with larger node density.
Table 1
 
Human Donor Posterior Sclerae Subjected to Inflation Testing
Table 1
 
Human Donor Posterior Sclerae Subjected to Inflation Testing
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 R2 Correlation Coefficient
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 R2 Correlation Coefficient
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 R2 Correlation Coefficient
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 R2 Correlation Coefficient
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
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
Supplement 1
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×