Dynamics Modeling of Corneal Deformation under Air Puﬀ with Linear and Nonlinear Viscoelastic Model

Background: The aim of the study is to model the corneal dynamic deformation under an air puﬀ excitation. The deformation response of the cornea was modeled by using linear and nonlinear viscoelastic models. The corneal deformation responses generated from the linear and nonlinear viscoelastic model were correlated with the clinical results, which were obtained from Corneal Visualization Scheimpﬂug Tonometer (Corvis ST) to evaluate the comparable biomechanical parameters of the cornea. Methods: A prompt deformation occurs when the external force applied to the cornea. Then a continuous deformation follows. A simple mass, spring and dashpot system were used to model human eyeball. Results: In linear viscoelastic model, the corneal elastic stiﬀness commanded behavior of the corneal deformation and its maximum, when the viscous component aﬀected for its lateral shifting and marginally alter the magnitude.Whereas, in the nonlinear viscoelastic model, the corneal material nonlinearity commanded the behavior and maximum of the corneal deformation, while the viscous component marginally contributed for its lateral shifting and demonstrated the minimum aﬀect on the magnitude and form. A multi-objective genetic algorithm-based optimization procedure was used to identify the material properties in the nonlinear viscoelastic model for 29 eyes of 20 normal people. Conclusion: The corneal deformation response model with nonlinear viscoelastic model showed to have a better ﬁt with the corneal dynamic deformation behavior under an air pulse excitation. The biomechanical properties of the cornea in vivo can be evaluated by using and analysing dynamic deformation of the cornea under an air puﬀ excitation model.


Background
Cornea accounts for 45 diopters of the total 60 diopters of ocular power [9,18] and is the most important sensory organs that convert light into electrochemical signals so that our brain cans interpret. The cornea is a viscoelastic material instead of an elastic material [28,20,44], it shows both elastic and viscous characteristics when undergoing deformation. The nonlinear elastic behavior has been reported in the literature [14,15,44] and the nonlinear behavior is attributed by the collagen arrangement in the corneal stroma layer [6,19]. The shape, thickness and biomechanical properties of the cornea can be altered by disease, surgery and injury, which can lead to serious changes in the visual performance of the eye [7,16,10,12,3]. A better understanding of the corneal biomechanical properties can help to improve the diagnosis of keratoconus [1,3], prediction of refractive surgery outcome [11], glaucoma development and progression prediction [24,25] and improve the accuracy of intraocular pressure measurement (tonometry) [8,26,36,40]. Corneal Visualization Scheimpflug Tonometer (Corvis ST; Oculus, Wetzlar, Germany) [2] is a clinical device, which is combination of a noncontact air puff tonometer and an ultrahigh speed Scheimpflug camera since 2011. Corvis ST emits a standardized and collimated air puff (3.05-mm diameter) and captures dynamic deformation of the cornea under the air puff at 4330 frames/sec. Certain parameters can be measured after the analysis of the captured images, which are related to corneal biomechanics. Those are the deformation amplitude, radius of curvature and time for the corneal deformation at highest concavity, first and second applanation length, velocity and time [21]. Amid those parameters, amplitude of the corneal deformation holds promise to provide relevant parameters related to corneal biomechanical properties [35]. The device is still under development and new parameters are being introduced and reported which is current not available in Corvis ST [42,23,31]. The corneal dynamic deformation response from the Corvis ST provides a very useful information for the study of the corneal biomechanics under air puff excitation. In this study, the corneal dynamic deformation response under air puff excitation is modeled using linear and nonlinear viscoelastic model. The proposed nonlinear viscoelastic model is then used to estimate the biomechanical properties of the cornea by using optimization with multi-objective genetic algorithm.

