Finite element modeling of the viscoelastic responses of the eye during microvolumetric changes

A linear viscoelastic finite element model was built to investigate factors that influenced the intraocular pressure (IOP) elevations due to micro-volumetric changes in the eye at three different rates. The viscoelastic properties of the cornea and the sclera, including the instantaneous modulus, equilibrium modulus, and relaxation time constants, parametrically varied to examine their effects on IOP elevations at different rates of volumetric changes. The simulated responses were in good agreement with the previously reported experimental results obtained from porcine globes, showing the general trend of higher IOP elevations at faster rates. The simulations showed that all viscoelastic properties influenced the profile of the dynamic IOP due to volumetric changes, and the relative significance of a specific parameter was highly dependent on the rate of change.


INTRODUCTION
Glaucoma is one of the leading causes of blindness worldwide [1]. It is known to be associated with elevated intraocular pressure (IOP), although the mechanisms of IOPinduced optic nerve damages are not well understood. The longstanding focus in the glaucoma field has been on the steady-state IOP [2], but it is also known that the IOP fluctuates during the day [3], as well as during normal physiological activities such as postural change [4,5], valsalva maneuver [6,7], water drinking [8], and eye movement [9]. Previous studies have shown that these short-term IOP fluctuations may be linked to glaucoma risk [8,10,11].
Very little has been done to understand the response of the eye to the short-term changes of volume and pressure. Recently, our group examined IOP elevations in porcine globes due to short-term micro-volumetric infusions. We found a strong rate-dependent response showing that a fast volumetric change induced a significantly higher IOP change in the same eye compared to slow volumetric changes [28]. The current study aims to build a finite element model of the whole corneoscleral shell (omitting the optic nerve head) to investigate how cornea and sclera viscoelastic properties influence the eye's response to microvolumetric changes and how each factor contributes to the rate-dependent response.
The whole eye models have been developed in the past for understanding ocular physiology. Srodka [29] utilized an axisymmetric finite element of the whole eye model to evaluate the accuracy of the standard clinical method for IOP measurement, the Goldmann applanation tonometry. Anderson et al. [30] utilized a 3D hyperelastic whole eye model to investigate the influence of the boundary conditions at the cornea and sclera conjunction on the apical rise of the cornea during inflation. Our current study incorporated the general formulation of these previous whole eye models such as using an axisymmetric geometry and omitting the ONH.

Finite Element Modeling
An axisymmetric model of the corneoscleral shell of a porcine eye was constructed in COMSOL (v4.3, Burlington, MA, USA) using a series of ellipses assuming symmetry along the optic axis. The ellipsoidal building blocks allowed for realistic representation of the thickness variations in the corneoscleral shell. The central cornea thickness (CCT) was set to 0.96 mm [31], with a radius of curvature of 7.5 mm, an assumed eccentricity of 0.5, and a white-to-white (WTW) diameter of 13.5 mm. The scleral radius of curvature was set to 12.8 mm, the anterior thickness to 0.89 mm, the equatorial thickness to 0.58 mm, and the posterior thickness to 1.12 mm. Because of the oblate spheroidal shape of the porcine sclera, the ratio between the posterior radius and the equatorial radius was set to 0.86. These dimensions were the average measurements in porcine eyes (SiouxPreme Packing Co, Sioux City, IA) acquired for other studies in our laboratory.
The cornea and sclera were divided by a radial line going through the center of the sclera ellipse. A small transition zone at the cornea-sclera junction was implemented with a finer mesh to avoid convergence difficultties at the boundary of material property change. A filling material was used to fill the inner space of the ocular shell and to couple the intraocular volumetric change to the mechanical responses of the shell (as described subsequently). Figure 1(a) shows the constructed geometry of the whole eye as described above, and Figure 1(b) illustrates the typical finite element mesh.
The mesh was constructed using quadratic triangular elements. A displacement mesh convergence test was performed with five different mesh densities that monitored five different locations to ensure mesh accuracy. The mesh density chosen for this study had approximately 11,900 elements for a total of 54,000 degrees of freedom (Figure 1(b)). A time-dependent solver in COMSOL was utilized for solving the constitutive equations.

