In Vivo Measurements of Prelamina and Lamina Cribrosa Biomechanical Properties in Humans

Abstract

Purpose: To develop and use a custom virtual fields method (VFM) to assess the biomechanical properties of human prelamina and lamina cribrosa (LC) in vivo.

Methods: Clinical data of 20 healthy, 20 ocular hypertensive (OHT), 20 primary open-angle glaucoma, and 16 primary angle-closure glaucoma eyes were analyzed. For each eye, the intraocular pressure (IOP) and optical coherence tomography (OCT) images of the optic nerve head (ONH) were acquired at the normal state and after acute IOP elevation. The IOP-induced deformation of the ONH was obtained from the OCT volumes using a three-dimensional tracking algorithm and fed into the VFM to extract the biomechanical properties of the prelamina and the LC in vivo. Statistical measurements and P values from the Mann-Whitney-Wilcoxon tests were reported.

Results: The average shear moduli of the prelamina and the LC were 64.2 ± 36.1 kPa and 73.1 ± 46.9 kPa, respectively. The shear moduli of the prelamina of healthy subjects were significantly lower than those of the OHT subjects. Comparisons between healthy and glaucoma subjects could not be made robustly due to a small sample size.

Conclusions: We have developed a methodology to assess the biomechanical properties of human ONH tissues in vivo and provide preliminary comparisons in healthy and OHT subjects. Our proposed methodology may be of interest for glaucoma management.

Glaucoma is an ocular disease that is often, but not always, associated with elevated intraocular pressure (IOP). It is characterized by irreversible and progressive damage to the retinal ganglion cell axons within the optic nerve head (ONH). At the ONH, the axons converge and exit the eye through the lamina cribrosa (LC), a porous connective tissue structure. Studies have suggested that elevated IOP (or changes in other loads, such as cerebrospinal fluid pressure (CSFP) and optic nerve traction) deforms the LC, which may lead to structural and physiological changes,^1–4 and potentially cause axonal death.^5,6 However, there is a wide range of patients’ susceptibilities to IOP elevations,^7,8 suggesting that the biomechanics of the ONH is important for understanding the responses of the ONH tissues to IOP elevation. 

ONH biomechanics is the interaction between IOP, the ONH geometry, and the ONH biomechanical properties, the first two of which can be clinically assessed.^9,10 ONH biomechanical properties have been shown to alter with the occurrence of glaucoma, and are regarded as a potential biomarker for glaucoma diagnosis and prognosis. For example, the stiffness of the sclera has been found to increase with age and with exposure to chronic IOP elevation in monkey^11,12 and human^13,14 eyes. However, ONH biomechanical properties have only been measured using ex vivo mechanical experiments, such as uniaxial and biaxial testing of tissue strips,^15–17 inflation tests of eye globes combined with displacement tracking, and inverse finite element methods (IFEM).^13,18–20 

We have recently proposed a technique to extract the biomechanical properties of the ONH tissues in vivo from optical coherence tomography (OCT) measurements, using the virtual fields method (VFM).²¹ VFM is a computational method designed to extract biomechanical properties from full-field displacement measurements.^22,23 VFM was shown to be faster and more robust than the gold-standard IFEM, and the accuracy of using VFM to extract ONH biomechanical properties was also verified.²¹ 

The aim of the current study was to apply this technique to derive the biomechanical properties of the LC and prelamina in healthy, ocular hypertensive (OHT), and glaucoma eyes. 

Methods

Subject Recruitment and Examination

A total of 76 Chinese subjects were recruited at the Singapore National Eye Center. The study was approved by the SingHealth Centralized Institutional Review Board and adhered to the tenets of the Declaration of Helsinki. Written informed consent was obtained from all subjects. 

