purpose. To determine the relative strength, independently and in interaction, of the influence of factors representing the geometry and mechanical properties on the IOP-induced stresses and strains within the optic nerve head (ONH).

methods. A computational model of the eye was developed such that 21 factors could be varied independently or simultaneously. A fractional factorial screening analysis was used to identify the factors and interactions with the largest influences on the lamina cribrosa (LC) and prelaminar neural tissue (PLNT).

results. Nine factors and their interactions accounted for the majority of the variance in the responses (between 95% and 99.8%). These factors were: the properties of the sclera (modulus, eye radius, and shell thickness), LC (modulus and radius), PLNT (modulus and compressibility), and optic nerve (modulus), and IOP. The interactions were stronger on the PLNT than on the LC (up to 16.4% and 9.0% of the response variances, respectively). No factor was the most influential on all the responses or sufficient to ensure high or low levels of strain or stress. Although the modulus of the sclera was among the most influential factors, its effects could be outweighed by other factors.

conclusions. There were strong interactions between and within the geometry and mechanical properties of the tissues of the ONH. This suggests that to ascertain individual susceptibility to IOP it may be necessary to determine several properties of the eye, as well as their interactions. The influential factors and their covariances should be better characterized.

^{ 1 }Although elevated intraocular pressure (IOP) is the primary risk factor for the development of the disease, the mechanisms by which elevated IOP eventually leads to damage and loss of neural function are still unclear and remain controversial.

^{ 2 }

^{ 3 }The controversy is due in part to the wide range of individual sensitivities to elevated IOP. Several recent publications have explored the hypothesis that IOP-induced forces and deformations (stresses and strains) of the tissues of the optic nerve head (ONH) and, in particular, the lamina cribrosa (LC) contribute to the pathogenesis of the disease.

^{ 2 }

^{ 4 }

^{ 5 }

^{ 6 }

^{ 7 }Within this framework, the range of sensitivities to IOP are proposed to be due, at least in part, to differences between individuals in geometry and mechanical properties of the tissues of the ONH. Therefore, to understand individual susceptibility to IOP it is necessary to determine how ONH biomechanics is influenced by the tissue geometry and mechanical properties.

^{ 5 }

^{ 8 }

^{ 9 }and computational

^{ 6 }

^{ 10 }

^{ 11 }models have been developed for this purpose. A previous study examined the influence on ONH biomechanics of 21 factors spanning the geometry, mechanical properties, and loading, and found that ONH biomechanics were influenced most strongly by scleral stiffness and thickness.

^{ 12 }Other factors such as eye size and LC mechanical properties were also influential, but to a lesser degree. Equivalent results were obtained across several aspects of the ONH response, even when the factor ranges were varied. Nevertheless, the analysis has some important limitations because it was performed with a one-factor-at-a-time (OFAT) technique. In OFAT, starting from a baseline configuration, factors are varied one at a time, while keeping all others constant. OFAT designs are common in computational biomechanics due to the simplicity in preparation and analysis.

^{ 13 }OFAT studies, however, privilege the baseline, limit the combinations of factors studied, and provide no information on the possible interactions between factors (i.e., on how the effects of factors depend on each other). This limitation is potentially important because in biological systems factors are often related, vary together, and have effects that depend on each other. Factor interactions have been shown to be important in various areas of biomechanics,

^{ 14 }

^{ 15 }

^{ 16 }

^{ 17 }and the eye is likely not an exception. To the best of the author’s knowledge there are no studies of the role of factor interactions in ONH biomechanics.

^{ 12 }

^{ 18 }The scripts were adapted so that factors could be varied independently and simultaneously. Briefly, a 3-D axisymmetric model with geometry chosen to represent an idealized generic human eye was developed. The ONH was modeled in some detail, whereas the rest of the eye was modeled as a spherical shell of constant thickness. Five tissue regions were defined: corneoscleral shell, LC, prelaminar neural tissue (PLNT, including the retina and choroid), postlaminar neural tissue (ON, including the optic nerve), and pia mater. Although the models were axisymmetric, they represented a 3-D geometry, which should not be confused with a 2-D model.

^{ 12 }

^{ 19 }

^{ 20 }The factors are illustrated in Figure 1and their ranges listed in Table 1 . All tissues were assumed linearly elastic, isotropic, and homogeneous.

^{ 12 }

