purpose. Until now, there has been no comprehensive mathematical model of the nonlinear viscoelastic stress-strain behavior of extraocular muscles (EOMs). The present study describes, with the use of a quasilinear viscoelastic (QLV) model, the nonlinear, history-dependent viscoelastic properties and elastic stress-strain relationship of EOMs.

methods. Six oculorotary EOMs were obtained fresh from a local abattoir. Longitudinally oriented specimens were taken from different regions of the EOMs and subjected to uniaxial tensile, relaxation, and cyclic loading testing with the use of an automated load cell under temperature and humidity control. Twelve samples were subjected to uniaxial tensile loading with 1.7%/s strain rate until failure. Sixteen specimens were subjected to relaxation studies over 1500 seconds. Cyclic loading was performed to validate predictions of the QLV model characterized from uniaxial tensile loading and relaxation data.

results. Uniform and highly repeatable stress-strain behavior was observed for 12 specimens extracted from various regions of all EOMs. Results from 16 different relaxation trials illustrated that most stress relaxation occurred during the first 30 to 60 seconds for 30% extension. Elastic and reduced relaxation functions were fit to the data, from which a QLV model was assembled and compared with cyclic loading data. Predictions of the QLV model agreed with observed peak cyclic loading stress values to within 8% for all specimens and conditions.

conclusions. Close agreement between the QLV model and the relaxation and cyclic loading data validates model quantification of EOM mechanical properties and will permit the development of accurate overall models of mechanics of ocular motility and strabismus.

^{ 1 }Collins et al.,

^{ 2 }and Simonsz

^{ 3 }have made measurements using strain gauges attached to isolated human EOMs during strabismus surgery that related uniaxial force and length to characterize EOM mechanical properties. Such studies characterized EOMs as if they were one-dimensional linear springs, providing results suitable for simple biomechanical models making the same assumption. In fact, however, EOMs have areas of variable thickness along their lengths, and appreciable variations in the properties as functions of time and length. One prominent mechanical feature of EOMs is

*hysteresis*, a history-dependent variation in mechanical behavior that, for example, might differ for loading versus unloading.

^{ 4 }

^{ 5 }Velocity-dependent variation mechanical behavior is termed

*viscosity*, a prominent feature of EOMs that results in energy dissipation in orbital tissues. Substantial hysteresis that is obvious in classical EOM length-tension data has nevertheless often been neglected in descriptions of even one-dimensional EOM mechanical properties.

^{ 3 }

^{ 6 }In some finite element analysis (FEA) studies, such as those by Schutte et al.,

^{ 6 }only static properties of tissues underlying ocular motility were considered despite the highly dynamic or time-dependent nature of deformation experienced by orbital tissue. Until now, there has been no effort to develop a constitutive relation that would predict passive mechanical properties of a generalized EOM with arbitrary three-dimensional size and shape.

^{ 7 }Nonlinear viscoelastic properties are significant in physiological functions of biosolids. In the field of biomechanics, theoretical models have been developed to describe various tissues, including artery,

^{ 8 }brain,

^{ 9 }

^{ 10 }

^{ 11 }

^{ 12 }

^{ 13 }

^{ 14 }heart muscle,

^{ 15 }

^{ 16 }

^{ 17 }kidney,

^{ 18 }

^{ 19 }

^{ 20 }liver,

^{ 21 }

^{ 22 }

^{ 23 }ligament and tendon,

^{ 24 }

^{ 25 }

^{ 26 }skin,

^{ 27 }

^{ 28 }

^{ 29 }and spinal cord.

^{ 30 }

^{ 31 }Despite differences among tissues in material properties and structures, gross mechanical behavior is generally influenced by mutually interacting structural components such as collagen.

^{ 23 }

^{ 26 }: continuum, microstructural, and rheological. Several theories of continuum models have been proposed. For example, rather than discretizing a system, DeHoff

^{ 32 }adapted a continuum mechanics-based constitutive equation to describe soft tissues. Unlike continuum models, microstructural models generalize the known mechanical responses of the components, often from discretization of the solid being investigated, to describe the gross tissue mechanical behavior.