Included subjects were older than 50 years old with no history of intraocular surgery. Subjects did not have visual field loss related to diabetic retinopathy or any nonglaucomatous optic neuropathy. Among the subjects, 20 were healthy controls, 20 were OHT, 20 had high-tension primary open-angle glaucoma (POAG), and 16 had primary angle-closure glaucoma (PACG). Glaucoma was defined by vertical cup-disc ratio >0.8, narrowing of the neuroretinal rim, visual field defects in standard automated perimetry (glaucoma hemifield test outside normal limits, a cluster of ≥3, nonedge, contiguous points on the pattern deviation plot, not crossing the horizontal meridian with a probability of <5% being present in age-matched normals, and pattern standard deviation <0.05; test results repeatable on two separate occasions with fixation losses <20%, false-negative error <33%, and false-positive error <33%), and IOP measurement higher than 21 mm Hg at least once after diagnosis. OHT was defined by IOP measurement higher than 21 mm Hg with no glaucomatous optic neuropathy or visual field loss. The glaucoma subjects were newly diagnosed or under pharmacologic treatment; the OHT subjects were under observation; the PACG eyes were treated with laser peripheral iridotomy prior to recruitment. 

Both eyes of each subject underwent a series of clinical examinations, including measurement of visual acuity, measurement of refraction (RK-5 autokeratometer, Canon, Tokyo, Japan), slit-lamp biomicroscopy (model BQ-900; Haag-Streit, Köniz, Switzerland), Goldmann applanation tonometry (AT900 D, Haag-Streit, Köniz, Switzerland), dark-room four mirror gonioscopy (Ocular Instruments Inc., Bellevue, WA, USA), standard automated perimetry (SAP, SITA-Standard 24-2 program; Humphrey Field Analyzer II-750i, Carl Zeiss Meditec, Dublin, CA, USA), and IOP measurement with a Tonopen AVIA applanation tonometer (Reichert Inc., Depew, NY, USA). We selected the test eye as the eye that fulfilled the earlier mentioned criteria, and if both eyes were eligible, one eye was randomly selected. 

Acute Elevations of IOP and OCT Imaging

At the beginning of the experiment, each test eye was treated with tropicamide 1% (Alcon, Puurs, Belgium) for pupillary dilation. The baseline IOP was then measured using a Tonopen, and the first OCT volume of the ONH (baseline volume) was acquired at baseline IOP using spectral-domain OCT (Spectralis; Heidelberg Engineering GmbH, Heidelberg, Germany).²⁴ After acquisition of the first OCT scan, an ophthalmodynamometer (spring-loaded indenter) was used to apply gentle indentation on the temporal side of the lower eyelid (perpendicular to the sclera) to increase IOP. With the ophthalmodynamometer maintained in place to hold the IOP constant, the elevated IOP was measured again with a Tonopen, and a second OCT volume (deformed volume) was acquired. The applied force to increase the IOP was consistently 0.64 N (82.5 g) as calibrated using a uniaxial tensile tester (Instron-5848; Instron Inc., Norwood, MA, USA) and aimed to an IOP value of >35 mm Hg (a preliminary assessment performed on 20 healthy eyes was 34.35 ± 5.21 mm Hg).²⁴ There was no difference in the IOP elevation values (19.1± 5.3 mm Hg) across different diagnostic groups. 

For each test eye, two OCT volumes and two IOP measurements were acquired and used to extract the biomechanical properties of the LC and of the prelamina. Each OCT volume consisted of 97 serial horizontal B-scans (axial resolution 3.8 µm; transverse resolution approximately 12 µm; approximately 30 µm distance between B-scans; 384 A-scans per B-scan; 20X B-scan averaging) that covered a rectangular area of 15° x 10° centered on the ONH. The eye tracking and enhanced depth imaging modalities of the Spectralis were used during image acquisition. Adaptive compensation (postprocessing) was also used to improve ONH tissue visibility at high depth and below blood vessels, and increase contrast to aid the segmentation of the tissues.^25,26 

Geometry Reconstruction and Displacement Tracking