^{ 18 }All tissues, other than the PLNT, were assumed incompressible. Accordingly, tissue mechanical properties were defined by the Young’s modulus of each of the five tissues and the compressibility (Poisson’s ratio) of the PLNT. The choice of mechanical properties and their consequences are addressed in the Discussion section. In this work, stiff and compliant are used to describe high and low Young’s moduli, respectively. In this sense, stiffness is equivalent to the tissue’s mechanical property and is independent of the geometry. IOP was represented as a homogeneous force on the interior surfaces. IOP was included as a factor to allow evaluation of the magnitude of its effects relative to other factors. The apex of the region representing the cornea was constrained in all directions to prevent displacement or rotation.

^{ 19 }Once sufficient element resolution was determined for a particular geometry, the resolution was quadrupled (element side length divided by 2 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. In every case, it was.

^{ 12 }

^{ 18 }the analysis of the PLNT was limited to the region within 5° of the axis of symmetry to focus on the ONH.

^{21-12}design requiring 512 configurations was selected.

^{ 21 }

^{ 22 }The factor configurations formed an orthogonal array, which means that all factors were distributed in a balanced manner, with an equal number of occurrences of low and high levels for each factor. For example, there were as many configurations with low (256) LC modulus as with high (256) LC modulus (Fig. 2) . The order in which the configurations were preprocessed, simulated, and analyzed was randomized. Ten replicates of a central point configuration (all factors at the midpoint between low and high levels) were added to the design and randomly inserted in the run sequence to check for pure error. Pure error is a measure of the variance intrinsic to the method. In deterministic studies, such as this one, pure error should be zero because of the perfect repeatability of the simulations. Pure error could be nonzero if there had been any drift in the predictions from one simulation to the next. For all the responses, pure error was zero.

^{ 15 }

^{ 21 }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, providing a measure of influence, as is usual in factor analysis.

^{ 21 }

^{ 23 }All factors influenced the responses to some degree, but some of these influences were either very small, or not statistically significant. To be deemed influential, a factor had to contribute at least 5% to the total variance of a response. Interactions had to contribute at least 5% of the variance in a response due to all interactions. To be influential, the contribution of a factor or interaction also had to be statistically significant (

*P*< 0.01) and greater than the residual. The residual was the portion of the corrected total sum of squares that was not accounted for by the factors considered, similar to the

*R*

^{2}measure in a correlation.

^{ 21 }

^{ 24 }The residual was also used to determine confidence intervals and the least-squares mean value for each factor level in interaction plots. Finally, factors could be influential by hierarchy: If an interaction between two factors was deemed influential, then both “parent” factors (i.e., the two interacting factors), were also considered influential,

^{ 22 }

^{ 25 }although this criterion never turned any factor into an influential one. For the contribution threshold, 5% was chosen based on observations of the interaction plots (see the results section). A stricter criterion would have increased the cutoff percentage, classifying as noninfluential some interesting interactions. Conversely, a less strict criterion would have led to more interactions being classified as influential, increasing the complexity of the analysis.

^{ 22 }

^{ 26 }A traditional Box-Cox analysis and plot method was used to determine the optimal transformation for each response. For all responses, it was found that the optimal transformation was a base 10 logarithm. This is equivalent to expressing the responses in decibels and is a common transformation for strain and stress in bioengineering DOE.

^{ 15 }For plotting, the responses were converted back to the original scale. The experiment was designed and analyzed with commercial software (Design-Expert, ver. 7; Stat-Ease, Inc., Minneapolis, MN).

^{ 12 }

^{ 18 }In contrast with previous studies, there were several responses for which other factors were more influential than the sclera, most notably the properties of the PLNT.

^{ 12 }However, there were some rather surprising results: Even when both the sclera and LC were compliant, the magnitude of the strains could be low. Similarly, strains could be high when only the sclera or the LC were stiff. Only when both tissues were simultaneously stiff were the strains guaranteed low. This was due to the influence of interactions whereby the effects of the scleral and LC moduli depend on each other and on other factors. Note that the high magnitudes of strain for some factor combinations are unlikely to be physiologic. This suggests that some of the configurations considered in the experimental design may be outside of physiologic ranges. This is addressed later in the discussion.

^{ 22 }Response ranges were chosen so as to make the interactions clearest. The interaction plots were grouped according to the response: tensile strains in Figures 6 and 7(separated depending on whether the interactions involve the sclera or not), compressive strains in Figure 8 , and stress in Figure 9 .

^{ 2 }

^{ 5 }

^{ 12 }

^{ 20 }have speculated about factor interactions, to the best of the author’s knowledge, this study is the first to focus on factor interactions on the ONH. It is hoped that this information will provide insights into the origin of the range of sensitivities to elevated IOP and that it will help focus future modeling and experimental studies of ocular biomechanics.

^{ 12 }

^{ 18 }

^{ 27 }

^{ 28 }

^{ 29 }

^{ 30 }If the deformation experienced by the sclera for a given level of IOP changes, then what it transmits to the ONH also changes. Under a given load, the deformation experienced by the sclera is a combination of its mechanical properties (modulus) and geometry (mainly its thickness and the eye radius). Since decreases in eye radius or increases in shell thickness may reduce the deformation of the sclera, these factors could be said to change the effective, or structural, scleral stiffness.

^{ 31 }

^{ 32 }

^{ 33 }The origin of this result is still not understood. An often-mentioned hypothesis is that perhaps CCT is a surrogate of other properties of the sclera or LC, such that a structurally weak (either compliant or thin) cornea could correspond to a structurally weak sclera or LC, leading to a mechanically sensitive ONH.

^{ 2 }

^{ 33 }

^{ 34 }

^{ 35 }

^{ 36 }

^{ 37 }However, despite efforts,

^{ 35 }

^{ 36 }

^{ 38 }the relationship between cornea and sclera or LC is still unclear. Oliveira et al.

^{ 38 }found no correlations between CCT and anterior sclera thickness or axial length. Jonas and Holbach

^{ 36 }found that in nonglaucomatous human globes, CCT did not correlate with LC or peripapillary sclera thickness, and concluded that “susceptibility to glaucoma cannot be explained by an anatomic correspondence between corneal thickness and histomorphometry of the optic nerve head.” Wells et al.

^{ 35 }found that corneal hysteresis, but not CCT, or other of the corneal parameters they measured, was associated with increased deformation of the optic nerve surface during elevations of IOP. These studies are valuable, but they may be more fruitful if factor interactions are considered. According to the interactions shown herein, the influence of the sclera thickness and axial length decreases dramatically as the sclera modulus increases, which could mean that sclera thickness and axial length correlate with IOP-induced changes in the ONH surface in eyes with a compliant sclera, but not in eyes with a stiff sclera.

^{ 39 }It is common to show plots in which more than one factor is varied. Investigators, however, generally do not explore the interactions, missing the insight that these may provide, or leaving the readers with the impression that these are intractable complications. This study shows that, despite the complexity of ONH biomechanics, it is possible to study and quantify factor interactions in a systematic way. Computational models, even when simplified, provide an ideal platform to explore ocular biomechanics and identify the key factors and their interactions to inform experimental design and analysis.

^{ 40 }

^{ 41 }

^{ 42 }

^{ 43 }However, analysis of this possibility requires models incorporating CSF.

^{ 12 }

^{ 20 }These results are not contradictory; it is an example of the benefits of two improvements of this work compared with previous work: First, adding the compressive strains as a response of interest revealed influential factors that had been missed because they did not influence the tensile strains. Second, use of a fractional factorial sampling strategy allowed exploration of a much larger region of factor space and consideration of a compliant and compressible PLNT. However, the rationale for allowing PLNT tissue compressibility must be noted. Compressibility in this work was intended to model the changes in PLNT volume associated with IOP increases, for example, due to altered vascular volume and perfusion or axoplasmic flow.

^{ 2 }

^{ 44 }

^{ 45 }

^{ 46 }

^{ 47 }

^{ 48 }Recall that the PLNT in this work included the choroid, which at elevated IOP may be unable to retain volume. Studies of brain tissue have found it to be “nearly incompressible.”

^{ 49 }The differential effects of PLNT properties on tensile and compressive strains could be important because experiments on neurons and astrocytes have established that cell sensitivity to mechanical stimulation depends on the type, magnitude, and temporal profile of the stimulus.

^{ 50 }

^{ 51 }

^{ 52 }

^{ 53 }

^{ 54 }

^{ 55 }

^{ 56 }The magnitude of the effects of the properties of the PLNT was somewhat surprising, especially because so little is known about them.

^{ 11 }

^{ 18 }of the factors and their ranges,

^{ 12 }

^{ 20 }and of the responses analyzed,

^{ 7 }

^{ 12 }

^{ 20 }and so these will not be discussed at length again. Instead, below is a summary of earlier discussions, with a focus on the limitations and considerations most relevant to this work.

^{ 57 }

^{ 58 }

^{ 59 }

^{ 60 }nonlinear,

^{ 27 }

^{ 61 }

^{ 62 }

^{ 63 }and anisotropic.

^{ 27 }

^{ 60 }

^{ 63 }

^{ 64 }

^{ 65 }

^{ 66 }

^{ 27 }

^{ 61 }

^{ 62 }

^{ 63 }This brings about a shift in the influences of the various factors that interact with the scleral modulus. In this work, the use of linear mechanical properties decoupled the effects of IOP and scleral modulus. Consequently, no interactions involving IOP were found. As more complex mechanical properties are incorporated into the models, an effort already in progress, IOP is expected to interact with other factors. The lessons learned from these simplified models may help in understanding models with more complex and physiologically accurate properties.

^{ 53 }

^{ 67 }

^{ 68 }In this work, as in previous studies,

^{ 6 }

^{ 12 }

^{ 18 }the highest magnitudes of strain were predicted in the neural tissue regions, not in the LC where damage to the retinal ganglion cells (RGCs) seems to initiate in glaucoma.

^{ 69 }

^{ 70 }

^{ 71 }

^{ 72 }

^{ 73 }This could be because with a homogeneous LC, the models presented herein could not account for the effects of the LC microarchitecture, which may amplify the levels of strain (Downs J, et al.

*IOVS*2007;48:ARVO E-Abstract 3301). It is also possible that the RGC axons are not equally sensitive to, or equally capable of recovering from, mechanical stimulation in all regions of the ONH.

^{ 19 }

^{ 74 }

^{ 75 }

^{ 76 }The biomechanical effects of the details of the geometry have been studied using OFAT techniques and individual-specific models based on human donor eyes.

^{ 6 }

^{ 7 }

^{ 11 }

^{ 20 }It was found that the details of the geometry had a modest influence compared with that of mechanical properties.

^{ 21 }

^{ 22 }A full factorial analysis (all the combinations) on 21 factors at two levels would have required 2

^{21}= 2,097,152 models. The 2

^{21-12}fractional factorial design enabled the study to be performed with a subset of only 512 configurations (plus centers), a 1/4096th fraction, saving considerable time, but at a cost: It was not possible to resolve all possible factor effects. Those unresolved become confounded, or aliased, with other effects.

^{ 21 }

^{ 22 }A careful design allowed the confounded effects to be only those of higher order interactions. The smaller the fraction of factor combinations studied, the more factor effects are confounded—a property called the design “resolution.” In a resolution V design, like the one used in this work, no main effects or two-factor interactions are confounded with each other.

^{ 21 }

^{ 22 }It was found during the analysis that no three-factor interaction had a strong and statistically significant influence, and therefore the discussion elsewhere was limited to independent factors and their two-factor interactions. Herein, 21 factors were studied, but the methodology can be extended to many more (even hundreds

^{ 77 }). There are other factorial designs requiring fewer runs, such those of Taguchi or Plackett-Burman,

^{ 26 }but it was not the purpose of this study to identify the minimum number of configurations necessary for the analysis. The design was chosen because it has an amenable mix of efficiency, robustness, and clarity.

^{ 22 }A general linear model (equivalent to a multivariate linear regression) was fit to each of the responses by using the parameters estimated by the ANOVA,

^{ 21 }

^{ 22 }

^{ 26 }with excellent correspondence (

*R*

^{2}> 0.97 and

*P*< 0.0001 for all responses). However, a lack-of-fit test at the center configuration was significant (

*P*< 0.01 for all responses). This means that the factorial analysis in this work successfully established overall factor influences, but also that it was unable to determine accurately the response curvatures—not surprising considering the curvature in responses observed in previous OFAT studies.

^{ 12 }

^{ 18 }

^{ 20 }A second-phase analysis with more factor levels, which allows a better fit of the curvature, is already under way.

^{ 12 }

^{ 20 }The choice of ranges and factorial analysis allowed some combinations of factors that produced levels of strain that appear too high to be physiologic. Without more information about the mechanical properties of the tissues of the ONH, of how they vary with other factors, and of which levels of strain or stress are unrealistic, it is impossible to determine a priori which combinations of factors are unrealistic. For this reason, and because the intention was in part to help guide future experimental work, these configurations were not excluded from the analysis.

^{ 49 }

^{ 78 }

^{ 79 }

^{ 80 }

^{ 81 }

^{ 82 }Previous works

^{ 7 }

^{ 83 }have reported and discussed the complexity of the mechanical response of the ONH to acute increases in IOP, the need to differentiate between tensile and compressive strains, as well as the need to distinguish the strain from the stress. This work was focused on the LC and PLNT because previous studies have shown that they exhibit the largest median and peak strains, respectively.

^{ 7 }

^{ 11 }

^{ 12 }Also, the LC is of interest since it is where insult to the RGC axons is believed to initiate.

^{ 69 }

^{ 70 }

^{ 71 }Unfortunately, however, the physiologically relevant response is still unknown.

^{ 2 }

^{ 7 }Therefore, the choice of responses, although based on an understanding of the physiology of the ONH, was ultimately arbitrary. Nevertheless, the author believes that the responses chosen are indicative of some of the fundamental aspects of the mechanical effects of IOP on the ONH, whether or not they represent the actual path or biological effect leading to the disease.

**Figure 1.**

**Figure 1.**

Factor | Code | Units | Low −1 | High +1 |
---|---|---|---|---|

Factors defining the geometry of the eye and ONH | ||||

Eye size (scleral shell internal radius) | A | mm | 9.6 | 14.4 |

Scleral thickness (at canal wall) | B | mm | 0.32 | 0.48 |

Laminar thickness at axis | C | mm | 0.24 | 0.36 |

Retinal thickness (shell) | D | mm | 0.16 | 0.24 |

Scleral thickness (shell) | E | mm | 0.64 | 0.96 |

LC anterior surface radius | F | mm | 0.76 | 1.14 |

Pia mater thickness | G | mm | 0.048 | 0.072 |

LC depth at axis below rim | H | mm | 0 | 0.2 |

Cup-to-disc ratio | J | — | 0.1 | 0.5 |

Canal wall angle to the horizontal | K | deg | 48 | 72 |

Optic nerve angle to the horizontal | L | deg | 64 | 96 |

Peripapillary sclera tapering | M | — | 0 | 1 |

Peripapillary rim height | N | mm | 0.24 | 0.36 |

Cup depth | O | mm | 0.26 | 0.4 |

Loading | ||||

IOP | P | mm Hg | 20 | 30 |

Factors defining the mechanical properties of relevant optic tissues | ||||

PLNT compressibility (Poisson ratio) | Q | — | 0.4 | 0.49 |

Pia mater stiffness (Young’s modulus) | R | MPa | 1 | 9 |

LC stiffness (Young’s modulus) | S | MPa | 0.1 | 0.9 |

Sclera stiffness (Young’s modulus) | T | MPa | 1 | 9 |

PLNT stiffness (Young’s modulus) | U | MPa | 0.01 | 0.09 |

ON stiffness (Young’s modulus) | V | MPa | 0.01 | 0.09 |

**Figure 2.**

**Figure 2.**

**Figure 3.**

**Figure 3.**

**Figure 4.**

**Figure 4.**

**Figure 5.**

**Figure 5.**

**Figure 6.**

**Figure 6.**

**Figure 7.**

**Figure 7.**

**Figure 8.**

**Figure 8.**

**Figure 9.**

**Figure 9.**

Lamina Cribrosa | Prelaminar Neural Tissue | |||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Tensile Strain | Compressive Strain | Von Mises Stress | Tensile Strain | Compressive Strain | Von Mises Stress | |||||||||||||||||

50th | 95th | 50th | 95th | 50th | 95th | 50th | 95th | 50th | 95th | 50th | 95th | |||||||||||

Factors | 94.1 | 87.1 | 95.3 | 93.5 | 96.7 | 95.5 | 91.5 | 8.7 | 87.0 | 82.4 | 86.1 | 87.0 | ||||||||||

Interactions | 3.7 | 9.0 | 2.2 | 3.8 | 1.8 | 2.8 | 6.6 | 16.4 | 12.8 | 15.7 | 13.7 | 11.2 | ||||||||||

Factors + interactions | 97.8 | 96.1 | 97.5 | 97.3 | 98.5 | 98.4 | 98.2 | 95.0 | 99.8 | 98.1 | 99.8 | 98.2 |

*H*]. 2003;217:341–348. [PubMed]