Linear viscoelastic model
When the corneal is loaded by the external air puff, the cornea demonstrates instantaneous deformation (pure elastic deformation) followed by progressive deformation (time-dependent viscoelastic behavior) [30,13]. These behaviors can be modeled with a simple mass, spring and dashpot system as shown in Fig. 1. The spring represents the pure elastic resistance and the dashpot represents the time-dependent viscous resistance to an external load, respectively.The cornea is supposed to be an isotropic membrane in the form of sphere with uniform thickness, whereas the deformation of the cornea is supposed to be axisymmetric. The equation of motion below can be used to model the dynamic apex deformation y(t) of the cornea under the air pulse field F(t).
Here c is the corneal damping coefficient, m is the mass of the cornea, k stands for corneal elastic coefficient (corneal stiffness).Multiplication of the air puff area and the experimental air puff pressure result by Ariza-Gracia et al. [4] is used to model the external force (Corvis ST air puff nozzle diameter = 3.05 mm and the area of air puff = π(3.05 2 /4) = 7.306 mm 2 ) as shown in Fig. 2.
The experimental air puff data is fitted with a second order Gaussian model, where a 1 , a 2 , b 1 , b 2 , c 1 and c 2 are the fitting coefficients of the experimental air puff force profile in equation (2) and their corresponding values are a 1 = 0.1737, a 2 = 0.08099, b 1 = 0.01591, b 2 = 0.01028, c 1 = 0.004393 and c 2 = 0.003599, respectively. In the model, the corneal mass is the equivalent corneal mass that involved in the corneal deformation and is estimated to be 15 mg [38]. The nominal value of viscous damping coefficient is assumed to be 0.05 Ns/m [22]. For comparison, the viscous damping coefficient c was varied from 0 Ns/m to 0.25 Ns/m to examine the influence to the corneal dynamics deformation. The value of the corneal stiffness is estimated to be in the range of 50 N/m and 200 N/m depending on the intraocular pressure (IOP) [32,33]. The stiffness of the cornea is defined by the applied force versus the corneal displacement along the direction of the applied load in Goldmann Applanation Tonometry (GAT) [32]. The GAT equivalent corneal stiffness is a constant while increasing the corneal displacement under small deformation and is linear with respect to IOP [32,33,27]. Goldmann Applanation Tonometer was introduced by Goldmann and Schmidt in 1957 [17] and is considered to be a gold standard in tonometry. It is the most widely accepted method of measuring intraocular pressure, the IOP is determined by measuring the force that used to flatten a predetermined applanation contact area of 7.35 mm2 (diameter of 3.06 mm) based on the Imbert-Fick law as shown in Fig. 3 [17]. The law states that the internal fluid pressure that acts on an infinitesimally thin, elastic and homogenous membrane sphere is equal to the pressure that is used to flatten a small portion of the membrane.
where F GAT is the force to applanate the anterior corneal surface, A is the applanated area at the posterior corneal surface and r is the radius of the applanated area, which is equal to 1.53 mm in GAT principle. The corneal apex displacement (y) along the direction of the applied load in GAT can be calculated by using Pythagorean Theorem as shown in equation (4).
where R is the corneal radius of curvature and is taken to be equal to the average human corneal radius of curvature of 7.80 mm by in this study [5]. The corneal stiffness can be obtained by equation (5)

Nonlinear viscoelastic model
The linear spring in the equation (1) are replaced with a nonlinear spring and the spring material nonlinearity is modeled with a nonlinear stress-strain relationship proposed by Woo et al. [44], where ε is the air puff indentation strain and can be calculated using ε= y/R [44], α and β are the corneal material nonlinearity coefficients. α represents the scaling factor and β represents the nonlinear component between stress and strain. Woo et al. used finite element analysis method to determine the nonlinear stress-strain relationship for clamped human corneas under inflation test of the form of equation (6), with α= 5400 Pa and β = 28.0 [44]. The corneal motion model is shown in the following equation, where f c is the nonlinear spring resistance force which is multiplication of the spring stress and the area of the applanation at the back corneal surface, The equation of corneal motion under the resistance force from the nonlinear spring can be obtained by substituting equation (6) and (8) into equation (7) and becomes, The equivalent corneal stiffness (k e ) with respect to Goldmann Applanation Tonometer (GAT) of the nonlinear spring can be calculated using equation (10) and is used to compare with the solution from linear viscoelastic model.

Estimation of corneal parameters from clinical data
A corneal dynamics deformation response (Fig. 4) are extracted from Koprowski et al. [31] and the extracted pixel information is then converted to length with the value of 18.6 µm/pixel [37]. The corneal apex dynamic deformation profile is used to compare with the deformation response generated from the linear and nonlinear viscoelastic model, and to estimate the corresponding corneal biomechanical parameters. 29 healthy eyes of 20 normal people were included in this prospective and observational study. Central corneal thickness (CCT), average IOP, corneal radius of curvature of the patients were measured. Global visual field parameters, including median deviation (MD) and standard deviation (SD), were recorded. Patient demographics and ocular characteristics were summarized in Table 1. 20 normal population were recruited with a mean age of 57.95±8.25 years, lower quartile 53.5 years, upper quartile 61 years, and range 41 to 77 years. There was a greater proportion of females 55% (11) than males 45% (9).

Linear viscoelastic model
The corneal dynamic response y(t) can be obtained by solving equation (1) as shown in Fig. 5. The comparison between the corneal deformations y(t) and different corneal elastic stiffnesses (k), viscous damping coefficients (c) and corneal mass (m) are showed in Fig. 6-10, respectively.

Nonlinear viscoelastic model
The corneal dynamic responses y(t) obtained by linear viscoelastic model and nonlinear viscoelastic model are compared in Fig. 11. The comparison between the corneal deformations y(t) and different corneal material nonlinearity coefficients (α ± 0.2α and β ± 0.2β) and viscous damping coefficients (c) are shown in Fig. 12-16 respectively. The corneal dynamics deformation response extracted from Koprowski et al. [31] is compared with the deformation response generated from the linear and nonlinear viscoelastic model as shown in Fig. 17 and the corresponding estimated values of the corneal biomechanical parameters are shown in Table 2.