To reconstruct the geometry of the LC, the prelamina and the Bruch's membrane opening (BMO), each OCT volume was manually segmented (Figs. 1a, 1b) using Amira (version 5.6, FEI, Hillsboro, OR, USA) following previously established segmentation protocols.²⁵ A custom MATLAB algorithm was written to mesh the region of interest using eight-node hexahedron elements (Fig. 1c). For each test eye, the deformed OCT volume was oriented to best-align with the baseline OCT volume in Amira, to reduce rigid body transformations. These can occur as a result of head and eye movements between different OCT scan sessions. A three-dimensional (3D) tracking algorithm^25,26 was then applied to each set of OCT volumes (baseline volume and deformed volume) to track tissue displacements induced by acute IOP elevation. Displacements of regions in the prelamina and LC shadowed by the blood vessels in the OCT images were filtered out and interpolated from neighboring regions. 

Figure 1.

 

(a) Manual segmentation of the prelamina and LC using the OCT images. (b) Reconstruction of the geometries of the prelamina and LC in 3D. (c) Meshing the prelamina and LC using eight-node hex elements.

Constitutive Models

The LC and the prelamina were modeled as neo-Hookean hyperelastic materials. The Cauchy stress tensor for neo-Hookean material is defined as:  

\begin{eqnarray}{{\sigma }} = \frac{{{\mu }}}{{{J}}}\left( {{{\bf B}} - {{\bf I}}} \right) + {\rm{p}}{{\bf I}}\end{eqnarray}

(1)

where μ is the shear modulus, J is the determinant of the deformation gradient tensor, B is the left Cauchy-Green deformation tensor, I is the identity matrix, and p is the hydrostatic pressure to ensure quasi-incompressibility. The shear moduli of the prelamina and the LC were unknown and to be extracted. 

Extraction of the ONH Tissue Biomechanical Properties

The biomechanical properties of the ONH tissues were extracted in vivo from OCT measurements using VFM. VFM is based on the principle of virtual work, which can be written as:  

\begin{eqnarray}\mathop \int \nolimits_{{V}}^{} {{\sigma }}\left( {{{\mu }},{{\bf u}}} \right):{{{\varepsilon }}^{{*}}}{{dV}} = \mathop \int \nolimits_{{S}}^{} {{\bf T}} \cdot {{{\bf u}}^{{*}}}{{dS}}\ \ \ \ \ \ \ \ \end{eqnarray}

(2)

where the Cauchy stress σ is a function of unknown constitutive parameters (the shear modulus μ) and the experimentally measured displacement field u; V denotes the tissue volume; T is the external force (here the IOP); S denotes the tissue boundary surface; u* is the virtual displacement, and ε* is the corresponding virtual strain. The latter can be expressed as:  

\begin{eqnarray}{{{\varepsilon }}^{{*}}} = \frac{1}{2}\left( {\nabla {{{\bf u}}^{\rm{*}}} + {\nabla ^{{T}}}{{{\bf u}}^{\rm{*}}}} \right)\end{eqnarray}

(3)

A cost function φ(IOP, u*,  u) was constructed based on the discretization of Equation 2 (see Appendix for derivation). The cost function required the IOP measurement, the displacement field u from the tracking and constructed virtual field u* (in Appendix) as inputs, and the only unknowns were the shear moduli of ONH tissues. The cost function was minimized iteratively using differential evolution, which is a genetic global optimization algorithm.²⁷ The optimization was regarded as converged when the variance of the cost function became 0 for at least 50 consecutive generations (from differential evolution). 

Statistical Analysis

Mean and standard deviations of the extracted shear moduli of the prelamina and the LC were calculated. Effect size (Cohen’s d) was used to calculate the statistical power between diagnosis group (Healthy, POAG, PACG, OHT) comparisons for age, IOP, and extracted shear moduli. For those comparisons whose effect size was larger than 0.5 (considered as medium effect size), a nonparametric statistical hypothesis test, that is, the Mann-Whitney-Wilcoxon test, was performed in R (version 3.0.2, R Foundation, Vienna, Austria). A P value <0.05 was considered statistically significant. 

