**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.

^{1}is characterized by a particular pattern of irreversible damage to the retinal ganglion cell axons.

^{2}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}

^{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.

^{30}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).

^{30}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.

^{31}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.

^{30}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.

^{35}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.

^{36}

^{34}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.

^{37}We also did not account for the retinal vessels' exit from the nerve.

^{30}The geometric factors are illustrated in Figure 1 and their ranges listed in Table 1.

^{30}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.

^{30}There is little information about the stiffness of dura mater. Based on previous measurements of the Young's modulus of human spinal cord,

^{38}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,

^{38}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.

^{39–41}After the analysis on IOP, CSFP, and BP, we also tested the overall influence of the translaminar pressure difference (TLPD = IOP − CSFP).

^{10,11,30,32}and with an optic nerve fully constrained at the end.

^{42}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.

^{43}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.

^{30}A complete list of all responses is given in Supplementary Table S1.

^{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.

^{46}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.

^{46}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),

^{47}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}

^{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).

^{50}

^{32}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).

^{33}A total of 7316 combinations were produced into models, simulated, and analyzed. The responses were fitted by polynomial functions of the form

*x*'s are the factors,

*β*'s are the regression coefficients to be estimated, and

*ε*is the residual.

^{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).

^{50}

^{30}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.

^{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.

^{16}and Wang et al.

^{17}Yang et al.

^{54}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.

^{55}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.

^{10,11,30,32}or fully constrained

^{42}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.

^{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.

^{60,61}

^{32,33,47}

^{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.

^{58}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.

^{63}

^{10,11,32–34,64}and monkey

^{58,65}LC, and also to some recent measurements in in vivo

^{56}and ex vivo

^{66}human and ex vivo porcine

^{67}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.

^{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.

**Y. Hua**, None;

**A.P. Voorhees**, None;

**I.A. Sigal**, None

*. 2006; 90: 262–267.*

*Br J Ophthalmol**. 2007; 48: 3632–3644.*

*Invest Ophthalmol Vis Sci**. 2005; 46: 2663–2670.*

*Invest Ophthalmol Vis Sci**. 1999; 18: 39–57.*

*Prog Retin Eye Res**. 1991; 109: 1090–1095.*

*Arch Ophthalmol**. 2002; 120: 1268–1279.*

*Arch Ophthalmol**. 2000; 19: 297–321.*

*Prog Retin Eye Res**. 2012; 96: 597–603.*

*Br J Ophthalmol**. 2010; 90: 203–209.*

*Exp Eye Res**. 2016; 57: 1901–1911.*

*Invest Ophthalmol Vis Sci**. 2016; 139: 031003.*

*J Biomech Eng**. 2014; 55: 4105–4118.*

*Invest Ophthalmol Vis Sci**. 2017; 58: 2739–2745.*

*Invest Ophthalmol Vis Sci**. 2008; 115: 763–768.*

*Ophthalmology**. 2008; 49: 5412–5418.*

*Invest Ophthalmol Vis Sci**. 2010; 117: 259–266.*

*Ophthalmology**. 2012; 119: 2065–2073.*

*Ophthalmology**. 2003; 87: 949–951.*

*Br J Ophthalmol**. 2013; 23: 139.*

*Eur J Ophthalmol**. 2007; 48: 5068–5084.*

*Invest Ophthalmol Vis Sci**. 2007; 48: 4597–4607.*

*Invest Ophthalmol Vis Sci**. 2011; 52: 9431–9437.*

*Invest Ophthalmol Vis Sci**. 1999; 32: 579–584.*

*J Biomech**. 2001; 23: 215–225.*

*Curr Eye Res**. 2006; 39: S385.*

*J Biomech**. 2007; 85: 312–322.*

*Exp Eye Res**. 2009; 8: 85–98.*

*Biomech Model Mechanobiol**. 2009; 8: 99–109.*

*Biomechan Model Mechanobiol**. 2016; 149: 40–47.*

*Exp Eye Res**. 2005; 46: 4189–4199.*

*Invest Ophthalmol Vis Sci**. 2010; 90: 70–80.*

*Exp Eye Res**. 2009; 50: 2785–2795.*

*Invest Ophthalmol Vis Sci**. 2011; 52: 5497–5506.*

*Invest Ophthalmol Vis Sci**. 2004; 45: 4378–4387.*

*Invest Ophthalmol Vis Sci**. 2014; 11: 16.*

*Fluids Barriers CNS**. 2005; 89: 22–34.*

*Morphologie**. Cambridge, UK: Cambridge University Press; 2007.*

*Introductory Biomechanics: From Cells to Organisms**Paper presented at the 2003 Summer Bioengineering Conference, June 25–29, Key Biscayne, Florida, United States*.

*. 1981; 20: 564–566.*

*Invest Ophthalmol Vis Sci**. Philadelphia: JB Lippincott; 1993: 37–50.*

*The Optic Nerve in Glaucoma**. 2013; 27: 97–106.*

*Saudi J Ophthalmol**. 2016; 57: 2452–2462.*

*Invest Ophthalmol Vis Sci**. Toronto, Canada: University of Toronto; 2006.*

*Human Optic Nerve Head Biomechanics: An Analysis of Generic and Individual-Specific Models Using the Finite Element Method*[doctoral dissertation]*. 2011; 52: 5457–5464.*

*Invest Ophthalmol Vis Sci**. New York, NY: Springer Science + Business Media, LLC; 2007.*

*Analysing Ecological Data**. 2012; 53: 4270–4278.*

*Invest Ophthalmol Vis Sci**. 2014; 33: 1381–1389.*

*IEEE Trans Med Imaging**. Boca Raton, FL: CRC Press; 2004.*

*RSM Simplified: Optimizing Processes Using Response Surface Methods for Design of Experiments**. Hoboken, NJ: John Wiley & Sons; 2008.*

*Design and Analysis of Experiments**R: A Language and Environment for Statistical Computing*. Vienna, Austria: R Foundation for Statistical Computing; 2013.

*. 1994; 36: 338–347.*

*Technometrics**. 2015; 12: 20141118.*

*J R Soc Interface**. 2017; 42: 568–574.*

*Curr Eye Res**. 2014; 55: 3067–3073.*

*Invest Ophthalmol Vis Sci**. 2017; 16: 33–43.*

*Biomech Model Mechanobiol**. 2016; 57: 5825–5833.*

*Invest Ophthalmol Vis Sci**. 2016; 57: 4979–4987.*

*Invest Ophthalmol Vis Sci**. 2011; 52: 9023–9032.*

*Invest Ophthalmol Vis Sci**. 2009; 88: 799–807.*

*Exp Eye Res**. 2010; 51: 3399–3404.*

*Invest Ophthalmol Vis Sci**. 1995; 113: 216–221.*

*Arch Ophthalmol**. 2010; 43: 254–262.*

*J Biomech**. Amsterdam, The Netherlands: Kugler Publications; In Press.*

*Biomechanics of the Eye**. 2016; 19: 591–602.*

*Comput Methods Biomech Biomed Engin**. 2010; 51: 295–307.*

*Invest Ophthalmol Vis Sci**. 2014; 55: 1–15.*

*Invest Ophthalmol Vis Sci**. 2016; 57: 2666–2677.*

*Invest Ophthalmol Vis Sci**. 2017; 53: 123–139.*

*Acta Biomater**, Proceedings of SPIE Vol. 10067 (SPIE, February 21, 2017), 100670B.*

*Optical Elastography and Tissue Biomechanics IV**. 2017; 58: 278–290.*

*Acta Biomater**. 2017; 58: 5336–5346.*

*Invest Ophthalmol Vis Sci**. 2015; 56: 2031–2042.*

*Invest Ophthalmol Vis Sci**. 2013; 12: 941–963.*

*Biomech Model Mechanobiol**. 2012; 44: 99–109.*

*Mech Mater**. 2009; 131: 051012.*

*J Biomech Eng**. 2008; 130: 041017.*

*J Biomech Eng**. 1991; 32: 2244–2258.*

*Invest Ophthalmol Vis Sci**. 2004; 78: 609–623.*

*Exp Eye Res**. 2017; 58: 735–744.*

*Invest Ophthalmol Vis Sci**. 2015; 6: 4705–4718.*

*Biomed Opt Express**. 2005; 46: 1286–1290.*

*Invest Ophthalmol Vis Sci**. 1972; 14: 29–39.*

*Exp Eye Res**. 1969; 67: 417.*

*Trans Am Ophthalmol Soc**. 1981; 99: 137–143.*

*Arch Ophthalmol**. 1991; 10: 877–888.*

*Curr Eye Res**. 1995; 36: 1163–1172.*

*Invest Ophthalmol Vis Sci**. 2003; 87: 777–781.*

*Br J Ophthalmol*