Simulated Infusions
Our laboratory has previously reported an experimental infusion study on porcine whole globes at three different rates (the details of the experiments are provided in the referenced article) [28]. In the present study, the simulations were performed following the design of the infusion experiments. Specifically, three different infusion rates (as shown in Table 1) were implemented to represent the typical time scales seen in short-term IOP fluctuations (seconds to minutes) [28]. The total volume change was set to 15 μL.

Using Thermal Expansion to Simulate Intraocular Volumetric Change
In order to simulate intraocular volume increments, it was computationally convenient to use a thermal expansion model although in reality there is no temperature dependence being studied. However, by assuming that the tissue and fluids inside the ocular shell are incomepressible, the volume change can be modeled as a simple uniform volumetric thermal expansion of an intraocular filling material, which can be implemented in the COMSOL software to conveniently provide mechanical coupling between the shell and the intraocular material. The thermal expansion model assumed that all mesh points within the enclosed filling material had a heat source and the total power (P) was computed by the following relation: (1) where dV/dt is the volumetric change rate, C p is the specific heat, ρ is the density and α V is the coefficient of thermal expansion for the filling material. The parameters for low density polyethylene (LDPE) was used to simulate the filling material because of its high coefficient of thermal expansion (α V = 4 × 10 −4 1/°C) [32]. The volumetric change rate (dV/dt) resulted from the thermal expansion was set equal to the experimental infusion rate assuming no leakage. A fixed point was placed at the center of the filling material to ensure that it expands equally in all directions. The filling material was thermally insulated at the outer boundary to ensure conservation of energy.

Constitutive Model
The cornea and sclera were modeled using a nearly incompressible linear viscoelastic generalized Maxwell model. The stress-strain relationship is defined as follows [32]: (2) where E equ is the equilibrium modulus, E i is the relaxation modulus of the i-th branch, τ i is the time constant of the i-th branch, ε is the strain, and ε̇ is the strain rate. A two-branch model (i.e., m = 2) was adopted [16]. The time constants (τ 1 and τ 2 ) are denoted as the short term time constant (τ s ) and the long term time constant (τ l ) in the later discussions for clear differentiation. The instanttaneous (E inst ) is defined as [14]: (3) The branch relaxation moduli were set to be equal (i.e., E 1 = E 2 ). Therefore, the branch relaxation moduli can be found as: (4) The instantaneous and equilibration moduli of the sclera were set to be 5 times of those of the cornea to ensure a stiffer response of the sclera commonly reported in experimental results and simulation studies [28,34,35]. All other parameters were allowed to vary independently.

Parametric Study
The values of the material parameters of ocular shell were varied within a range of the plausible values for the porcine eye (as shown in Table 2) to evaluate the effects that each parameter has on the short-term IOP elevations. To the best of our knowledge, experimental measurements of the viscoelastic properties of the porcine cornea and sclera including the instantaneous modulus, equilibrium modulus, and time constants have not been reported. The selected baseline values were estimated from approximating the model predictions with the average results obtained from the infusion experiments in porcine eyes [28].