Results

Subjects with poor LC visibility in OCT volumes (<60% of the total BMO area visible) were excluded from the study (Healthy: 3, OHT: 3, PACG: 1, POAG: 1). Subjects with central prelaminar thickness less than 100 µm (100 µm is the thickness assuring that at least two tracked points exist in the axial direction) were also excluded as the extraction required complete displacement field of the prelamina (Healthy: 2, OHT: 4, PACG: 2, POAG: 7). Finally, a total of 53 subjects were included (Healthy: 15, OHT: 13, PACG: 13, POAG: 12). 

PACG subjects were significantly older (68.4 ± 4.1 years) than healthy (57.2 ± 3.7 years), OHT (60.8 ± 4.6 years), and POAG (61.5 ± 6.4 years) subjects (P < 0.01 using the Mann-Whitney-Wilcoxon tests; Fig. 2a). The baseline IOP in the OHT subjects (20.7 ± 2.6 mm Hg) was significantly higher than that in healthy (17.1 ± 2.6 mm Hg) and PACG (16.9 ± 4.1 mm Hg) subjects (P < 0.05; Fig. 2b). The baseline IOP in POAG subjects was 18.9 ± 3.1 mm Hg. The effect sizes of comparisons for age and baseline IOP are shown in Table A1. 

Figure 2.

 

Box plots of the age and baseline IOP values across diagnostic groups (Healthy, OHT, PACG, and POAG). P values derived from the Mann-Whitney-Wilcoxon tests are shown. (a) The PACG subjects were significantly older than subjects from other diagnostic groups. (b) The baseline IOP values of OHT subjects were significantly higher than healthy and PACG subjects.

The extracted shear moduli of the prelamina were 56.0 ± 36.6 kPa for healthy subjects, 82.4 ± 40.1 kPa for OHT subjects, 57.5 ± 25.2 kPa for PACG subjects, and 62.3 ± 38.1 kPa for POAG subjects (Table 1). The shear moduli of the prelamina of healthy subjects were significantly lower than those of the OHT subjects (P = 0.019; Fig. 3a). 

Table 1.

 

Means and Standard Deviations of Extracted Shear Moduli of Prelamina and LC of Different Diagnostic Groups and all Subjects

 

Figure 3.

 

Box plots of the shear moduli of the prelamina and the LC across diagnostic groups (Healthy, OHT, PACG, and POAG). P values derived from the Mann-Whitney-Wilcoxon tests are shown. (a) The shear moduli of the prelamina of healthy subjects were significantly lower than those of the OHT subjects. (b) There was no statistically significant difference in the extracted shear moduli of the LC across diagnostic groups.

The extracted shear moduli of the LC were 57.4 ± 42.8 kPa for healthy subjects, 79.2 ± 45.6 kPa for OHT subjects, 80.2 ± 50.5 kPa for PACG subjects, and 78.1 ± 50.4 for POAG subjects (Table 1; Fig. 3b). 

Note that the effect sizes for comparing shear moduli (LC and prelamina) across groups are shown in Table 2. Due to small effect sizes (<50%), comparisons between healthy and glaucoma subjects could not be made robustly. 

Table 2.

 

Effect Sizes (Cohen’s d) for Comparisons between Diagnostic Groups

 

For all included subjects, the average shear modulus of prelamina was 64.2 kPa, whereas the average shear modulus of LC was 73.1 kPa (Table 1). 

Discussion

In this study, we used the VFM to extract the biomechanical properties of human ONH tissues in vivo, based on full field displacement data from clinical OCT images. Whenever feasible, we compared and reported the extracted shear moduli of the prelamina and LC across different diagnostic groups. 