^{ 27 }

^{ 28 }

^{ 33 }

^{ 34 }

^{ 35 }

^{ 36 }

^{ 37 }A main advantage of structural modeling is that it yields parameters connected with known tissue morphology, thereby relating structure to material properties. Rheological models use the simplest possible terms to describe gross tissue mechanical behavior.

^{ 23 }For instance, Viidik

^{ 38 }proposed a rheological model of viscoelastic tissues characterized by combinations of springs and dashpots. Barbenel et al.

^{ 39 }and Sanjeevi et al.

^{ 40 }also proposed rheological models exhibiting viscoelastic behavior. The most successful viscoelastic tissue model is the quasilinear viscoelasticity (QLV) theory of Fung,

^{ 7 }

^{ 23 }

^{ 41 }

^{ 42 }which can represent both linear viscoelastic and nonlinear elastic responses. The QLV models have proven helpful in describing several different tissues exhibiting viscoelastic characteristics.

^{ 16 }

^{ 43 }

^{ 44 }

^{ 45 }

^{ 46 }

^{ 47 }

^{ 48 }In biomechanical studies of articular cartilage by Woo et al.

^{ 48 }and periodontal ligament by Toms et al.,

^{ 45 }for instance, constitutive relations were well described by QLV.

^{ 49 }formulated a model having an infinite number of Voigt (parallel circuit of dampers and springs) and Maxwell (serial circuit of dampers and springs) elements. A nonlinear Kelvin model was proposed by Viidik

^{ 38 }as a sequence of springs of different natural lengths, so that the number of participating springs increases with increasing strain.

^{ 41 }

^{ 42 }to quantify the nonlinear, time-dependent stress-strain behavior under tensile loading for passive bovine EOM tissues, which are readily available in appropriate sizes and are similar in structure and function to those of humans. Although passive bovine EOM tissue lacks oxygenated perfusion and innervation, it is assumed that these deficits would have less influence on EOM mechanical properties under tensile loading because the EOM is not contracting. Incorporation of these constitutive properties into an FEM would enhance basic and clinical knowledge of EOM biomechanics.

^{ 7 }Preliminary testing determined the number of cycles and the range of strain for cyclic loading required to reach steady state. The appropriate range of displacement was obtained from simple tensile testing, in which the specimen was pulled to failure. Based on these experiments, preconditioning was chosen to be performed from stretch ratio (deformed length/original length), λ =1.0 to λ =1.1. A 10% upper limit of displacement in cyclic loading was taken to ensure that the strain remained within the linear region on the load versus deformation curve obtained from testing of a separate specimen to failure. Except for some control experiments without pre-conditioning, all specimens prepared for tensile testing, relaxation testing and cyclic loading testing were preconditioned.

*t*)/σ(

*t*∼ 250 ms)) function must be determined. Because the viscous effect is generally not evident before 250 ms

^{ 48 }and maximum strain rate achievable by the load cell was 550 mm/min, the sample was taken at maximum speed to approximately 30% elongation, which is equivalent to 92%/s strain rate. The specimen was then maintained for 1500 seconds during recording of the force required to maintain this deformation.

^{ 7 }or QLV, parsing the stress relaxation function into time-dependent and elastic portions to describe stress under constant strain as:

*t*) is the stress at any time

*t*,

*T*

^{e}(λ) is the elastic stress corresponding to a stretch ratio, λ, and

*G*(

*t*) is the reduced relaxation function representing the stress of the material divided by the stress after the initial ramp strain. In this way, the QLV model is linear with respect to relaxation but still accounts for nonlinear large deformation. The operator ∗ in equation 1represents the convolution of

*G*and

*T*

^{e}(ε), which is

*G*(

*t*) multiplied by

*T*

^{e}(ε) and then integrated over time.

*G*(

*t*) is defined as

*t*) is first Piola-Kirchhoff stress as a function of time, which relates forces in the present configuration with area in the reference (“material”) configuration. The complete stress history at any time

*t*is then the convolution integral:

*G*is the reduced relaxation function,

