Cerebrospinal Fluid Pressure: Revisiting Factors Influencing Optic Nerve Head Biomechanics

Abstract

Purpose: To model the sensitivity of the optic nerve head (ONH) biomechanical environment to acute variations in IOP, cerebrospinal fluid pressure (CSFP), and central retinal artery blood pressure (BP).

Methods: We extended a previously published numerical model of the ONH to include 24 factors representing tissue anatomy and mechanical properties, all three pressures, and constraints on the optic nerve (CON). A total of 8340 models were studied to predict factor influences on 98 responses in a two-step process: a fractional factorial screening analysis to identify the 16 most influential factors, followed by a response surface methodology to predict factor effects in detail.

Results: The six most influential factors were, in order: IOP, CON, moduli of the sclera, lamina cribrosa (LC) and dura, and CSFP. IOP and CSFP affected different aspects of ONH biomechanics. The strongest influence of CSFP, more than twice that of IOP, was on the rotation of the peripapillary sclera. CSFP had similar influence on LC stretch and compression to moduli of sclera and LC. On some ONHs, CSFP caused large retrolamina deformations and subarachnoid expansion. CON had a strong influence on LC displacement. BP overall influence was 633 times smaller than that of IOP.

Conclusions: Models predict that IOP and CSFP are the top and sixth most influential factors on ONH biomechanics. Different IOP and CSFP effects suggest that translaminar pressure difference may not be a good parameter to predict biomechanics-related glaucomatous neuropathy. CON may drastically affect the responses relating to gross ONH geometry and should be determined experimentally.

Glaucoma, the second leading cause of blindness worldwide,¹ is characterized by a particular pattern of irreversible damage to the retinal ganglion cell axons.² The damage is believed to initiate within the optic nerve head (ONH), where the axons pass through the lamina cribrosa (LC) and exit the eye.^3,4 IOP is considered to be the most important modifiable risk factor for the optic neuropathy of glaucoma, regardless of the level of IOP at which the neuropathy occurs.^5,6 Although the mechanism by which elevated IOP contributes to axon damage remains unclear, it is often considered that IOP affects the susceptibility to glaucoma by causing an altered biomechanical environment within the ONH. Within this framework, the forces from IOP distort the tissues of the ONH, and LC within, triggering events such as compromised axoplasmic flow, vascular perfusion, and astrocyte activation that eventually lead to glaucomatous optic neuropathy.^7–9 However, the ONH is exposed not only to IOP from within the globe, but also to cerebrospinal fluid pressure (CSFP) within the subarachnoid space and blood pressure (BP) within the central retinal vessels. CSFP and BP can also potentially influence the biomechanical environment within the ONH,^10–13 and thus contribute to determine the susceptibility to glaucoma. In fact, in recent years, evidence has been mounting that the susceptibility to glaucoma may be influenced by CSFP.^14–17 Evidence has also been presented that BP may influence the susceptibility to glaucoma.^18,19 

Whilst many studies have addressed the effects of IOP on the ONH using experimental,^20–22 theoretical,^23–25 and numerical^26–29 methods, the potential effects of CSFP or BP on the ONH biomechanical environment have not been studied in nearly as much detail. A better understanding of the effects of these pressures and the potential interactions between their effects is necessary to understand the etiology of glaucoma and the puzzling range of sensitivities to IOP. 

Our goal was to model the sensitivity of the ONH biomechanical environment to acute variations in IOP, CSFP, and BP. To do so, we extended a previously published numerical model of the ONH³⁰ to include a central retinal vessel and more detailed retrolaminar anatomy. We endeavored for this to be the most comprehensive analysis to date of the factors influencing ONH biomechanics. Hence, we studied the effects of 24 factors, including the three pressures and 21 other factors representing ONH and globe geometry and tissue mechanical properties, and the constraints on the optic nerve (CON). We simulated ONH biomechanics using a broad set of 98 responses, including pressure-induced local deformations (strains) and forces (stresses) as well as gross deformations (e.g., peripapillary sclera [PPS] rotation or bowing). 

Methods

The general strategy was to produce a large set of models representing a diversity of ONHs with varying tissue anatomy and mechanical properties and optic nerve constraints, and then use finite element (FE) modeling to simulate the effects on each of the models of acute changes in IOP, CSFP, and/or BP. Due to the large number of factors and responses of interest, we split the analysis into two phases. In a first screening phase, we predicted the main effects and interactions of 24 factors on 98 responses. We identified the 16 most influential factors and 10 responses representative of the whole response set. The representative responses were selected in a process informed by dimensionality reduction techniques and principles of mechanobiology. In a second phase, we focused on the most influential factors from phase 1 to predict factor effects in detail. We then inspected the results using archetypal analysis to identify ONHs representative of the diversity of potential ONH biomechanical responses to the pressures. The steps are described in detail below. 