The prelaminar tissue is commonly regarded as an extension of the retina into the scleral canal opening. In the literature, the retina is considered to be soft, even though few studies have investigated the biomechanical properties of retina. Chen et al.²⁸ performed tensile testing of human retina strips and reported Young's modulus of the retina to be 10 to 20 kPa at 0.02 strain level. Wollensak and Spoerl²⁹ did similar experiments with porcine retina and reported Young's modulus of retina to be 100 kPa. Worthington et al.³⁰ tested the compression resistance of porcine retina and reported the modulus to be approximately 10 kPa. The accuracy of these ex vivo tensile studies could be compromised by the nonphysiological boundary conditions. Attaching the delicate retina to the experimental device was also reported to be difficult.²⁸ These factors could be the reason behind the wide range of reported retinal stiffness values. Qian et al.³¹ extracted the cat's retinal stiffness in vivo by combining OCT technique with IFEM. They reported a shear modulus of 75.9 kPa, which is similar to our results (64.2 ± 36.1 kPa). 

We found that the modulus of the prelamina was on average smaller than that of the LC, but surprisingly, both were of the same order of magnitude. This was consistently observed across all eyes. Because the prelaminar tissue is mostly composed of cells, but the lamina contains a significant amount of collagen and elastin, we should have expected a much lower modulus for the prelaminar tissue.³² However, it is important to note that the central retinal vessel trunk was also segmented and considered as part of the prelaminar tissue. It is highly plausible that the central retinal vessel trunk (also containing collagen and pressurized owing to the presence of blood) could provide resistance to neural tissue compression, so that the prelaminar tissue stiffness would appear higher than what it actually is. Segmenting the central retinal vessel trunk as a separate class could be considered in future analyses, although we may require an OCT device with a higher resolution to fully map both the stiffness of the prelaminar tissue and central retinal trunk. 

There are few direct experimental measurements of LC stiffness. Spoerl et al. estimated the Young's modulus of the LC to be approximately 100 kPa using tissue strips (Spoerl et al. Biomechanical behaviour of the lamina cribrosa influenced by various substances. Invest. Ophthalmol. Vis. Sci. 2003;44(13):3318. doi:https://doi.org/). Sigal et al.³² calculated the Young's modulus of the LC to be 77 to 405 kPa by fitting an LC model to the monkey ex vivo LC displacement measurements. Edwards and Good³³ compared their mathematical model with the response of LC to elevated IOP, and estimated the Young's modulus of the LC to be 140 to 380 kPa. Our extracted shear moduli of the LC (73.1 ± 46.9 kPa) are consistent with these reported LC stiffness values. 

Our results showed that the OHT subjects had relatively stiffer prelamina compared with the healthy subjects. A stiffer prelamina might provide the optic nerve fibers with more resistance to IOP elevation, which is consistent with the biomechanical theory of glaucoma.⁶ The retina/prelaminar tissues were reported to exhibit slightly nonlinear behavior with low transition strain²⁸ or linear behavior.²⁹ As a result, higher baseline IOP values of OHT subjects might contribute to a stiffer prelamina. Furthermore, it has also been reported that there is an association between OHT and high blood pressure.³⁴ It is possible that higher blood pressure and associated stiffening of blood vessels³⁵ would contribute to a higher compressive resistance. 

Our results showed that healthy subjects had on average more compliant LCs compared with OHT subjects and glaucoma subjects, but we could not report on whether this was statistically significant due to a small sample size. Studies have reported that connective tissues (mainly sclera) in glaucoma eyes were stiffer than those in healthy eyes,^12,18,36 but not for the prelamina and the LC. More work will be required to confirm these trends if they exist. 