Linear viscoelastic model
The corneal dynamic response y(t) to the external air puff showed a small corneal deformation at the earlier deformation stage and showed a quite linear portion at the intermediate loading stage (Fig. 5). The corneal deformation response (Fig.  6) showed no differences with different corneal masses (m ± 0.2m; 15 ± 6 mg). The influence of mass in the air puff tonometry to the corneal deformation can be negligible and the result is consistent with the literature [38]. The magnitude of the maximum corneal deformation (ymax) in the Fig. 7 showed a significant decrease of 64.2% ± 2.4% (average for different value of viscous damping coefficients (c = 0, 0.05, 0.1, 0.15, 0.2 and 0.25 Ns/m)) with the increase of corneal elastic stiffnesses from k = 65 to k = 194 N/m. Fig. 8 showed that in the case of k = 65 N/m and c = 0.25 Ns/m, the corneal apex deformation cannot return back to the original position within 30 ms, and is consisted with the clinical observation that Keratoconic cornea with lower elasticity and higher viscosity take more time to fully recover to the original shape (position) [34,41,43]. The magnitude of the maximum corneal deformation (ymax) in the Fig. 9 showed a slightly decrease of -9.4% ± 6.2% (

Nonlinear viscoelastic model
It is shown in Fig. 11 that the corneal deformations are consistent between the linear and nonlinear viscoelastic model at the early deformation stage, however, it showed a significant divergence in the intermediate deformation stage. The spring in the nonlinear viscoelastic model is nonlinear and is subjected to strain stiffening as shown in equation (10). The higher the corneal deformation, the stiffer the cornea is. In the intermediate corneal deformation stage, the cornea is strain stiffened and the corneal stiffness become higher compared to the early stage, therefore, the cornea can withstand a larger magnitude of external load and hence, resulted in a relatively lower corneal deformation compared to the linear viscoelastic model. In the last stage of deformation, the cornea returned back and strain reduced, the cornea is softened and resumed to a smaller corneal stiffness value compared to the linear viscoelastic model, and therefore, the last stage of corneal deformations are consistent between the linear and nonlinear viscoelastic model. The magnitude of the maximum corneal deformation (ymax) in the Fig. 12-15 showed a significant decrease of 17.

Estimation of corneal parameters from clinical data
The comparison results of the corneal dynamics deformation response in Fig. 17 showed that the linear and nonlinear viscoelastic model can be used to fit the clinical corneal deformation response and to estimate the corneal biomechanical parameters. The corneal biomechanical parameters in Table 3 showed that the estimated corneal stiffness (k) and damping factor (c) are dependent on the model. The GAT equivalent corneal stiffness in the linear model (155 N/m) showed a larger value compared to the nonlinear model (66 N/m). The reason for the difference is that, the corneal stiffness here in this study is a GAT equivalent stiffness value. The corneal stiffness is a strain dependent value that it changes with strain (corneal deformation). In the GAT equivalent strain (deformation) stage (equation (4)), due to the nonlinear stress-strain behavior of the spring, the stiffness is lower at the early stage and increases exponentially as shown in Fig. 18 according to equation (10) at the intermediate deformation stage. Therefore, the estimated GAT equivalent corneal stiffness in the nonlinear viscoelastic model is low compared to linear viscoelastic model as shown in Fig. 14. The different in the estimated damping factor between linear and nonlinear viscoelastic model is due to the compensation of the strain stiffening corneal stiffness in the nonlinear model. Since the cornea is an viscoelastic material with nonlinear stress-strain relationship [44], nonlinear viscoelastic model should be used in order to have a better modeling of the corneal dynamics deformation behavior under an air puff excitation. A multi-objective genetic algorithm-based optimization procedure was used to identify the parameters of the nonlinear model, also the minimum difference between actual and optimized corneal deformation was identified. The procedure was implemented using MATLAB's global optimization. The parameters of the nonlinear model were the scaling factor (α) and the nonlinear component between stress and strain (β), viscous damping coefficient (c). Statistical analyses of the nonlinear model parameters, corneal elastic stiffness and the maximum corneal deformation were presented in Table 3.
A comparison of the clinical corneal deformation response and an optimized nonlinear viscoelastic models' corneal deformation is shown in Fig. 19. It showed minimized difference between clinical and modelled corneal deformation and that the nonlinear viscoelastic model can be used to determine the corneal biomechanical parameters.

Conclusions
In this study, viscoelastic models are presented to model the corneal dynamic deformation under the air puff excitation.

Consent for publication Not applicable
Availability of data and materials The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.