Modeling

We extended our previous numerical model of the ONH³⁰ to include a central retinal vessel and more detailed retrolaminar anatomy (Fig. 1). The model was then parameterized to allow independent and simultaneous variations in 24 factors. The factors and the ranges over which they were varied are listed in Table 1. For each factor, the range of admissible values was defined from the literature, when available, or from our own estimates, based on measurements on serial sections of the ONH from ostensibly healthy donor human eyes.³¹ The development, processing, simulation, and postprocessing of the FE models were as described elsewhere.^30,32–34 The rationale for using our previous numerical model as the basis for this study, rather than developing a completely new model, and its implications are addressed later in the Discussion. 

Figure 1

 

Input factor definitions. Only the optic nerve head region is shown. See Table 1 for input factor ranges. In addition to the input factors shown, the compressibility (Poisson's ratio) of the prelaminar and retrolaminar neural tissues and the stiffness of each tissue region were varied, for a total of 24 input factors. The illustration represents the model geometry at IOP of 5 mm Hg, BP of 50 mm Hg, and CSFP of 0 mm Hg.

Table 1

 

Factors and Their Ranges

 

Factors

Geometry

Three major changes were made to the model geometry from our previous model.³⁰ First, the dura mater around the optic nerve was included. The space between the dura and pia mater defined the subarachnoid space, holding the cerebrospinal fluid.³⁵ The dura mater was 0.75 mm thick at the junction with the sclera, tapering down to 0.375 mm at 3.0 mm posterior to the junction, remaining constant to the end of the model. The arachnoid was not considered independently in the models, under the assumption that it is extremely thin and soft.³⁶ 

Second, a central retinal vessel was incorporated to model the effects of BP, similar to model 2 in a previous study.³⁴ A single vessel is a simplification of the complex vasculature passing through the ONH. This simplification was intended to capture the main elements of the arteries and veins, modeling them as pressurized vessels with thick collagen-rich walls. The vessel was idealized as a circular cylinder with parameterized external diameter and wall thickness (see below). The vessel wall was assumed to be mechanically attached to the neural tissues of the ONH and the LC, such that they necessarily deformed together. 

Finally, a longer optic nerve was modeled. The optic nerve, including the retrolaminar neural tissue, the surrounding dura and pia mater and the central retinal vessel, was extended posteriorly 10 mm from the sclera. We modeled the optic nerve straight without considering any kinking or waviness in the nerve itself, as their effects would fall somewhere within the CON we considered (fully free and fully constrained, discussed below).³⁷ We also did not account for the retinal vessels' exit from the nerve. 

The model geometry was parameterized using 11 factors representing the ONH anatomy. Four of them were from the above-mentioned changes, which were the subarachnoid space width (the distance between the dura and pia mater), the dura thickness, and the vessel external diameter and wall thickness. Seven other factors were selected from the 16 considered in our previous models based on a preliminary multivariate sensitivity analysis. These factors were the eye radius, the depth, radius, and thickness of the LC, the pia mater thickness, the scleral thickness at the canal, and the scleral shell thickness. The other factors were not varied and were set at the baseline levels used in our previous work.³⁰ The geometric factors are illustrated in Figure 1 and their ranges listed in Table 1. 

Mechanical Properties

To maximize the comparability of this study with our previous ones on factors influencing IOP-related ONH biomechanics, tissues were modeled using the same tissue mechanical properties.³⁰ Tissues were assumed homogeneous, isotropic, and linearly elastic. Tissue stiffness was defined by Young's modulus and compressibility by Poisson's ratio. All tissues other than the neural tissues were modeled as incompressible. Here we provide details only on the range definitions of the new tissues in this work. Details for other materials are provided elsewhere.³⁰ There is little information about the stiffness of dura mater. Based on previous measurements of the Young's modulus of human spinal cord,³⁸ the range of dura modulus was assumed as 1 to 5 MPa. To the best of our knowledge, there are no direct measurements of the central retinal vessel mechanical properties. Based on the Young's modulus of similar-sized vessels in other regions of the human body as well as the strengthening effects of connective tissues,³⁸ the range of vessel stiffness was assumed as 0.5 to 5 MPa. The material factors and their ranges are summarized in Table 1. We recognize that our choices in material properties will have important implications on the results and conclusions. These will be addressed in detail in the Discussion. 

Pressures

Pressures were modeled as a distributed load acting perpendicular on the tissue surface exposed to the pressurized fluid. IOP was applied to the anterior surface of the prelaminar neural tissue, CSFP to the subarachnoid space (the inner surface of the dura mater and the external surface of the pia mater), and BP to the inner luminal surface of the central retinal vessel. The base model was defined to represent a case with low IOP (5 mm Hg) and BP (50 mm Hg), and no CSFP (0 mm Hg). These pressures were then parameterized to represent increases of 15 mm Hg in IOP and CSFP and 30 mm Hg in BP. The range selected for IOP and CSFP corresponds to pressures going from a very low, borderline physiological level, to high levels, elevated but still normal, whereas that for BP represents normal arterial blood pulsation.^39–41 After the analysis on IOP, CSFP, and BP, we also tested the overall influence of the translaminar pressure difference (TLPD = IOP − CSFP). 

Optic Nerve Constraints

The degree to which the optic nerve is constrained is unclear. Elsewhere, modeling studies have been conducted both with a completely free optic nerve boundary^10,11,30,32 and with an optic nerve fully constrained at the end.⁴² Preliminary tests suggested that the choice could have important effects on the predicted ONH biomechanics. Hence, to avoid biasing this study we decided to consider CON as a categorical parameter with two levels: a completely free boundary or a boundary with fully constrained displacements. In the Discussion section we elaborate on the rationale for our choices and address potential consequences. 

Nodes at the equator were constrained to remain in a plane. For cases with the completely free boundary, one node at the sclera equator was also constrained in the anterior–posterior direction to preclude setting up an ill-defined problem. 

Numerical Details

Commercial FE software (ANSYS, ver. 8; ANSYS, Inc., Canonsburg, PA, USA) was used to develop and analyze the models. The process was scripted in Ansys parametric design language. A configuration could be produced, solved, and analyzed without user intervention, typically requiring less than a minute per configuration on a desktop workstation with 32 GB of RAM. 

All tissue regions were meshed with eight-node quadrilateral elements (PLANE 82 in Ansys). Optimal element size was determined in a preliminary mesh refinement study.⁴³ Once sufficient element resolution was determined for a particular geometry, the resolution was quadrupled (element side length divided by two in each direction) to allow for the higher resolution requirements of other configurations. After the study, cases with particularly high strain or stress levels were refined and solved again to verify that the default resolution was sufficient. 

Responses and Dimensionality Reduction

It remains unknown which aspects of ONH biomechanics determine the risk of glaucomatous optic neuropathy. Hence, while it would have simplified things substantially, we deemed it unwise to restrict attention to a small set of responses. We strived for comprehensiveness, and therefore we expanded substantially the ONH biomechanical responses analyzed, studying 98 responses, or outcome measures, compared with the 29 in our previous work.³⁰ A complete list of all responses is given in Supplementary Table S1. 

To deal with such a large number of responses without getting lost, we utilized the dimensionality reduction technique called principal component analysis (PCA).^44,45 PCA rests on the idea of utilizing the covariations between variables, responses in this case, to identify the common variation content in groups of responses. Specifically, PCA involves computing the eigenvectors and eigenvalues of the correlation matrix of the responses. The eigenvectors describe independent patterns in the variation of the responses. In PCA, these new variables are called the principal components (PCs), and are ordered according to the amount of variance they account for. In this sense, PC1 is the variable with the largest variance, PC2 has the second largest variance and is orthogonal to PC1 (i.e., PC1 and PC2 are uncorrelated), PC3 has the third largest variance and is orthogonal to PC1 and PC2 (i.e., PC3 is uncorrelated with both PC1 and PC2), and so on. 

We have demonstrated the application of PCA to the study of ONH biomechanics.⁴⁶ In that study, we found that four PCs captured over 96% of the variance in 25 responses to IOP. In this study, we used PCA on both phases. In the screening phase, we identified PCs accounting for, at least, 96% of the variance. We then identified the factors most strongly influencing these PCs. This process guaranteed that we would identify the most influential factors over all responses. The PCA analysis was repeated for the outcome measures in the second phase and used for finding archetypes and to interpret the results. 

A limitation of PCA is that the PCs do not necessarily have a clear interpretation.⁴⁶ Hence, guided by the PCA and our understanding of ONH biomechanics, for the second phase of this work we also selected 10 representative responses from the set of 98. These responses are not orthogonal and therefore have some degree of redundancy, but they have a clear interpretation and several have historically been of interest in ONH biomechanics research. Thus, results presented in this manuscript are focused on the PCs and 10 representative responses. These were the anterior–posterior lamina cribrosa displacement (LCD), the scleral canal expansion (SCE) at the canal opening, the displacement and rotation of the PPS at its anterior surface 1.7 mm from the axis of symmetry (to mimic the ring 3.4 mm in diameter centered in the canal), the rotation of the scleral canal wall (an in-plane rotation that was measured as the change in the canal wall angle, also referred to as PPS bowing for short),⁴⁷ the absolute and relative displacement of the prelaminar neural tissue (including the retina), the median tensile and compressive strains, and the von Mises stress within the LC. For full details of the definitions of these variables and the rationale for computing them, we refer readers to the papers where they were first introduced.^26,30,46 

For the analyses, the response variables were transformed with a base 10 logarithm, as indicated by a Box-Cox analysis,^48,49 to improve the normality of their distributions and of the residuals, to satisfy the requirements of analysis of variance (ANOVA), and to allow factor effects to be added in an unbiased fashion. The analyses were done in the software package R (v2.12.0).⁵⁰ 

Analysis of Factor Influences

For the screening phase we followed essentially the same approach as we have reported elsewhere.³² Briefly, a design of experiments approach following a two-level fractional factorial 2^24–14 design requiring 1024 configurations was used to sample a subset of the corners of the 24-dimensional factor space. ANOVA was used to determine the statistical significance of the factor and interaction effects. For each response, the percentage of the total sum of squares corrected by the mean was used to represent the approximate contribution of each factor and interaction to the variance of the response. To be deemed influential, a factor or interaction had to contribute at least 5% to the total variance of a response or a PC. In addition, the effects had to be statistically significant (P < 0.01). 

In the second phase, a denser sampling of the factor space was carried out to map in detail the nonlinear relationships between the influential factors and the representative responses. The combinations of factors were chosen using response surface method, again following the approach we have described elsewhere.³³ A total of 7316 combinations were produced into models, simulated, and analyzed. The responses were fitted by polynomial functions of the form  

\(\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}{\rm Response} = f\left( {{x_1},{x_2},...,{x_n}} \right) = {\beta _0} + \sum\limits_{i = 1}^n {{\beta _i}{x_i} + } \sum\limits_{j = 1}^n {\sum\limits_{i = 1}^j {{\beta _{ij}}{x_i}} } {x_j}{\rm{ + \ }}\sum\limits_{k = 1}^n {\sum\limits_{j = 1}^k {\sum\limits_{i = 1}^j {\beta {_{ijk}}} } } {x_i}{x_j}{x_k}{\rm{ + }}\varepsilon \end{equation}

where the x's are the factors, β's are the regression coefficients to be estimated, and ε is the residual.  

Archetypal Analysis

The thousands of models in the response surface analysis represented the wide diversity of ONH sensitivities to IOP, CSFP, and BP possible from the many factors we considered. We were interested in using these models to obtain a better understanding of what this diversity of responses tells us about ONH sensitivity to the three pressures. For this purpose, we used a technique called archetypal analysis.^51–53 Archetypes are cases selected in the multidimensional response space such that all other cases can be represented as convex combinations of the archetypes. In other words, the archetypes are examples of extreme cases. As such, archetypes illustrate aspects of the responses using extreme cases. Whereas a mean response is expected to be near the center of the response space, archetypes are the opposite, being located at extremes. To avoid selecting an outlier as an archetype, archetypal analysis also requires that archetypes themselves be convex combinations of the individual responses and limited to the boundary of the occupied response space. In practice, the archetypes are selected by minimizing the residual sum of squares (RSS) of a representation of all responses as a mixture of archetypes. Computing the number of archetypes is therefore a nonlinear least-squares problem, which is solved using an alternating minimizing algorithm. We used the implementation of nls function in R (v2.12.0).⁵⁰ 

Results

Screening—Gross Analysis of 24 Factors

The screening analysis showed that, among all 24 factors considered, 16 factors and their interactions accounted for between 94.4% and 99.7% of the variance in the responses. These factors were the pressures (IOP, CSFP, and BP), CON, the eye radius, the properties of the sclera (modulus and shell thickness), LC (modulus, depth, thickness, and radius), dura mater (modulus), pia mater (modulus), retina (modulus), optic nerve (modulus), and vessel (modulus). 

Response Surface—Detailed Analysis of the 16 Most Influential Factors

The influence strengths of the most influential factors on the PCs and representative responses are summarized in Figure 2. Overall, IOP was predicted as the most influential factor, with substantial effects on almost every response considered. Particularly strong were its contributions to the variations in the SCE as well as the stress and strain within the LC (more than half the variance). Following IOP, CON ranked as the second most influential factor, with effects mainly on those responses relating to gross ONH geometry. For example, the influence of CON on the LCD was 15 times more than that of IOP. The displacement and rotation of the PPS and the rotation of the canal wall were also more sensitive to CON than IOP. Sclera and lamina moduli were predicted as the third and fourth most influential factors. Sclera modulus affected most of the responses, except for retina displacements. Lamina modulus mostly affected the lamina, especially the stresses, and the retina (as support within the canal). Dura modulus was predicted as the fifth most influential factor, with influence on the displacement and rotation of the PPS three and seven times more than that of IOP. CSFP was predicted as the sixth most influential factor overall, surpassing the lamina radius (seventh) and scleral thickness (eighth). The strongest influence of CSFP, more than twice that of IOP, was on the rotation of the PPS, although it also affected the strains within the LC. BP ranked as the 15th most influential factor, only above the stiffness of the vessel wall, with overall influence on the responses 633 times smaller than that of IOP. 

Figure 2

 

Strength of factor influences as determined by response surface methodology. Columns 1 through 4 present the top four principal components (PCs). Columns 5 through 14 present the 10 representative responses. Columns 15 and 16 present the maximum and average influence of factors on the 10 representative responses (columns 5–14). Rows 1 through 16 present the 16 most influential factors, sorted from highest (top) to lowest (bottom) average influence. Cells are colored according to the strength of a factor influence (row) on a response (column). These were computed as the percentage of a response variance due to each of the factors, with strong influences shown in red and weak influences in blue. Strengths of factor interactions were calculated, but are not shown.

Our models predicted that the overall influence of TLPD on the responses would be 28 times smaller than that of IOP, weak compared with the overall influence of CSFP, which was 16 times smaller than IOP. The strongest influence of TLPD was on the median tensile and compressive strains within the LC, which was still 10 times smaller than that of IOP. 

PCA—Comprehensive View of Many Responses

PCA in the second phase indicated that four PCs accounted for 96% of the variance in 98 responses. Biplots of the top four PCs were used to demonstrate the relationship between PCs, responses, and factors (Fig. 3). PC1 accounted for 45% of the variance and corresponded to a wide range of responses, including the stress and strain within the LC, the SCE, and the displacement and rotation of the PPS. PC1 was most strongly influenced by IOP and the sclera modulus. PC2 explained 37% of the variance and represented the LCD, the rotation of the canal wall, and the displacement and rotation of the PPS. PC2 was most strongly influenced by CON and the sclera modulus. PC3 accounted for 10% of the variance and represented the stress within the LC and the LCD. PC3 was most strongly influenced by the lamina and dura moduli. PC4 accounted for 5% of the variance and represented the stress and strain within the LC. PC4 was most strongly influenced by the LC modulus. 

Figure 3

 

Biplots of the top four principal components (PCs). Left: PC1 and PC2; right: PC3 and PC4. The top four PCs accounted for over 96% of the total variance. A biplot shows two-dimensional projections of the responses (black lines) and factors (red lines). The angle between lines represents the strength of the correlation between variables. Strongly correlated variables are parallel (0°) or antiparallel (180°), and independent variables are orthogonal (90°). All lines have a length of 1 in a 98-dimensional space. Line length in a biplot is the variance accounted for by the two PCs. The factors were not included when computing the PCs and are shown only as covariates to illustrate their relationship with the responses and the PCs. (Readers unfamiliar with principal component analysis or biplots may refer to our previous publication.⁴⁶)

Archetypal Analysis—Illustrative Examples of Factor Influences

To understand the diversity and distribution of responses it would be desirable to do a scatter plot of the model responses. However, it is impossible to show explicitly the 98-dimensional response space. As an alternative, we projected the model responses onto the planes defined by the top four PCs (Fig. 4). Since the PCs account for the vast majority of the variance in responses, this visualization guarantees the maximal spread of the model responses. As expected, the model responses clustered. 

Figure 4

 

Scatter plot of model responses on the top four principal components (PCs). The axes are the same as in Figure 3. Each dot represents the response of one of the 7316 models in the second-phase response surface methodology analysis. The large red numbers are the five archetypes. As expected, the archetypes are spread on the periphery of the response cloud.

The RSS computed using a single archetype was 0.028, decreasing monotonically to 0.007 as the number of archetypes increased to five. The convex space defined by five archetypes, thus, encompassed more than 99.9% of the response space. Further increasing the archetypes from five to eight barely decreased RSS to 0.005, and resulted in some perceived redundancy in the archetypes themselves. Hence, for presentation herein we used the five archetypes listed in Table 2. 

Table 2

 

Factors Selected for Phase 2 and Their Corresponding Levels for the Five Archetypes

 

The five archetypes are highlighted in the PC scatter plots (Fig. 4). As expected, the archetypes are extreme examples and are therefore spread on the periphery of the response cloud. To illustrate the extreme biomechanical behavior of the archetypes, Figure 5 shows how the maximum principal strain distribution in the five archetypes changed as IOP and CSFP increased independently and simultaneously. An interpretation of each archetype is given in the fifth column of Figure 5. 

Figure 5

 

Contour plots of maximum principal strain in the five archetypes. Rows represent the five archetypes. Columns represent the strain response of each archetype to elevations in IOP only, CSFP only, and both IOP and CSFP. The interpretation of each archetype is listed in the rightmost column. Deformations are shown exaggerated five times for clarity. Recall that, by definition of archetype, the responses of all other ONHs are linear combinations of these five archetypes. Note how in all of these cases, the effects of IOP and CSFP did not balance out; they added up. Although the largest retrolaminar strains were also accompanied by a large enlargement of the subarachnoidal space, it was still possible to have substantial retrolaminar deformations without much enlargement.

Discussion

Our goal was to model the sensitivity of the ONH biomechanical environment to acute variations in IOP, CSFP, and BP. Four main predictions arose from this work: First, IOP and moduli of the sclera and lamina are among the most influential factors on the biomechanical environment within the ONH. Second, retrolaminar factors, including CSFP, the dura modulus, and CON, have important influence on ONH biomechanics. Third, IOP and CSFP affect different aspects of ONH biomechanics, and these effects do not balance one another. Fourth, BP has only modest effects on the biomechanics of the ONH. Below we discuss each of these predictions in detail. 

Our sensitivity study revealed that IOP and moduli of the sclera and lamina were among the most influential factors on the biomechanical environment within the ONH. This prediction is consistent with previous results obtained from ONH models without detailed retrolaminar factors.³⁰ The effects of IOP and moduli of the sclera and lamina on ONH biomechanics have been extensively discussed elsewhere,^{26–28,30,32–34,46} and we will not discuss them here. Note that although in this study IOP was predicted to be the most influential factor, this was not the case in some of our previous studies.^30,32 There are three reasons for this. First, some of our previous studies did not consider the interactions of IOP with other factors, as we do here. Not considering such interactions will underestimate the influence of IOP. Second, in this study we monitored a broad set of 98 responses, many more than the 29 of the previous one. For example, we predicted that IOP would have substantial effects on the displacement of the retina (Fig. 2), a response that was not included in the previous study. Third, while some factors influence a few responses, IOP is a consistently influential factor on the majority of the responses. Hence, as more factors and responses are considered, the rank influence of IOP increases. 

Our models also predicted that retrolaminar factors, including CSFP, the dura modulus, and CON, may have important influence on ONH biomechanics (Fig. 2). In fact, these factors were more influential than some previously identified influential factors, such as the scleral thickness and lamina radius. The importance of CSFP has also been identified by two recent computational studies that conducted parametric analysis to investigate the effects of CSFP on ONH biomechanics.^10,11 Both studies predicted that increasing CSFP would induce large deformation within the ONH, especially in the retrolaminar neural tissue. The importance of CSFP conforms to its association with susceptibility for optic neuropathy. Berdahl et al.^14,15 retrospectively reviewed medical records of over 50,000 patients and compared CSFP in subjects with and without glaucoma. They found that CSFP was significantly (P < 0.0001) lower in subjects with normal-tension glaucoma (8.7 ± 1.16 mm Hg) and primary open-angle glaucoma (9.1 ± 0.77 mm Hg) than in the control group (11.8 ± 0.71 mm Hg). Similar observations were found in prospective studies by Ren et al.¹⁶ and Wang et al.¹⁷ Yang et al.⁵⁴ found that chronic reduction of CSFP in monkeys led to decreased retinal nerve fiber layer thickness and neuroretinal rim area of the ONH, features of progressive optic neuropathy. Despite the associations, the mechanistic relationship between CSFP and glaucoma, or other optic neuropathies, is still not fully understood, and further studies are needed. 

Dura modulus was predicted as the fifth most influential factor in ONH biomechanics, even more influential than CSFP, although this varied between responses (Fig. 2). Despite its importance, there is little information about the mechanical properties of the dura mater, especially the portion surrounding the optic nerve. Raykin et al.⁵⁵ recently characterized the mechanical properties of porcine dura mater in vitro. We analyzed their results and calculated a dura modulus of approximately 4 MPa, within the range considered in this study, that is, 1 to 5 MPa. Considering the predicted importance of dura modulus in ONH biomechanics, characterization of the mechanical properties of human dura mater is worthful. 

The movement of the optic nerve would be constrained at the point of orbit exit, but the degree and exact nature of the constraints, and how these are transmitted to the ONH region, remain unclear. Elsewhere simulations have assumed completely free^10,11,30,32 or fully constrained⁴² optic nerves. Acknowledging this uncertainty, and to avoid a potentially biased decision, we considered CON as a categorical parameter with two levels: a completely free boundary and a boundary with fully constrained displacements. These two constraints represent two extremes and the true physiological situation is likely somewhere in between. The fully free condition is also important to study because it mimics the boundary conditions of most ex vivo inflation tests. Surprisingly, our models predicted that CON would rank as the second most influential factor in ONH biomechanics, with effects mainly on those responses relating to gross ONH geometry. Whether CON vary between individuals or may even change with aging or disease is unknown, but seems unlikely. In this sense CON would not be considered as much a risk factor, but as a key parameter that must be determined experimentally and incorporated into biomechanical models. The importance of CON also indicates that it is essential to carry out further studies to better understand the boundary conditions of the optic nerve at the orbit exit, as well as other nerve characteristics that affect how these boundary conditions may interact with the ONH and globe. These include optic nerve tortuosity and tissue incompressibility. Interestingly, recent work shows that changes in gaze may result in optic nerve exerting forces on the ONH and PPS.^42,56,57 Further experiments and nonaxisymmetric models are needed to understand this. 

Our model predictions also showed that IOP and CSFP affected different aspects of ONH biomechanics and that these effects did not balance one another (Fig. 5). This can be explained by the distinct manner in which IOP and CSFP deform the ONH, directly and indirectly. IOP affects the ONH directly by “pushing” the cup and LC posteriorly, and indirectly by deforming the sclera, causing expansion of the scleral canal, which in turn “pulls” the lamina taut from the sides.^58,59 The CSFP also has direct and indirect effects on the LC. We note three mechanisms (Fig. 6). The magnitude of each effect and the overall response when also under elevated IOP will depend on the specific anatomy and mechanical properties of the eye. Considering how different mechanisms of action of IOP and CSFP are, it came as no surprise that their effects generally did not balance out. 

Figure 6

 

Schematic description of three mechanisms by which increases in CSFP cause ONH deformations. Undeformed ONH is shown with continuous lines, and deformed ONH with dashed lines. (a) CSFP acts inwardly compressing the pia mater and the retrolaminar neural tissue within. Due to the Poisson effect, lateral compression may cause expansion in the axial direction, increasing retrolaminar pressure⁸⁶ “pushing” anteriorly on the lamina and causing clockwise rotation of the PPS. The extent of this effect depends on the compressibility of the retrolaminar tissue, which is still not well characterized. (b) CSFP acts outwardly on the dura mater away from the pia mater, causing the known distension of the dural sheath,⁸⁷ rotating the PPS counterclockwise, and displaces the periphery of the lamina posteriorly. (c) CSFP “pushes” the PPS anteriorly, causing flattening of the globe and clockwise rotation of the PPS, and displacing the periphery of the lamina anteriorly. The magnitude of each of these effects will depend on different factors. For example, (a) will depend on the stiffness and thickness of the pia mater, as well as the stiffness and compressibility of the retrolaminar neural tissues; (b) is influenced by the stiffness of the dura and flexibility of the sclera (a combination of its stiffness and thickness). Hence, the various mechanisms will add up or cancel out in various proportions in a given eye.

Our model predictions showed that BP would have only modest effects on the biomechanics of the ONH. This is counter to our original expectation that BP would be influential, considering that the increase in BP (30 mm Hg) was twice that of IOP and CSFP (15 mm Hg). Note that herein we studied the effects of normal variations in BP on overall ONH biomechanics. Variations in BP may also have local biomechanical effects on neighboring neural tissues and on the LC that should be studied and better understood. These studies are important in light of the evidence for a role of ocular perfusion pressure in glaucomatous optic neuropathy.^60,61 

To the best of our knowledge, this is the most comprehensive study modeling the sensitivity of the biomechanical environment within the ONH, in terms of parameters considered and responses examined. Understanding of ONH biomechanics requires parameterized models that incorporate a wide range of anatomic and material factors, pressures, and other boundary conditions. Similarly, it is necessary to study the effects of these factors on a wide range of responses. Not doing so risks missing important aspects of ONH biomechanics. We recognize that considering many responses can also be problematic, as it is possible for results to point to factors that are influential on those responses but that may not have a major role in the neuropathy. However, until the mechanistic link between ONH biomechanics and the neuropathy has been established, we believe that it is better to be comprehensive. We have demonstrated in this work that considering many responses is possible using PCA and archetypal analysis. Note that the use of PCA does not imply that the system biomechanics were linear. PCA can represent nonlinear relationships precisely. This is important because the relationships between responses as identified in this study were nonlinear, consistent with those reported elsewhere.^32,33,47 

Archetypal analysis revealed a wide diversity of sensitivities to IOP and CSFP. Some ONHs were highly sensitive to one pressure and insensitive to the other. From a clinical perspective, a potentially troubling finding is that sometimes the pressure effects were mainly in the retrolaminar region (of IOP in archetype 3, and of CSFP in archetype 5). This means that the effects might be difficult to observe and measure, even with current swept-source optical coherence tomography systems. 

Parametric modeling of the kind we present in this work, and which we have published elsewhere,^29,32 serves to obtain a general understanding of how all eyes work. These models are not intended to represent any specific eyes. There is value in pursuing specimen-specific models that can be inverse fit, or validated against experimental tests, which we have also done.^27,28,62 Specimen-specific models provide excellent information on the particular eyes, but generalizing to a population is problematic and can be highly misleading. We have illustrated potential problems with those generalizations and how parametric modeling can help prevent some of those misunderstandings.⁵⁸ Carefully done, parametric modeling helps provide fundamental new insight into the mechanical behavior of the posterior pole of the eye that would be otherwise unobtainable. A more detailed discussion of the role of parametric modeling in posterior pole biomechanics is outside the scope of this work, and is available elsewhere.⁶³ 

Peak strains predicted by the models in this study were slightly above 5%, which is similar to those predicted by comparable models of human^{10,11,32–34,64} and monkey^58,65 LC, and also to some recent measurements in in vivo⁵⁶ and ex vivo⁶⁶ human and ex vivo porcine⁶⁷ eyes. For comparison with experiments it is important to consider that the strains predicted by our models assume the LC to be homogeneous. We, and others,^10,11,64 have followed this approach, as it is a reasonable approximation of the large-scale behavior of the tissues. As the resolving power of imaging technique increases, experimental studies of ONH biomechanics have reported higher levels of strain at the microscale within the LC, which sometimes exceeded 10%.^56,66–69 Elsewhere we studied the relationship between model detail and predicted LC strain by developing models with a detailed microarchitecture of the beams and pores of the LC.^70,71 We found that models with detailed LCs predicted higher strains, particularly in the pores adjacent to the sclera. However, when observed at a larger, mesoscale resolution, the models predicted LC strains between 2% and 4%, similar to the levels we reported here. 

To compare with our previous studies and extend the lessons and predictions, we adopted the same model simplifications; that is, the geometries were axisymmetric and the mechanical properties were isotropic and linear. A thorough discussion of the limitations of this modeling approach can be found in our previous studies.^{26,29,30,34,46,47,59} While these simplifications may not capture some of the complex behavior of the ONH, they provide a reasonable first approximation. We also aim to inspire others to do more comprehensive analysis of the powerful models they develop. Work is ongoing within our lab^70,71 and others^72–74 to create improved computational models that capture the anisotropy,^75–80 nonlinearity,^{75,76,79,81,82} and inhomogeneity^{79,80,83–85} of the ONH. Another limitation of this work is that the pressure variations were all within the normal range. Given the linear mechanical properties of the model, we find it best to limit the change in pressures to a small range. At elevated or abnormal pressures, the nonlinear material properties would more strongly influence the mechanical behavior of the ONH system. Finally, although we based the factor ranges on the literature and on reasonable assumptions, the choice of factor ranges may affect factor influences and outcome sensitivities. Hence, it is important to interpret the results as an estimate of the factor influences and not take the factor ranking as precise. 

In conclusion, our models predicted that IOP and CSFP are the top and sixth most influential factors on the biomechanical environment within the ONH. IOP and CSFP may affect different aspects of ONH biomechanics, explaining why the overall influence of TLPD was substantially smaller than that of either IOP or CSFP. This suggests that TLPD alone will not be sufficient to predict biomechanically induced glaucomatous neuropathy. CON may drastically affect the responses relating to gross ONH geometry and thus should be accurately determined through experiments. Due to the substantial model simplifications, our results should be considered as an approximation to understand the complex biomechanical environment within the ONH under the simultaneous effects of IOP, CSFP, and BP. 

Acknowledgments

The authors thank Jonathan L. Grimm for assisting with programming, modeling, and plotting. 

Supported in part by National Institutes of Health Grants R01-EY023966, R01-EY025011, P30-EY008098, and T32-EY017271 (Bethesda, MD, USA) and the Eye and Ear Foundation (Pittsburgh, PA, USA) and Research to Prevent Blindness (New York, NY, USA). 

Disclosure: Y. Hua, None; A.P. Voorhees, None; I.A. Sigal, None