Our study has several limitations. First, the viscoelasticity of the ONH tissues was not taken into consideration. In the experiment, IOP was only elevated for 2 to 3 minutes, which might not be long enough for the ONH to fully deform in some subjects. As a result, the extracted stiffness parameters might be overestimated for some subjects, depending on the respective viscoelastic properties. However, the experiment protocol was kept consistent throughout this study. Second, the prelamina and the LC were assumed to be homogeneous, thus potential local variations in the biomechanical properties could not be revealed. This limitation notwithstanding, the assumption of homogeneity reduced the number of unknowns, making the extracted result concise as a potential biomarker for glaucoma. Third, the significance of this study was limited by the small group sizes of the subjects. Low LC visibility and severe prelaminar atrophy in some subjects further reduced the sample size. Hence a large-scale clinical study is warranted to further investigate the intergroup differences. Fourth, the glaucoma subjects in this study had received IOP reduction treatment. The pharmacologic influence on the biomechanical properties of the ONH tissues was not studied. Fifth, only measured axial displacements were used in the extraction of ONH biomechanical properties, and the virtual displacement field used in this study was mainly defined in the axial direction (Equation A3 in Appendix). We decided to only use the measured axial displacements, as the OCT axial resolution is much higher than that in transverse directions. As a consequence, the displacement tracking errors in the radial and the circumferential directions are much higher than that in the axial direction. We found that including the radial and the circumferential displacement components as part of the VFM could compromise the accuracy of the extracted parameters, and we only focused on the axial displacement component. Sixth, the stress-free geometry of the ONH was not available for clinical experiments, and the ONH geometry at normal IOP state was used instead. Future studies considering prestress in the tissues are warranted. Seventh, the current OCT penetration depth is often not good enough to reveal the complete morphology of the LC in every subject. For those LCs in which the posterior surface was not visible, the posterior-most visible surface was treated as the LC posterior surface, and the extracted stiffness only represented the stiffness of that visible portion of the LC. Eighth, the estimated time for IOP measurement and OCT scan was 2 minutes. We applied a constant external force (0.64 N) on every study eye. The IOP may reduce during the OCT scanning due to a change in aqueous humor dynamics. However, the scanning time was short (on average 2 minutes) and the force was applied on the sclera (through the eyelid) rather than on the cornea. Thus we believe that the IOP was approximately constant for individual eyes throughout the OCT scans. 

Conclusions

We extracted the biomechanical properties of human ONH tissues in vivo using a VFM. This study further highlights the potential of in vivo characterization of ONH biomechanics for use as a biomarker for glaucoma, which needs to be explored in a larger patient cohort. 

Future Research

The accuracy of any given inverse method can only be as good as that of the input displacement data. The accuracy of displacement tracking methods is currently limited by the imaging capability of OCT devices. Larger penetration depth, higher lateral resolutions, and faster B-scan acquisition may be essential to provide accurate 3D full-field displacement measurements and subsequent in vivo stiffness extractions. These advancements are likely to come in the future and would also help reduce the variations (owing to error) in extracted stiffness parameters. A large cohort of subjects would also be needed to confirm our results. With more recruited subjects, associations between ONH stiffness, diagnosis, age, and IOP could be investigated using standard statistical methods. 

Acknowledgments

Supported by Singapore Ministry of Education Academic Research Funds Tier 1 (R-397-000-294-114 [MJAG]); Singapore Ministry of Education Academic Research Funds Tier 2 (R-397-000-280-112 and R-397-000-308-112 [MJAG]); National Medical Research Council (Grant NMRC/STAR/0023/2014 [TA]); Cardiff University International Research Leave Fellowship (CB); National Institute for Health Research, Biomedical Research Centre, Moorfields Eye Hospital National Health Service, Foundation Trust and University College London, Institute of Ophthalmology (NGS). The sponsor or funding organization had no role in the design or conduct of this research. The views expressed are those of the author(s) and not necessarily those of the National Health Service, the National Institute for Health Research, or the United Kingdom Department of Health. 

Disclosure: L. Zhang, None; M.R. Beotra, None; M. Baskaran, None; T.A. Tun, None; X. Wang, None; S.A. Perera, None; N.G. Strouthidis, None; T. Aung, None; C. Boote, None; M.J.A. Girard, None 

Appendix

Because the LC and BMO are circular shaped, u* was defined in the cylindrical coordinate system (r, ω, h):  