*t*starts at 0 rather than negative infinity. As shown in equation 5 , a decaying exponential equation is chosen as the reduced relaxation function as in Toms et al.

^{ 45 }and Wills et al.

^{ 50 }:

*a*,

*b*,

*c*,

*d*,

*g*, and

*h*are all experimentally determined. The instantaneous stress response is assumed to be represented by the nonlinear elastic relationship:

*A*and

*B*, are determined from uniaxial tensile testing, and the rest of constants

*a*,

*b*,

*c*,

*d*,

*g*, and

*h*are determined from stress relaxation testing.

*A*and

*B*were determined to be 575.1 Pa and 12.3, respectively.

*G*(

*t*) and plotted against time (Fig. 8) . As seen in Figure 8 , most stress relaxation occurred during the first 30 seconds; however, if the parameters for the reduced relaxation function are extracted from a prematurely terminated experiment, the limiting value

*G*(∞) becomes erroneous.

^{ 7 }Thus, curve fitting was performed on the entire 1500-second time series. Using the nonlinear least square method, parameters for equation 5were extracted from the average relaxation curve and are shown in Table 1 , along with constants extracted for equation 6 .

^{ 45 }the model was less accurate for the unloading region and for minimum stress values, where good agreement is not demanded for QLV validity (Figs. 9 10 11) .

^{ 45 }

^{ 48 }The QLV model based on parameters derived from nonlinear least square fitting of uniaxial strain and relaxation data predicted peak cyclic loading results within 7.4% of experimental values for multiple strain rates and ranges. The present study thus demonstrates that QLV theory is accurate in predicting peak stresses of bovine EOM specimens during cyclic loading and that the quantitative description of material properties is general. It is important to emphasize that the two components of the QLV model were derived from different types of experiments, yet their combination into the model yielded values within 7.4% of experimental data collected under still different conditions. The QLV model not only predicted the first peak stress value after application of cyclic loading, it captured the dynamic relaxation of the peaks over succeeding cycles until steady state was achieved and captured stress during loading phases. This broad range of dynamic agreement suggests that the QLV model accurately characterizes EOM mechanical properties.

^{ 24 }

^{ 26 }

^{ 48 }The present study avoided this variability by making use of a series of preliminary tests to optimize environmental conditions. Presumably, temperature control, humidity control, and precise specimen preparation reduced experimental variability. Nevertheless, slight variations remained in uniaxial tensile test and relaxation data, as seen in Figures 5 and 7 . These differences may relate to collagen fiber orientation as a function of anatomic location along the direction of loading.

*elastic response*of EOMs and the relatively simple mathematical formulations of the QLV model. Indeed, more complex mathematical formulations (such as more elaborate exponential equations, for example) for elastic and reduced relaxation functions might better characterize the cyclic unloading region.

^{ 7 }

^{ 48 }More important, observed stress values lower than predicted values suggest that the EOM relaxation occurs more quickly than predicted. This phenomenon deserves further investigation.

^{ 45 }

^{ 48 }EOMs relax occurred much earlier, at approximately 30 seconds. The duration of time at which the relaxation occurred was similar to that of the myocardium, as reported by Miller et al.

^{ 46 }This shorter relaxation was nevertheless roughly comparable to a whole orbital time constant of approximately 10 seconds, observed in the anesthetized primate using a different technique.

^{ 51 }To validate further the EOM biomechanical model, it would be desirable to derive a QLV model through the characterization of creep, the other characteristic of viscoelasticity. Given that creep is inverse to the relaxation characteristic, a creep experiment and the associated mathematical formulation based on QLV could provide an additional empiric check on the validity of a biomechanical model of EOMs. Creep testing will require additional experimental work.

**Figure 1.**

**Figure 1.**

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

a | 2.8073 |

b (1/s) | 1.57156 |

c | 0.86452 |

d (1/s) | 0.00014 |

g | 0.33755 |

h (1/s) | 0.017 |

A (Pa) | 575.1 |

B | 12.3 |

**Figure 9.**

**Figure 9.**

**Figure 10.**

**Figure 10.**

**Figure 11.**

**Figure 11.**