|Home | About | Journals | Submit | Contact Us | Français|
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.
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.
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.
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.
Because the extraocular muscles (EOMs) are manipulated mechanically during surgery for strabismus, their mechanical properties have potential therapeutic importance. Investigators such as Robinson et al.,1 Collins et al.,2 and Simonsz3 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.
Extraocular muscles should be considered to have both hyperelastic and viscous properties. Hyperelastic materials are defined as nonlinearly elastic, isotropic, and incompressible, and they have properties generally independent of strain rate. Filled elastomers and biological tissues are often modeled as hyperelastic. The viscosity of a substance reflects a strain rate dependence. Viscoelastic substances, unlike purely elastic materials, dissipate mechanical energy when subjected to cyclic loading.
Computational tools such as finite element models (FEMs) are now being more widely applied to complex anatomic structures. As the name suggests, FEA is a numerical technique for finding approximate solutions by discretization of a solid into a finite number of small elements. Given that each of the numerous elements responds to overall applied force in a manner that depends strongly on assumed constitutive tissue properties, simulated behaviors of overall FEMs of complex biomechanical structures do not substantiate the validity of the models or the validity of assumed constitutive tissue properties. Existing uniaxial length-tension measurements are insufficient descriptors of tissue constitutive properties, leaving actual biomechanical properties of orbital tissues undefined. It is, therefore, imperative to characterize constitutively the hyperelastic materials within the orbit, including EOMs and orbital connective tissue.
Most biological tissues exhibit time and history dependence of stress-strain behavior, resembling that of viscoelastic materials.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–14 heart muscle,15–17 kidney,18–20 liver,21–23 ligament and tendon,24–26 skin,27–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
Three major categories of viscoelastic models have been proposed for soft biological tissues26: continuum, microstructural, and rheological. Several theories of continuum models have been proposed. For example, rather than discretizing a system, DeHoff32 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–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, Viidik38 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–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.
Many experimental studies of soft tissue have been framed in the linear theory of viscoelasticity, relating stress and strain using Voigt, Maxwell, or Kelvin models. Buchthal and Kaiser49 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 Viidik38 as a sequence of springs of different natural lengths, so that the number of participating springs increases with increasing strain.
It is reasonable to expect that linear viscoelasticity theory should adequately describe tissue mechanical responses to small perturbations of perhaps up to 1% strain (elongation). However, much larger strains, in the range of tens of percent, are typical of the oculomotor system. For such physiological situations, the nonlinear stress-strain characteristics of the tissues must be considered. Assumptions of unrealistic tissue constitutive properties could lead to serious errors in, for example, FEMs of EOM and connective tissue relationships. Unfortunately, the classical mechanical measurements available for EOMs have only incompletely described constitutive properties and, more important, have assumed linearity in situations in which this assumption is hazardous or entirely invalid.
In the present study, we used the QLV theory of Fung41,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.
Fresh heads of adult cows were obtained from a nearby abattoir. In the laboratory, orbits were carefully dissected for extraction of EOM and connective tissue. Transport time from abattoir to laboratory was approximately 30 minutes; the additional time elapsed to dissect the EOMs averaged 3 hours Extracted EOMs as prepared for experiments are shown in Figures 1A–F for the lateral rectus (LR), medial rectus (MR), inferior rectus (IR), and superior rectus (SR) and for the inferior oblique (IO) and superior oblique (SO) muscles, respectively.
After extraction, EOMs were maintained in lactated Ringer's solution at 37°C. To minimize axial damage to EOM fibers, each specimen was initially prepared in the shape of an approximately 8- to 10-cm long prism with a 2.5 mm × 2.5 mm cross-section and a length equal to that of the EOM. For consistency, samples were prepared from the transverse center of each EOM. Specimens were then divided into 25-mm lengths for biomechanical measurements. In preliminary testing, specimens were prepared with both ends cut slightly thicker to reduce the stress concentration anticipated to cause mechanical failure around the clamped ends. However, this precaution proved unnecessary because preliminary experiments showed no differences between Lagrangian stress (force/original cross-sectional area) and stretch ratio (stretched length/original length) or location of failure (specimen breakage) between specimens with uniform thickness and specimens with the thicker ends. Therefore, all the specimens for which data are reported were prepared at 25-mm length and at uniform thickness. Given that both ends had to clamped, the actual tested length was the 10-mm middle portion of each specimen and in every case avoided the terminal tendon. Specimen preparation time was approximately 45 minutes.
Experiments were conducted using a tabletop micro-tensile load cell (model 4411; Instron, Norwood, MA) with compatible software (Series IX; Instron). In this instrument, designed for testing a wide range of materials in tension and compression, the specimen was secured between the rigid base of a load cell and a vertically moving crosshead. The crosshead was mechanically displaced upward or downward to apply a tensile or compressive load to the specimen, respectively. A force transducer mounted between the specimen and crosshead measured the applied load. In the present study, a ±50 N load transducer with 0.001 N precision (2530–437; Instron) was used to perform tensile tests. Using the general-purpose interface bus interface (model 4411; Instron), the load cell was connected to a computer for automatic test control and data logging. Before carrying out experiments, a preliminary test was performed to determine equipment sensitivity. Measured voltage change during loading of the transducer was 1 V over a range of a 5-N force. Given instrument sensitivity of 2 mV/V, actual static force sensitivity was thus determined to be 0.01 mN—adequately precise for the purpose.
Environmental control is important in testing biological specimens. Given that preliminary tests indicated that humidity and temperature influenced the result of EOM testing, the load cell was enclosed in a plastic chamber within which a humidifier and heater maintained humidity and temperature at optimal levels. A humidity/temperature sensor (Fisher Scientific, Chino, CA) indicated 93% relative humidity and temperature 37°C adjacent to the specimen.
Biomaterial testing is considered to require proper preconditioning because the internal structure of tissue changes during cyclic loading. By repeated cycling, a steady state is eventually reached after which no further change occurs unless the cycling pattern is changed.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.
To reduce the number of specimens tested and to verify intuition that all six EOMs were mechanically similar, uniaxial tensile testing was performed on 12 EOM specimens, two of each EOM type. Each specimen was elongated or imposed with tensile loading to failure (rupture) to evaluate the instantaneous stress response. Results of tensile testing not only determined the upper limit in cyclic loading for preconditioning, they characterized the elastic response function in modeling.
To characterize viscoelastic or time-dependent behavior, the reduced relaxation (σ(t)/σ(t ~250 ms)) function must be determined. Because the viscous effect is generally not evident before 250 ms48 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.
To validate the model that was characterized using a combination of relaxation and tensile test data, cyclic loading tests were carried out at multiple strain rates (1.42%/s, 2.77%/s, and 3.14%/s) and extension ranges (21.1%–34.6%, 18.1%–31.2%, 4.1%–28.9%) for three specimens. The range of extension had to be well within the linear range of force versus displacement; therefore, maximum extensions for all three tests were within 40% of overall specimen length. Cyclic loading results were then compared with model predictions.
Constitutive equations based on the conventional theory of linear elasticity assume small strain. However, most soft tissues undergo large deformation that violates this assumption. To resolve this dilemma, Fung introduced the theory of quasilinear viscoelasticity,7 or QLV, parsing the stress relaxation function into time-dependent and elastic portions to describe stress under constant strain as:
where σ(t) is the stress at any time t, Te(λ) 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 1 represents the convolution of G and Te(ε), which is G(t) multiplied by Te(ε) and then integrated over time. G(t) is defined as
We then assume that the stress response to an infinitesimal change in stretch superposed on a specimen in state of stretch λ at instant of time τ can be written as equation 3:
σ(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:
where G is the reduced relaxation function,
represents the instantaneous elastic response for infinitesimal change in stretch ratio, and
is the strain history. In experimental settings 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:
In equation 5, coefficients 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:
Thus, in this one-dimensional QLV model, two constants, 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.
When a biomaterial is subjected to cyclically varying strain, the resultant stress exhibits a cyclic hysteresis loop. A large, variable initial period of hysteresis after a large strain change is behavior common to all biomaterials. Figure 2 shows how the hysteresis loop for a sample EOM decreases with succeeding cycles, eventually reaching steady state after three to five cycles. Imposition of cyclic stress until steady state is reached is known as preconditioning. The strain rate for preconditioning was chosen as 1.667%/s because different strain rates did not significantly alter the hysteresis loop. The range of deformation was chosen to be 10% to avoid plastic (irreversible) deformation.
From preliminary tests, it was verified that preconditioning at 10% deformation would yield consistent results. Therefore, before every test, each specimen was preconditioned with five cycles of 10% elongation. In addition, to verify that the preconditioning process does not substantially alter the mechanical properties of the specimens tested, three specimens were loaded to tensile failure without preconditioning and were compared with similar but preconditioned specimens (Fig. 3). Measured stress to failure at approximately 1.8 times initial length was similar with or without preconditioning.
Lagrangian stress (force/initial undeformed cross-section area) was determined in two specimens of each EOM. This method effectively normalizes against variations in cross-section among the rectus and oblique EOMs. Measured values of Lagrangian stress agreed well both between the duplicate EOMs of the same type (for example, two MR EOMs from two different orbits, data not shown) and among all six types of EOMs generally (Fig. 4).
Based on their similarity, data from two samples of each of the six EOM types was pooled and are illustrated in Figure 5. This figure demonstrates linear variation in stress in region A-B from 30% to 60% elongation, whereas maximum load was reached at C, above which the specimens failed (tore). The slope of the line in linear region A-B can be interpreted as “elastic stiffness.” On average, the maximum values exceed the mean by 14.2%, and the average minimum value was 13.6% lower (Fig. 5). An important result from tensile testing was that all six EOM types exhibited similar response to loading. This finding represents a strong indication that all bovine EOMs have similar viscoelastic properties.
Figure 6 shows the relationship between Lagrangian stress and strain for EOMs, plotting the observed mean of 12 specimens. Using the nonlinear least squares method, parameters for equation 6 were extracted. Values for A and B were determined to be 575.1 Pa and 12.3, respectively.
Stress relaxation tests at 30% deformation were carried out for 16 different EOM samples, consisting of two LR, three MR, three IR, two SR, three IO, and three SO specimens. The data indicate that all six EOM types have a similar time-dependent nature (similar decay to asymptotic value) and elastic properties (similar peak stress value on imposition of ramp stimulus). Because the behavior of all six EOM types was similar (Fig. 7A), the data for all 16 specimens were pooled in Figure 7B. On average, maximum values exceed the mean by 11.4%, and the average minimum value was 11.4% lower.
Because the QLV model requires a reduced relaxation function, the original relaxation function was first converted to 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 5 were extracted from the average relaxation curve and are shown in Table 1, along with constants extracted for equation 6.
From consistent results of relaxation and tensile loading tests, it was evident that all six EOMs exhibited similar mechanical properties. Viscoelastic materials were defined to have similar mechanical properties when time-dependent and elastic characteristics were same. Because it was established that time-dependent nature and elastic response were similar for all bovine EOMs, it was assumed that they behaved similarly under cyclic loading, a test therefore not performed for every EOM type. Three experiments were conducted to validate the assembled QLV model. Each of three EOM specimens was tested for 10 cycles in a different cyclic loading condition. The aim of the model was to characterize accurately the dynamic change in peak cyclic stress over time. Figure 9 shows data from a trial in which the strain rate was 3.14%/s and the strain range was 4.1% to 28.9%. As seen in Figure 9, the experimental data agreed well with QLV model predictions. There was similar good agreement between predicted and observed peak stress values for separate experiments (Figs. 10 and and11)11) performed at 1.42%/s and 2.77%/s, with strain rages of 21.1% to 34.6% and 18.1% to 31.2%, respectively. The average difference between experimental and predicted peak values in all three trials over 10 cycles was 7.4%. The error in the first three peak values was even less, 3.6%. Although the accuracy was excellent in predicted peak stress, thus validating the QLV for EOMs,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–11).
A mechanical model of material properties is considered well defined when loading and peak stress values agree well with experimental values.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.
High variability in tensile testing is often reported for biosolids.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 and7.7. These differences may relate to collagen fiber orientation as a function of anatomic location along the direction of loading.
Despite accurate QLV predictions of loading phase and peak stress values in cyclic loading, predictions were less accurate during unloading phases and minimum stress values (Figs. 9–11). These disparities may be related to the rapid stretch rate needed for obtaining the experimental 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.
The extracted EOM specimens tested here were in a non-physiological state, lacking in innervation and oxygenated blood perfusion, and not subject to active tonus. Such physiological conditions are not achievable using the present apparatus. However, the goal of the present study was to formulate a model capturing the passive mechanical behavior of EOMs under tensile loading. The present model does not attempt to capture changes in mechanical properties caused by active EOM contraction that were not the focus of the study. That worthy endeavor will be facilitated by the present demonstration of the validity of QLV for EOMs.
It must be recognized that plasticity or permanent deformation might occur during EOM tensile testing. Unloading behavior could have been influenced by plasticity that might have been avoided if the experiments had been conducted in vivo. In vivo experiments in adult cows were technically impossible using the present apparatus. Even so, unloading behavior is of lesser importance that loading behavior in evaluating a QLV model, which is said to be well defined as long as predicted peak values are in good agreement with experimental data.45
The present data also indicate that the tensile stress relaxation of EOMs is rapid, as typical of a material with a short memory. As can be seen in Figure 8, unlike other biomaterials such as articular cartilage, in which most relaxation occurred in the first 15 minutes,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.
An ultimate goal of EOM biomechanical modeling is to produce a comprehensive, quantitative mechanical description of EOM viscoelastic properties and the way in which these properties change during active contraction. The present demonstration that QLV is valid for EOMs subjected to passive loading suggests an improved approach to study mechanical behavior using FEM. Once a valid constitutive relationship is implemented into an FEM model, the model can be used to study and compare the responses of normal, diseased, surgically operated, or aging EOMs under various loading conditions.
The authors thank Manning Beef, LLC (Los Angeles, CA) for the generous contribution of bovine specimens; Jose Martinez, Claudia Tamayo, and Ramiro Carlos for assistance with specimen preparation; and Andrew Shin for assistance with photography.
Supported by National Eye Institute Grants EY08313 and EY00331. JLD is Leonard Apt Professor of Ophthalmology.
Disclosure: L. Yoo, Manning's Beef, LLC (F); H. Kim, Manning's Beef, LLC (F); V. Gupta, Manning's Beef, LLC (F); J.L. Demer, Manning's Beef, LLC (F)