\begin{eqnarray}\left\{ {\begin{array}{@{}l@{\quad}l@{}} r = \sqrt {{{{x}}^2} + {{{z}}^2}} & r \in \left[ {0,{{R}}} \right]\\ \omega = atan2\left( {{{z}},{{x}}} \right) & \omega \in \left[ {0,\ 2{{\pi }}} \right)\\ h = y & h \in \left[ {{{{H}}_{{{\min}}}},{{{H}}_{{{\max}}}}} \right] \end{array}} \right.\end{eqnarray}

(A1)

where R, H_min, and H_max are the boundary values for nodes located in the region of interest (as shown in Fig. A1a). In this case, H_min is y-coordinate of the center point on the posterior-most visible LC surface; H_max is the y-coordinate of the anterior-most point on the anterior prelaminar surface, as illustrated in Figure A1. atan2 is the arctangent function with two arguments:  

\begin{eqnarray} &&{{atan}}2\left( {{{y}},{{x}}} \right) \nonumber \\ &&\quad = \left\{ \begin{array}{@{}l@{\quad}l@{\;}l@{}} {{arctan}}\left( {\frac{{{y}}}{{{x}}}} \right) + {{\pi }} & {{if}} & {{x}} > 0, \\ {{arctan}}\left( {\frac{{{y}}}{{{x}}}} \right) + 2{{\pi }} & {{if}} & {{x}} < 0\ {{and}}\ {{y}} \ge 0, \\ {{{arctan}}\left( {\frac{{{y}}}{{{x}}}} \right)}& {{if}} & {{x}} < 0\ {{and}}\ {{y}} < 0, \\ {\frac{{3{{\pi }}}}{2}}& {{if}} & {{x}} = 0\ {{and}}\ {{y}} > 0, \\ {\frac{{{\pi }}}{2}\ }& {{if}} & {{x}} = 0\ {{and}}\ {{y}} < 0, \\ {undefined} & {if} & {x} = 0\ {and}\ {y} = 0. \end{array}\right.\nonumber \\ \end{eqnarray}

(A2)

 

Figure A1.

 

(a) Schematic of a human ONH and the definitions of geometric parameters in the constructed virtual field. (b) Illustrations of window function values used to extract the shear modulus of the prelamina (left plot) and to extract the shear modulus of the LC (right plot).

The virtual displacements used to extract the stiffness parameters were expressed as:  

\begin{equation*}{{{\bf u}}^{\rm{*}}} = \left( {{{u}}_{{\rm r}}^{\rm{*}}{{{\bf e}}_{{\bf r}}} + {{ru}}_{{\rm w}}^*{{{\bf e}}_{{\bf w}}} + {{u}}_{{\rm h}}^{\rm{*}}{{{\bf e}}_{{\bf h}}}} \right){\cal L}\end{equation*}

 

\begin{eqnarray} {{u}}_{{\rm r}}^* = 0;\;u_{\rm w}^* = 0;\;u_{\rm h}^* = - \frac{{\sin \left( {\frac{{\pi (R - r)}}{{2R}}} \right)(h - {H_{\min }})}}{{{H_{\max }} - {H_{\min }}}} \nonumber \\ \end{eqnarray}

(A3)

where \({\cal L}\) is a scalar window function to enforce the required boundary conditions for the virtual fields. The window functions were defined as unity for most of the targeted tissues in each step, and gradually decreased to zero at required regions and surfaces (as shown in Fig. A1b). 

The distribution of T may not be known at every part of the ONH. For instance, for the eye, the only measurable external load is IOP; the forces exerted on the posterior surfaces of the peripapillary sclera and LC, for example, from the CSFP or orbital fat, and the hoop stresses acting on the ONH boundaries are typically unknown. Fortunately, these unknown boundary forces can be excluded from the calculation of the external virtual work by applying zero virtual surface displacement at the corresponding locations. This is explained using the following equation: 

\begin{eqnarray} \int_{{{V}}}^{\;}{{{{\sigma }}:{{{\varepsilon }}^*}{{dV}}}} &=& \int_{{{S}}}^{\;}{{{{\bf T}} \cdot {{{\bf u}}^*}{{dS}} = {{IOP}}}} \int_{{{{{S}}_{{{\rm IOP}}}}}}^{\;}{{{{\bf n}} \cdot {{{\bf u}}^*}_{{{\rm IOP}}}{{dS}}}} \nonumber \\ && +\, \int_{{{{{S}}_{{\bf O}}}}}^{\;}{{{{{\bf T}}_{{\bf O}}} \cdot {{\bf u}}_{{ O}}^*{{dS}}}}\end{eqnarray}

(A4)

where S_IOP is the inner surface of the prelamina to which IOP is applied (Fig. A1a), S_O represents the surface of the ONH other than S_IOP (including the posterior surface of the visible portion of the LC, and the boundary between the area of interest, and all other tissues as illustrated in Fig. A1), and S = S_IOP∪S_O; n is the unit surface normal vector, T_O represents the unknown boundary forces on S_O, u*_IOP is the virtual displacement at S_IOP, and u*_O is the virtual displacement at S_O. If we define u*_O = 0 over S_O, the second right-hand term of Equation A4 vanishes and T_O does not contribute to the calculation of the virtual work, as shown here:  

\begin{eqnarray}\int_{{{V}}}^{\;}{{{{\sigma }}:{{{\varepsilon }}^*}{{dV}} = {{IOP}}}}\int_{{{{{S}}_{{{\rm IOP}}}}}}^{\;}{{{{\bf n}} \cdot {{\bf u}}_{{{\rm IOP}}}^*{{dS}}}}\end{eqnarray}

(A5)

 

As a result, IOP is the only external force that needs to be known to extract the constitutive parameters. Note that u*_IOP should not be zero, otherwise the external virtual work would be zero, in which case only the relative relationship between the constitutive parameters of different tissues can be derived. To be more rigorous, for those LCs whose posterior surface was defined as the posterior-most visible surface, the traction acting on S_O is not CSFP, but the internal traction acting on that surface. Nonetheless, Equations A4 and A5 still hold true in this case. 

A cost function (Equation A6) was constructed based on the discretization of Equation A5:  

\begin{eqnarray}{{\phi }}({{IOP}},{{\bf u}^*},{\bf u}) &=& \bigg|\sum\nolimits_{{{i}} = 1}^{{{\textit{N}}_{\bf e}}} {{{\sigma }}_{{i}}}({\bf \mu },{\bf u}):{{\varepsilon }}_{{i}}^*({{\bf u}^*}){{{V}}_{{i}}} \nonumber \\ && - {{IOP}}\sum\nolimits_{{{j}} = 1}^{{{{N}}_{{\rm s}}}} {{{\bf n}_{{j}}} \cdot {\bf u}_{{j}}^*{{S}_{{j}}}} \bigg|\end{eqnarray}

(A6)

where N_e is the number of elements within the region of interest, N_s is the number of surface elements exposed to IOP, n_j is the unit normal surface vector of surface element j. Two-dimensional four-point Gaussian quadrature and 3D eight-point Gaussian quadrature were used to calculate the surface integral of external virtual work and the volume integral of internal virtual work respectively. 

Table A1.

 

Effect Sizes (Cohen’s d) for Comparisons Between Diagnostic Groups

 

A two-step procedure was used here: first, the shear modulus of the prelamina was extracted by applying zero virtual displacements in the LC region (enforced by applying zero window function values in the LC, as shown in Fig. A1b); second, the shear modulus of the LC was extracted, with the shear modulus of the prelamina known. In each step, only one unknown parameter existed, which reduced the computational cost. The accuracy of using the VFM to extract the biomechanical properties of the ONH tissues in a multistep manner has been numerically verified in our previous article.²¹ 

For simplicity, all tissues were assumed to be incompressible. Therefore the Jacobian J was set to 1 for all elements, thus eliminating the need to include the volumetric stress in the virtual work calculation.