RESULTS
A comparison of the simulated IOP/infusion volume response (solid lines) and the average experimental responses in the porcine eyes (individual markers) is presented in Figure 2. All viscoelastic properties were set at the baseline values. Selected volume levels (11 levels evenly distributed between 0 to 15 μL) were used for plotting the experimental response for visual clarity. Most of the markers fell on the corresponding curves, suggesting a good agreement between the model simulation and the experimental results could be simultaneously achieved at all three rates. Figure 3 presents the effects of instantaneous modulus and equilibrium modulus on IOP elevations during micro-volumetric changes in the eye at different rates. The instantaneous moduli of the cornea and the sclera were varied proportionally (i.e., 1:5) while all other parameters were kept at their baseline. In general, IOP elevations were higher at higher modulus regardless of the infusion rates. Figure 3(a) shows that the IOP elevations increased almost linearly with the instantaneous moduli of the cornea and sclera at all rates. The slope however was different for different rates with a larger slope at the higher infusion rate. Figure 3(b) shows the effects of the equilibrium modulus. The response was less linear for the fast and intermediate infusion rates, particularly at the lower modulus range.
Treating each data point in Figure 3(a) as an individual eye, the data points form a group of simulated eyes that only differ in instantaneous modulus. The relationships of the IOP elevations at different infusion rates in the group of eyes are plotted in Figure 4(a), showing that the IOP elevations at different infusion rates are predicted to be linearly correlated in a group of eyes that only differ in the instantaneous modulus. A similar plot for the eyes that only differ in equilibrium modulus (using data points in Figure 3(b)) is shown in Figure  4(b). The relationships of IOP elevations at different infusion rates become nonlinear for this case. Figure 4(c) shows the experimental results obtained from 11 porcine eyes, showing that the IOP elevations at different infusion rates were approximately linearly correlated in the tested eyes (R = 0.99, p < 0.001, for fast and intermediate correlation; R = 0.89, p < 0.001, for fast and slow correlation) [28].  (Figure 5(b)). A larger long term time constant results in larger IOP elevations for the slow and intermediate rates. Figure 6 shows the comparative influence of each of the corneal and viscoelastic properties on the IOP elevations, separately presented for each infusion rate. The x-axis shows the value for all parameters linearly scaled from its minimum value ( 1) to maximum value (+1). It is found that the instantaneous modulus has the largest influence for the fast and intermediate rates, but its influence diminishes for the slow infusion rate, for which the equilibrium modulus takes the dominant role. The time constants generally have smaller effects than the moduli, but their influence could be comparable to one of the two types of moduli depending on the rates. Figure 7 shows the distribution of the von Mises stresses and principal strains in the cornea and the sclera at the end of the fast infusion. The stresses and strains were apparently different in the cornea as compared to those in the sclera, due to the discrepancy in the material properties. The stresses and strains also showed variations through the thickness as contrasted to what a thin-shell model would predict [22].

DISCUSSION
In this work, we have built a finite element model of the whole ocular shell to investigate factors influencing the short-term volume/pressure changes in the eye. We found that a simple linear viscoelastic model with two branches simulating the short-term and long-term relaxation processes could adequately reproduce the experimental data obtained from the infusion experiments performed in porcine eyes (Figure 2). Although the three different infusion rates used in this study are all within the time scale of short-term IOP fluctuations seen in the eye (i.e., a time scale shorter than a few minutes), there was a clear ratedependent response suggesting that the eye responds differently to micro-volumetric changes that occur within seconds, 10's of seconds, or minutes. The larger the rate of change within this time scale, the higher the IOP elevations observed in both experiments and simulations. This rate-dependence confirmed the viscoelastic responses of the ocular shell, as well as the relevant range for the time constants of the relaxation processes in the ocular shell (i.e., seconds to minutes), which were found using stress relaxation tests on tissue strips [16]. Both experiments and simulations reported a linear ΔIOP/volume relationship, which was also experimentally observed by Pierscionek et al. in fresh porcine eye [36].
Based on the finite element model, a higher instantaneous modulus, i.e., a stiffer response to a step load, would result in a larger final ΔIOP for a given infusion volume regardless of the rate of change (Figure 3(a)). A similar positive relationship between ΔIOP and equilibrium modulus was also shown (Figure 3(b)). The relationship between instantaneous modulus and ΔIOP relationship was linear (i.e., a constant slope) while there was a degree of nonlinearity in the relationship between equilibrium modulus and ΔIOP, although that was mostly seen in the fast rate response.
In our previous experimental study, we found a strong correlation between the IOP elevations measured at different rates across different eyes (Figure 4(c)), that is, an eye that had a larger ΔIOP at the fast rate would also have a larger ΔIOP at the slower rates. Interestingly, this correlation was strongly linear (R = 0.99 or 0.89) according to the experimental data. This linear correlation could be better predicted by the variance in instantaneous modulus (Figure 4(a)) rather than the variance in the equilibrium modulus which corresponded to a more nonlinear response (Figure 4(b)).
The short term time constant had a positive relationship with ΔIOP during fast infusions, but a flat response for intermediate and slow infusions ( Figure 5(a)). This time constant represents the time scale of the fast relaxation processes in the tissue. Because of the selected range of this constant in the present study (0.05 to 0.65 sec), it is expected that it does not influence the intermediate and slow rates whose time scales were far beyond this range. The long term time constant represents the slower relaxation process. As shown in Figure 5(b), the ΔIOP in the intermediate and slow infusions were noticeably affected by this time constant, but the ΔIOP in the fast infusion had a flat response because its short duration.
As discussed above, all viscoelastic properties affected the volume/pressure response, and the relative significance of each specific parameter was rate-dependent. For example, at the faster rates (e.g., the fast and intermediate infusion rates used in the present study), the instanttaneous modulus had the largest influence on the magnitude of ΔIOP (Figures 6(a) and (b)). As the rate decreased, the equilibrium modulus started to dominate (Figure 6(c)). In addition, the short-term time constant had minimal effect on intermediate and slow infusions, but significant effect on fast infusion; while the long-term time constant showed the opposite trend.
Previous computational work has shown that the posterior sclera is a major biomechanical structure that strongly influences the mechanical states of the optic nerve head where the glaucomatous injuries occur [12]. The results from the present study showed that the viscoelastic properties of the corneoscleral shell could also influence the dynamic profile of IOP. For example, a larger instantaneous modulus or equilibrium modulus (i.e., a stiffer corneoscleral shell) would result in higher IOP fluctuations due to micro-volumetric changes in the eye during daily activities. Although the pathophysiolgical consequences of IOP fluctuations are not fully understood, it is generally believed that they are detrimental to the delicate tissue structures in the eye, especially when they are at high magnitudes and fast acting [37,38]. It is therefore important to understand and characterize the viscoelastic properties of the ocular shell and their relationship to the dynamic IOP. This work has several limitations. First, the eye is actually a "leaky" shell with a pressuredependent outflow [39]. For simplicity, this pressure-dependent outflow was not implemented in the current model. Instead, a net flow rate was used as the infusion rate. This could have introduced certain inaccuracy in the calculations of the ΔIOP; however, this error was likely small due to the short time duration of the simulated experiments. Second, the relaxation moduli for the branches in the Maxwell model were set to the same values. Although this assumption has been used in the past [14], future experimental studies as well as modeling work are needed to verify this assumption and investigate its implications. Third, a linear viscoelastic model was used and the reported nonlinearity in the stress-strain relationship of the cornea and sclera was neglected. Because the corneoscleral shell is primarily a collagenous tissue, mechanical nonlinearity is expected and has been confirmed in the past. Constitutive modeling that accounts for the nonlinear responses would thus be more desirable [40][41][42][43][44]. However, in the present study, the IOP was always above the baseline of 15 mmHg during the simulated infusion studies, which likely positioned the eye close to or beyond the toe region of the corneoscleral response [28]. Therefore, the linear approximation, which simplifies the model, is likely sufficient in the present study. Fourth, previous studies have shown the importance of collagen microstructure (for example, anisotropic collagen fiber alignment) in affecting the stress and strain distributions in the cornea and sclera [20,[45][46][47][48]. Future work should investtigate how heterogeneous microstructural and material properties could affect the eye's response to dynamic microvolumetric changes.

CONCLUSION
In conclusion, this study utilized finite element models of a porcine corneoscleral shell to investigate the effects of the shell's viscoelastic properties on the IOP elevations due to volumetric changes in the eye. The major finding was that the instantaneous and equilibrium moduli, as well as the short-term and long-term time constants, had an effect on IOP elevations, and the relative significance of the specific parameters was highly dependent on the rate of change.     The effects of (a) the short term time constant and (b) the long term time constant on IOP elevations at different infusion rates.  Color maps of von Mises stress (a) and first principal strain (b) in the cornea and the sclera at the end of the fast infusion (15 μL/s for 1 second).