## Abstract

Experimental data on the passive mechanical properties of the ventral interior lateral muscle of the tobacco hornworm caterpillar, *Manduca sexta*, are reported. The stress–deformation response of the *Manduca* muscle is shown to be nonlinear pseudo-elastic, capable of large deformations and subject to stress softening during initial loading cycles. The muscle passive mechanical properties also depend on multiple time-dependent processes. In particular, we show new experimental data from cyclic loading tests of an unstimulated muscle with constant maximum stretch and different, constant engineering strain rates. Then, on the basis of these data a constitutive model is derived to reproduce the main characteristics of this behaviour. In formulating the constitutive model, we consider the muscle as a complex macromolecular structure with fibrous components at numerous size scales. The model uses a phenomenological approach to account for different mechanisms by which passive force changes during applied deformation and how the muscle properties recover after unloading.

## 1. Introduction

Animals moving in natural environments must cope with complex three-dimensional surroundings that change over different time scales and are never identical at all locations. The motor systems that have evolved to adapt to such variability are diverse in overall structure but remarkably similar at the level of tissues and molecular mechanisms. All animals use neural commands to direct the contractile state of muscles, which are themselves composed of interacting proteins. Cross-bridge cycling between actin and myosin is the fundamental building block of muscle contraction, but each type of muscle is tuned to specific applications by other molecular and mechanical processes. Hence, cardiac muscle is distinguished from skeletal muscle by cross-striations, unique ion channels and enzyme isoforms and by its myogenic activity and force–length characteristics. Similarly, there are significant differences between striated and smooth muscles and between adult insect muscles and those of the larval stages.

Most studies focus on the *active* properties of these tissues; that is, their role as motors that are activated to develop tension or shortened to perform work. However, there is a growing appreciation of the importance of the *passive* properties of muscles together with their roles as struts, beams, brakes and dampers. Unlike either amorphous or crystalline materials, muscles are complex composites with profound anisotropy. Each muscle fibre contains aligned actin and myosin filaments within an amorphous matrix material composed of proteins, lipids and polysaccharides. In addition, muscles are now known to contain molecularly massive elastic proteins such as titin (Ziegler 1994; Tskhovrebova *et al*. 1997; Tskhovrebova & Trinick 2002, 2003; Granzier & Labeit 2004), with folded domains that result in asymmetric, time-dependent stretching and shortening during load cycling. The contribution of these material responses to muscle function and animal behaviour is not well understood. This is partly because most commonly studied muscles function in combination with tissues (tendons, ligaments, bones and joints; Marsh 1999; Biewener 2003) that have very different material properties. In these cases, characterization of the material properties of isolated muscles is insufficient to describe *in vivo* function. In addition, it is difficult to interpret the responses of many skeletal muscles that are composed of both fast and slow cell types with differing physiological and biochemical properties.

In response to this limitation, the approach taken in the current work is to characterize and model the properties of muscles in a soft-bodied animal, the caterpillar *Manduca sexta*. In contrast to most vertebrate muscles, or even those of most adult arthropods, *Manduca* muscles do not act upon large attachment structures or jointed skeletons. Instead, they must function much like tongues, trunks and vertebrate hearts (Kier & Smith 1985). The muscles are attached to the body wall such that muscle activity has a direct action on overt movements (Belanger & Trimmer 2000; Mezoff *et al*. 2004; Dorfmann *et al*. 2007; Trimmer & Issberner 2007). Furthermore, because *Manduca* changes shape, the muscles undergo very large deformations (Woods *et al*. 2007) and their passive properties are expected to be particularly important. Structurally and biochemically, larval *Manduca* muscle cells appear to be of one type and there is no differentiation into fast and slow fibres.

In a previous work (Dorfmann *et al*. 2007), the force-extension characteristics of *Manduca* muscles were measured under constant rate loading and unloading, both in a passive state and during tetanic stimulation. *Manduca* muscles were found to be elastomeric in both states, dissipating work with each strain cycle (pseudo-elasticity) and undergoing stress softening on repeated cycling (the Mullins effect). These properties have been described in detail for natural rubber reinforced with carbon black filler particles (Dorfmann & Ogden 2004), and they can be modelled using constitutive equations based on nonlinear continuum mechanics principles (Fung 1993). This approach contrasts with the more common mathematical descriptions of muscle that rely on mechanical theories developed for small deformations occurring in relatively rigid materials. In their simplest form, these models are variations of Voigt and Maxwell dampened spring systems or more sophisticated versions that incorporate viscoelasticity (Speich *et al*. 2006). Other models are based on known and postulated biochemical (Holmes & Geeves 2000; Huxley 2000) or structural processes with different time and strain dependencies (Campbell & Lakie 1998). The advantage in developing constitutive equations based on nonlinear solid mechanics theories is that they capture the characteristics of a material without restrictive assumptions on its behaviour and treat muscle as an anisotropic three-dimensional body capable of large deformations. It is therefore possible to model muscle during biologically realistic elongations and forces.

An outstanding issue in passive muscle properties is the mechanism by which passive force changes during applied strain and how it recovers after unloading (Proske & Morgan 1999; Mutungi & Ranatunga 2000). This behaviour includes the thixotropic short-range elastic component (Hill 1968) of frog's striated and smooth muscles. Both actin–myosin cross-bridge breakage and reformation, and the unfolding of gap-filament proteins (e.g. titin) have been proposed as probable mechanisms (Granzier & Wang 1993*a*,*b*; Campbell & Lakie 1998; Bagni *et al*. 2004; Granzier & Labeit 2004). Intramuscular collagenous structural elements (Gosline *et al*. 2002) and muscle junctions (Lieber *et al*. 2000) may also contribute to the properties observed in experimental muscle preparations. Complicating the matter still further is the possible variation in the relative contribution of different structural elements in different locomotory muscles, even in the same species (Prado *et al*. 2005). These different molecular mechanisms are important because they influence the assumptions and validity of most mechanical models used to describe muscle behaviour. In contrast, the constitutive approach described here describes changes in the effective contour length of fibrous molecules but makes no assumptions about the mechanism underlying this process (see §4.1). One extremely important aspect of muscle function is that its performance often depends on multiple time-dependent processes. It is difficult to incorporate these changes accurately using classical mechanical models but is theoretically attainable using constitutive approaches. In the work described in the current paper, we have taken the first step towards this goal by measuring and modelling the responses of passive muscle to constant rate deformations at differing strain rates. We show that, in addition to strain-dependent stress softening, the responses of *Manduca* muscle depend on the rate of loading and unloading. In preconditioned muscles, higher strain cycling rates result in higher stiffness and increased nonlinearity of the work loops. A constitutive model of *Manduca* muscle is described, which accurately reproduces the hysteresis and rate dependency of these muscles. This model is expected to be useful in the construction of numerical simulations of soft-bodied locomotion.

## 2. Mechanical behaviour of *Manduca* muscle

### 2.1 Experimental methods

In the absence of a stiff skeleton, *Manduca* larval muscles are attached directly to apodemes, infoldings of the cuticle wall. The muscles are organized in repeated segments corresponding to the body segments, in layers. Test data were gathered for the ventral interior lateral (VIL) muscle of the third abdominal segment (A3). A3 VIL is one of the largest larval muscles, comprising 14 fibres (Dorfmann *et al*. 2007), the most interior of the layer of muscles spanning the proximal terminus of the first pair of prolegs. A3 VIL was chosen because its rate-independent pseudo-elastic properties have been characterized and modelled, both with and without stimulation (Dorfmann *et al*. 2007), and because during crawling VIL reaches its highest stress during the lengthening portion of its strain cycle (Woods *et al*. in press), suggesting that its passive properties play an important role in its biological function.

Experimental methods have been previously reported in detail by Dorfmann *et al*. (2007). Briefly, muscle length was determined prior to the dissection of a cold-anaesthetized larva in the second day of its fifth instar using previously established external markers of muscle attachment points. The muscle was dissected out in physiological saline (Trimmer & Weeks 1989) along with a small portion of attached cuticle at each end. The preparation was transferred to a horizontal saline bath, with one end pinned to an elastomer platform in the bath and the other secured by a hook to an Aurora 300B-LR lever arm ergometer that administered strain cycling while measuring force (Aurora Scientific, Inc., Aurora, Ontario). The departures in the current experimental procedure from methods described by Dorfmann *et al*. (2007) are that (i) values of engineering strain rates 0.0144, 0.072, 0.36 and 1.8 s^{−1} are imposed and (ii) all measurements are of unstimulated muscle.

A3 VIL is capable of large nonlinear elastic deformations and shows both hysteresis and stress softening during repeated strain cycling in both the passive and active conditions (Dorfmann *et al*. 2007). For active muscle, part of the stress softening is due to the time course of the contractile force developed by the muscle during stimulation. Owing to the diminishing of contractile force over time, data for stimulated muscle at lower strain rates, which would require very lengthy stimulation, are not biologically representative or experimentally feasible.

During horizontal crawling, A3 VIL experiences strain rates between 0.2 and 0.3 s^{−1}. During exploratory movements, however, strain rates can be far lower, while during defensive *strike* behaviour (Walters *et al*. 2001) strain rates can be far higher, with 2 in some instances (Michael Simon 2007, unpublished data).

### 2.2 Rate-dependent stress–deformation response

To assess the effect of rate dependence on the mechanical response of an unstimulated A3 VIL muscle, several series of periodic loading–unloading uniaxial extension tests were carried out at a constant temperature of 25°C.

The periodic loading, unloading and reloading tests were performed using constant strain rates representative of those encountered in nature. The initial distance of the pinned connections at each end of the muscle was used to determine the longitudinal strain. The reference length of the muscle was found to be 5.5 mm. Changes in the distance between these connections were measured with an accuracy of 1 μm. The tensile force was measured with an accuracy of less than 0.3 mN. The nominal stress was determined as the ratio of the axial force to the cross-sectional area, measured in the reference configuration. For a detailed description of the microscopic analysis to determine the cross-sectional area, we refer to Dorfmann *et al*. (2007). The magnitude of the cross-sectional area of the muscle tested and used to determine the nominal stress was 0.4 mm^{2}.

During each of the tests, the muscle was subjected to five cycles of preconditioning, with constant strain rate up to a pre-selected extension with stretch *λ*=1.18. The preconditioning was performed in order to monitor the progression of stress softening and to determine the ultimate stress–deformation response for stretches up to *λ*=1.18. Figure 1 shows the nominal stress versus stretch *λ* for the muscle in unstimulated state, with engineering strain rates of 0.0144, 0.072, 0.036 and 1.8 s^{−1}. Figure 2 compares the response of the muscle during the fifth loading–unloading cycle for the different loading rates. The following observations are made.

The initial configuration, corresponding to the natural configuration of the muscle in the animal, is not stress free. We take the corresponding geometric configuration to be the reference configuration from which to measure any subsequent extension generated by the application of mechanical loads.

The stresses corresponding to the same strain level increase with the loading rate during loading and unloading.

There are large differences in the stresses corresponding to the same strain level under loading and unloading during the first cycle in periodic tests with a fixed strain amplitude. The differences increase with the loading rate.

There is a reduction in the stress at a given strain on each successive loading. The reduction is largest on the first and second loading–unloading cycles and becomes rather small after approximately five cycles. Again, the reduction in the stress increases with the loading rate.

After five preconditioning loading–unloading cycles, the stress–stretch responses are essentially repeatable and additional stress softening is negligible.

Energy lost during cyclic loading is often presented as an

*efficiency value*or*work loop efficiency*. Using data from the fifth loading–unloading cycle to evaluate the work loop efficiency, no dependence on the loading rate is observed. The efficiency values for cyclic loading with engineering strain rates of 0.0144, 0.072, 0.36 and 1.8 s^{−1}are 0.0025, 0.0029, 0.0027 and 0.0027 Nmm, respectively.

We note that the degree of stress softening and the ultimate stress–deformation responses are dependent on time and electrical stimulus (Dorfmann *et al*. 2007). However, in this study we are not concerned with phenomena such as muscle stimulation, recovery time and viscous effects.

## 3. Constitutive modelling

The mechanical behaviour of hard tissues, such as bones and teeth, can be described using infinitesimal deformation theories. Under this condition, no distinction is made between the current deformed configuration and the undeformed natural reference state. In contrast, soft tissue can sustain finite deformations and therefore nonlinear theories must be used to describe the structural response. Existing finite deformation models may be distinguished based on a micromechanical or phenomenological approach. Micromechanical models (Cowin (2004) used the term mechanistic models) are less developed, rely on biological mechanisms, provide insight into chemical processes and interaction and may also describe the rearrangements of the microstructure. Phenomenological models, on the other hand, provide a descriptive and predictive capability for the mechanical response of soft tissue. They are based on a continuum mechanics approach and, in the simplest situation, are capable of simulating the nonlinear elastic behaviour of isotropic materials. The theory has been extended to incorporate more complicated effects such as pseudo-elasticity, anisotropy and rate-dependent stress–strain responses. We will first provide an overview of the basic equations of the finite deformation theory and then develop a constitutive phenomenological model for the rate-dependent stress–strain behaviour of the *Manduca* muscle.

### 3.1 Basic equations

Consider a body as a set of points referred to as particles (or material points), which have a one-to-one correspondence with points of a region in a three-dimensional Euclidean point space. We denote such a region as a *configuration* of the body. In the reference configuration, the body occupies the region _{r} and a generic point of is identified by the position vector ** X** relative to an arbitrary chosen reference configuration. For biological tissue, in general, the reference configuration

_{r}will not be stress free as a result of non-uniform tissue growth that has taken place up to this time, i.e. there are residual stresses in

_{r}. These residual stresses have an influence on the overall mechanical and biological interactions as widely reported during the past two decades. We take the fixed geometric configuration

_{r}to be the reference configuration from which to measure any subsequent deformations generated by the application of mechanical loads. Suppose the application of mechanical loads deform the body so that the point

**occupies the new position**

*X***=**

*x***(**

*Χ***) in the time-independent deformed configuration, which we denote by . The vector field**

*X***, which is a one-to-one, orientation-preserving mapping with suitable regularity properties, describes the deformation of the body.**

*Χ*The deformation gradient tensor relative to _{r} and its determinant are(3.1)where Grad denotes the gradient operator with respect to ** X** and wherein the notation

*J*is defined. The Cartesian components of are given by , where

*x*

_{i}and

*X*

_{α}are the Cartesian components of

**and**

*x***and**

*X**i*,

*α*=1, 2, 3. Denoting d

*V*and d

*v*as volume elements in

_{r}and , respectively, we have the relation(3.2)and for a volume-preserving (isochoric) deformation,

*J*=det =1.

The unique polar decompositions of the deformation gradient are(3.3)where is a proper orthogonal tensor and and are the positive definite and symmetric, respectively, the right and left stretch tensors. It follows that det =det =det . The tensors and have the spectral forms(3.4)where *λ*_{i}>0, *i*=1, 2, 3 are the principal stretches; the unit vectors *u*^{(i)} and *v*^{(i)} indicate the principal directions in the reference and current configurations, respectively; and ⊗ denotes the tensor product. For incompressible materials,(3.5)Using the polar decompositions (3.3), we may define the following measures of deformations:(3.6)which are, respectively, the right and left Cauchy–Green deformation tensors. The three principal invariants of , equivalently , are defined by(3.7)where tr is the trace of a second-order tensor.

The motion ** Χ**, with =Grad

**, deforms an infinitesimal material line element d**

*Χ***at**

*X***into a spatial line element d**

*X***(which consists of the same material as d**

*x***) at**

*X***. Using equations (3.3), we have(3.8)**

*x*For a detailed discussion on the kinematics of continua, we refer to, for example, Ogden (1997) and Holzapfel (2001).

### 3.2 Isotropic hyperelasticity

A hyperelastic material is defined as an elastic material whose nominal stress is given by(3.9)where the scalar function *W*(), called the *strain energy function* (or *stored energy function*), is a function of only. The strain energy function represents the work done per unit volume at ** X** in changing the deformation gradient from to . The Cauchy (or true) stress tensor has the from(3.10)

The elastic stored energy is a scalar-valued function and therefore indifferent to a superposed rigid body rotation after deformation, i.e. it is required to be objective. Therefore, for all rotations and each deformation gradient , we have(3.11)Using the polar decomposition (3.3) and the choice =^{T} in (3.11) gives(3.12)Thus, *W* depends on only through the stretch tensor or equivalently through the Cauchy–Green deformation tensor , where we recall that =^{2}.

We now restrict attention to isotropic materials, for which the mechanical response is, by definition, unaffected by a rigid body rotation in the reference configuration prior to loading. Then, for all rotations we have the condition(3.13)Bearing in mind that the 's appearing in (3.11) and (3.13) are independent, the combination of these two equations yields(3.14)for all rotations . Equation (3.14) states that *W* is an isotropic function of . It follows from the spectral decomposition (3.4) that *W* depends on only through the principal stretches *λ*_{1}, *λ*_{2}, *λ*_{3}. Equivalently, we may regard *W* as a function of the principal invariants of (3.15)

### 3.3 Incompressibility

Incompressible (constrained) material behaviour is a convenient idealization and is used frequently in practical applications of biological materials. As a matter of fact, many biological tissues do not change their volume for deformations within the physiological range. The incompressibility constraint is given by equation (3.5), equivalently *I*_{3}≡1. The expressions (3.9) and (3.10) for the nominal stress and Cauchy stress are modified and are now given, respectively, by(3.16)where *p* is a Lagrange multiplier associated with the incompressibility constraint (3.5) and represents an arbitrary hydrostatic pressure. Equation (3.16) states that the stress in an incompressible material is determined by only within an arbitrary hydrostatic pressure. The dependence of the strain energy function *W* on *I*_{3} can be removed, hence we write *W*=*W*(*I*_{1}, *I*_{2}).

### 3.4 Transversely isotropic materials

A transversely isotropic material is composed of an elastic matrix and an embedded family of fibres, attributing the composite strong directional properties. In general, the stiffness is much higher along the fibre direction (preferred direction) compared with directions orthogonal to the fibres. The preferred direction in the reference configuration _{r} is given by the unit vector ** M** and we note that the material response is unaffected by an arbitrary rotation around

**and by replacing**

*M***by −**

*M***. For transversely isotropic materials, the stress at a material point**

*M***depends in addition to the deformation gradient on the fibre direction**

*X***.**

*M*The material considered here is said to be transversely isotropic if the strain energy function *W* is an isotropic function of the two tensors and ** M**⊗

**. Following the work by Spencer (1972, 1984) and Ogden (2001), the form of**

*M**W*is reduced to dependence on the principal invariants

*I*

_{1},

*I*

_{2},

*I*

_{3}of the right Cauchy–Green deformation tensor as defined by equation (3.7) and the additional invariants

*I*

_{4}and

*I*

_{5}given by(3.17)Note that the invariant

*I*

_{4}represents the square of the stretch in the fibre direction

**. This implies that the strain energy for an unconstrained transversely isotropic material may be written as(3.18)and for an incompressible material subject to the constraint (3.5), equivalently**

*M**I*

_{3}≡1,

*W*is regarded as a function of the remaining four invariants (

*I*

_{1},

*I*

_{2},

*I*

_{4},

*I*

_{5}).

### 3.5 Constitutive equations

Restricting attention to incompressible materials with *I*_{3}≡1, the strain energy of a transversely isotropic material is a function of the four invariants (*I*_{1}, *I*_{2}, *I*_{4}, *I*_{5}). The response of the constrained transversely isotropic material, with fibre direction ** M** in the reference configuration

_{r}, is given by the nominal stress and the Cauchy stress shown in equation (3.16). To write the explicit forms of and , we need formulae for the derivatives of the invariants with respect to the deformation gradient . These are(3.19)(3.20)

A direct calculation then leads to(3.21)(3.22)where , *i*=1, 2, 4, 5. When the dependence on *I*_{4} and *I*_{5} in equations (3.21) and (3.22) is omitted, the associated expressions for an isotropic material are obtained.

For fibre-reinforced materials, it is common to write the strain energy function as the sum of two terms, one associated with the isotropic properties of the base matrix and the other to introduce transverse isotropy in the mechanical response. We follow this tradition and consider a strain energy function given by(3.23)where the term *W*_{iso} represents the isotropic matrix material and *W*_{fib} accounts for the directional reinforcement, the latter also known as the *reinforcing model* (Qiu & Pence 1997; Merodio & Ogden 2003).

### 3.6 Pseudo-elasticity

The notation of *pseudo-elasticity* has been introduced by Fung (1980) to describe the hysteretic response of biological tissues during cyclic loading, i.e. the loading–unloading responses do not coincide, even though the body may return to its original state. Pseudo-elasticity describes an elastic material during loading and a different elastic material during unloading. A theory describing the mechanical response of these materials has been developed by Ogden & Roxburgh (1999) and used by Dorfmann & Ogden (2003, 2004) to account for stress softening and residual strains in carbon black-reinforced elastomers.

Following Dorfmann *et al*. (2007), we use the term *pseudo-elasticity* to describe the hysteretic response of a preconditioned *Manduca* muscle during a loading–unloading cycle (see experimental data in §2.2). The theory, developed by Ogden & Roxburgh (1999), modifies the elastic strain energy function *W*() by incorporating an additional variable *η*. Thus, we write(3.24)The inclusion of *η* provides a means of continuously changing the form of the strain energy function and as a consequence the stress–strain response during the deformation process. In general, the response of the material is then no longer elastic and *W*(, *η*) is referred to as a *pseudo-energy* function. In this section we provide an overview of the main ingredients of the general theory of pseudo-elasticity and appropriate specifications will be made in §4.2.

The internal variable *η* may be inactive or active; activating *η* introduces a change in the material properties. A change from inactive to active may be induced, for example, when unloading is initiated.

If the variable *η* is inactive, we set it to the constant value unity and write(3.25)for the resulting strain energy function. In (3.25) and in what follows the subscript ‘0’ is associated with the situation in which the variable *η* is inactive. For an incompressible material, the nominal stress associated with an inactive *η*, denoted _{0}, is given by(3.26)

If *η* is active, we take it to depend on the deformation through the gradient . The nominal stress is then given by(3.27)Following Ogden & Roxburgh (1999), we take *η* to be given implicitly by the constraint(3.28)which uniquely defines *η* in terms of . We may write the solution to equation (3.28) formally as(3.29)Then, the expression of the nominal stress (3.27) simplifies to(3.30)whether or not *η* is active, where when *η* is active the right-hand side is evaluated for *η* given by (3.29). It is convenient to introduce the notation *w* for the resulting (unique) strain energy function. Thus,(3.31)and the nominal and Cauchy stress tensors for incompressible materials are given by the standard relations(3.32)This general framework allows for substantial flexibility with regard to the dependence of *W* on *η* and on the form of the function *η*_{e}() in (3.29).

## 4. A model for the *Manduca* muscle

In this section we develop the constitutive model to describe the nonlinear mechanical response of a preconditioned *Manduca* muscle in the passive state. The muscle consists of a microfibrillar matrix material with an embedded fibrous network of elastic proteins, and the proposed model accounts for its nonlinear mechanical properties, including finite deformation pseudo-elasticity, anisotropy, energy dissipation associated with hysteresis and rate-dependent stress–strain response. We further assume incompressibility, since changes in volume for deformations within the physiological range are small and therefore negligible. The model is fitted to available data and its predictions are assessed.

It should be emphasized that in the natural configuration the muscle is not stress free and if released from the body wall its length will shorten. In this paper, for simplicity, our model does not account for the residual stresses in the natural configuration. We assume instead that our reference state, from which the deformation is measured, coincides with the geometric configuration of the preconditioned muscle in the passive state where no residual stresses remain (see experimental data in §2).

For an elastic, incompressible and transversely isotropic material, the total strain energy *W*, corresponding to equation (3.23), derives contributions from the isotropic microfibrillar matrix and the fibrous network. For the particular application here considered and in order to reduce the number of material parameters, it is convenient to eliminate the dependence of *W* on the invariants *I*_{2} and *I*_{5}. We consider a reduced form of equation (3.23) given by(4.1)which still provides sufficient flexibility to capture the experimentally observed response of the *Manduca* muscle.

### 4.1 Fibrous network behaviour

Muscles are complex macromolecular structures with fibrous components at numerous size scales. Thick and thin filaments composed mainly of myosin and actin, respectively, form the primary fibrous components. In striated muscles, these are arranged in more or less parallel alternating strands oriented along the primary axis of the muscle fibres, with each thick filament surrounded by six to eight thin filaments. They are linked to one another by reversible cross-bridges (which drive the passive-to-active state transition) and attached to the complex proteinaceous Z-line (or disc) that spans the short axis of each muscle fibre. In vertebrates, and most adult insect skeletal muscles, these Z-lines are aligned into contractile elements called sarcomeres and bundles of muscle fibres are also aligned to give the *striated* appearance. During activation, the actin and myosin filaments slide past one another and the Z-lines constitute an ultimate barrier limiting the maximum effective shortening. Massive elastic proteins of the titin family connect the thick filaments to the Z-line and are proposed to contribute significantly to the passive properties of resting muscle. Titin varies in size depending on the muscle type and appears to be folded to different degrees at slack sarcomere length. For example, vertebrate cardiac muscle expresses the shortest titin isoform and it is packed between the thick filament and the Z-line, where it acts as a bidirectional spring (Granzier & Labeit 2004). In contrast, much longer titin isoforms with repeating domains are found in soleus muscle and titin is thought to straighten and then unfold as the muscle is stretched, thereby providing strain- and rate-dependent spring-like behaviour (Tskhovrebova *et al*. 1997; Tskhovrebova & Trinick 2003). Force-extension experiments on single molecules generate a sawtooth-like pattern of equally spaced peaks. Every downstroke of the response reflects the mechanically induced unfolding of a domain, the weakest unfolds first and the strongest domain last (Rief *et al*. 1997; Bullard *et al*. 2006).

Although larval *Manduca* muscles are striated, they also have features that are reminiscent of tonic or smooth muscles. The Z-lines are indistinct and poorly aligned and each thick filament is surrounded by 10–12 thin filaments (Rheuben & Kammer 1980). It is not known if invertebrate titin-like molecules such as kettin are expressed by *Manduca* muscles; however, they are certainly present in Lepidoptera (Suzuki *et al*. 1999) and other insects groups and appear important for muscle passive properties (Kulke *et al*. 2001). Our current results do not distinguish between the different molecular mechanisms of strain-rate dependency, but it is interesting to consider some of the possibilities. Test data have shown that the unfolding of coiled domains affects the macroscopic stress–strain response (Qi *et al*. 2006; Bertoldi & Boyce in press). Our test data are consistent with such unfolding of domains during the initial loading, and with partial refolding when the muscle returns to the reference configuration upon unloading. Furthermore, the mechanical contributions of these different domains are expected to be highly strain rate dependent, perhaps contributing to the changes in stiffness observed in *Manduca* muscle. Although we have not systematically studied recovery time, it appears that muscles return to an unconditioned state within several minutes of repeated strain cycles. It will be interesting to compare such recovery with that of titin, in which the refolding of a single domain occurs within seconds (Carrion-Vazquez *et al*. 1999).

We do not attempt to formulate a theory based on the worm-like chain (WLC) model to capture the distinct fall-offs in load associated with the unfolding of individual domains. The WLC model, originally introduced by Kratky & Porod (1949), is described in detail in the appendix G of Flory (1988). Since the molecular architecture of the *Manduca* muscle is not sufficiently known, we prefer a phenomenological approach to describe the effect of unfolding and the associated increase in the contour length of the macromolecular chain. This is accomplished by introducing the material parameter *c*_{0} with the dimensions of stress, which serves as a measure of the passive strength of the oriented proteins in the molecular structure of the muscle. Unfolding of individual domains increases the contour length of the macromolecular chain and makes the response less stiff (more compliant). We propose to use the square of the macroscopic stretch *I*_{4} in the fibre direction as a representative measure of domain unfolding, such that(4.2)where is a stress-like material parameter. An appropriate choice of the dimensionless parameters *α*_{1} and *α*_{2} enables one to describe the increase in compliance observed during unfolding.

### 4.2 Hysteretic response

The constitutive model of the *Manduca* muscle must be able to describe the experimentally observed hysteretic response and the associated energy dissipation during loading–unloading cycles. Following the development by Dorfmann *et al*. (2007), we again use the general theory of pseudo-elasticity outlined in §3.6 and select a specific form of the pseudo-energy function (3.24). For the *Manduca* muscle, considering an incompressible and transversely isotropic material, we propose a pseudo-energy function given by(4.3)where the subscript ‘0’ refers to the loading paths on which the variable *η* is inactive and equal to unity. The function *ϕ*, which depends only on *η*, accounts for the energy dissipation during a loading–unloading cycle and for consistency with equation (3.25), it must satisfy the condition *ϕ*(1)=0. Equation (3.28) then gives the important relation(4.4)which for a specific strain energy formulation, *W*_{0} determines *η* implicitly in terms of the invariants *I*_{1} and *I*_{4}. Unloading, which may be initiated from any point on the loading path, is taken as a signal to activate *η*. We emphasize that the magnitude of *η* derived from equation (4.4) depends on the maximum values of *I*_{1} and *I*_{4} attained on the loading path, which we denote and . Since *η*=1 at any point on the loading path from which unloading is initiated, it follows from equation (4.4) that(4.5)wherein the notation *W*_{m} is defined. This is the current maximum value of the energy achieved on the loading path. In accordance with the properties of *W*_{0}, *W*_{m} increases along a loading path. In view of (4.5), the function *ϕ* depends on the point from which unloading begins through the energy expended on the loading path up to that point.

When the material is fully unloaded and in the reference configuration, with *I*_{1}=3 and *I*_{4}=1, the magnitude of the variable *η* is a minimum. This is determined by inserting these values into equation (4.4) to give(4.6)where we assume that no elastic energy is stored in the reference configuration. The residual (non-recoverable) energy *ϕ*(*η*_{min}) is given by (4.3) and has the value(4.7)This may be interpreted as a measure of the energy dissipated in the muscle during the loading–unloading cycle and represents the area between the loading and the unloading curves. It is therefore appropriate to denote *ϕ* as a *dissipation function*. Following Dorfmann & Ogden (2003), we select the derivative of the dissipation function *ϕ* to have the form(4.8)where *r* and *m*/*μ* are dimensionless positive material parameters, *μ* being the shear modulus of the matrix material. The explicit form of *η* is obtained from equations (4.8) and (4.4) and has the form(4.9)The variable *η* assumes the minimum value *η*_{min} when *I*_{1}=3 and *I*_{4}=1, which corresponds to the reference configuration of the muscle. It is given by(4.10)Finally, integration of equation (4.8) gives the energy dissipation explicitly in terms of the variable *η* in the form(4.11)

### 4.3 Rate-dependent response

Experimental data of the rate-dependent behaviour of the *Manduca* muscle have been shown in §2. To model the rate-dependent loading and unloading responses, we propose a modified pseudo-elastic strain energy formulation that, in the absence of strain-rate effects, reduces to the formulation (4.3).

We propose that the pseudo-elastic energy function for an incompressible and transversely isotropic material has the form . Following Sweeney & Ward (1995) and Selvadurai & Yu (2006), the variable is an effective strain rate defined by(4.12)where the exponent *β* is a material parameter. The modified stretches satisfy the conditions(4.13)This implies that the strain-rate effect is activated only when the material is subjected to extension *λ*_{i}≥1. For uniaxial extension in the one-direction, for example, we have *λ*_{2}=*λ*_{3}<1 and the effective strain becomes , i.e. it reduces to uniaxial strain.

Let the rate-dependent material response be associated with domain unfolding and refolding as described in §4.1. This implies that only the reinforcing term *W*_{fib} in (4.1) depends on the loading rate (see also Bertoldi & Boyce in press). We propose to replace the parameter *α*_{1} in equation (4.2) and the parameter *r* in equation (4.9) by(4.14)to affect, respectively, the rate of unfolding and refolding of domains during loading and subsequent unloading. The superscript on and *r*′ indicates dependence on the loading rate.

### 4.4 A specific model

The theory of pseudo-elasticity is very general and allows for substantial flexibility in the selection of *W*_{0} in equation (4.1). In this section we describe a specific form of the constitutive formulation to model the cyclic response of a preconditioned *Manduca* muscle. In particular, we propose the formulation(4.15)where the neo-Hookean model is used to determine the isotropic response of the microfibrillar base matrix and *μ*(greater than 0) is the associated shear modulus. The embedded fibrous network, which in the case of the *Manduca* muscle includes actin, myosin and other proteins, provides a strong stiffening effect during axial extension. This motivates the use of the exponential function in equation (4.15) with *c*_{0} defined by equation (4.2) and *c*_{1} being a dimensionless positive material parameter. The exponential part is similar to the formulation used for the mechanical response of collagen fibres in arterial walls as described in detail by Holzapfel *et al*. (2000). Note that the strain energy formulation of the fibrous network given by (4.15) differs from the one used previously by Dorfmann *et al*. (2007).

#### 4.4.1 Simple tension

It is convenient to specialize the above equations for simple extension in the fibre direction. We take the principal stretches *λ*_{2}=*λ*_{3} and use the notation(4.16)where *λ* denotes the stretch in the fibre direction and the incompressibility condition (3.5) is automatically satisfied. From equations (3.7)_{1} and (3.17)_{1}, we have(4.17)and the pseudo-strain energy (4.3) then depends on *λ* only. This is given by(4.18)and defines the notation . The superscript on *η*′ again indicates dependence on the loading rate. The variable *η*, corresponding to (4.9), for a rate-dependent response has the specific form(4.19)with the variable *r*′ defined in equation (4.14)_{2}. We write the strain energy (4.15) as a function of *λ* and, using equation (4.18) with *η*′=1, obtain the expression(4.20)which describes the loading path on which the variable *η*′ is inactive. The rate-dependent stress-like material parameter , corresponding to equation (4.2), has the form(4.21)with the parameter given by (4.14)_{1}. The effective strain rate for uniaxial extension has the reduced form(4.22)where is the engineering strain rate used in §2.2. The associated Cauchy stress is then given by(4.23)where the subscript ‘0’ again indicates that this expression describes the loading path in simple tension and(4.24)Unloading may take place from any point on the loading path. The start of unloading is taken as the signal for *η*′ to be active and to change the form of the energy function as shown by equation (4.18). The Cauchy stress during unloading is then given by(4.25)which shows that stress softening during unloading requires that *η*′≤1, with equality only at the point where unloading is initiated.

## 5. Results

To examine the quality of the model, we perform a numerical evaluation of the stress–strain relations given by equations (4.23) and (4.25), respectively, for the loading and unloading responses. In addition to the stress–strain relations, we evaluate equations (4.19) and (4.21) together with (4.14) to account for the rate dependence of the material and the energy dissipation during cyclic loading.

To find the adjustable parameters of the governing equations, we use experimental data from cyclic loading tests of an unstimulated muscle with constant maximum stretch of *λ*_{max}=1.18 and imposed constant values of engineering strain rates 0.0144, 0.072, 0.36 and 1.8 s^{−1} (see §2.2).

For each value of , the muscle is subjected to a total of five loading and unloading cycles of preconditioning up to a pre-selected extension of *λ*_{max}=1.18. The preconditioning is performed to monitor the progression of stress softening and to determine the ultimate stress–deformation response for stretches up to *λ*_{max}. After the initial four loading–unloading cycles, the stress–deformation response is essentially repeatable and no additional stress softening occurs. After the preconditioning cycles are completed, the subsequent loading–unloading cycle is used to evaluate and compare the muscle response for =0.0144, 0.072, 0.36 and 1.8 s^{−1} up to *λ*=*λ*_{max} (figure 2).

The explicit formulations (4.14) of the parameters and *r*′ as a function of time, used respectively in equations (4.19) and (4.21), have not yet been specified. We found that the experimental data in figure 2 are best approximated by the phenomenological relations(5.1)where *α*, *α*_{d} and *r*, *r*_{d} are the adjustable material parameters, and =0.0144 s^{−1} and =1.8 s^{−1} are, respectively, the minimum and maximum values of the loading rate used during the experimental characterization of the muscle.

To determine the values of the material parameters, we start by first considering the experimental data of the loading path corresponding to the largest strain rate =1.8 s^{−1}. Equation (5.1)_{1} shows that the value of *α*_{d} does not contribute to the numerical response for . We fix the values of *α* and *α*_{2} and define some intervals [0,*μ*_{max}], and [0,*c*_{1max}], where the best-fit parameters of *μ*, and *c*_{1} are assumed to be located. The best-fit values of these parameters are subsequently determined iteratively by applying the least-squares method to increasing values of *μ*, and *c*_{1} within these intervals. Next, the initial values of *α* and *α*_{2} are replaced and the procedure repeated until the numerical results of the stress–deformation response approximate the experimental data (figure 3). Then, the loading responses for 0.0144, 0.072 and 0.36 s^{−1} are used to determine the value of *α*_{d}.

To determine the value of *r*′ and *m*, we start with the unloading response corresponding to for which *r*′=*r* (see equation (5.1)_{2}). The least-squares method is used to determine the best-fit values of *r* and *m*. Then, the value of the remaining parameter *r*_{d} is found by matching the unloading responses corresponding to 0.072, 0.36 and 1.8 s^{−1}. The values of the adjustable parameters used to fit the experimental data are listed in table 1. Notably, the best-fit value of the shear modulus of the oriented protein structure is equal to 49.1 kPa. The shear modulus *μ* of the matrix material, including intracellular microtubules and other intramuscular connective tissues, is given by 37.7 kPa (Suga *et al*. 2003; Prado *et al*. 2005).

Circles in figure 3 indicate the experimental loading and unloading data corresponding to an engineering strain rate of 1.8 s^{−1}. The results of the numerical calculations for the same strain rate are given by continuous curves. It can be seen that the fit of the model to the data is good. Figure 4 shows the numerical results for engineering strain rates of 0.0144, 0.072, 0.36 and 1.8 s^{−1}. Comparing the data in figure 2 with the numerical results in figure 4 shows that the fit of the model to the experimental data of the *Manduca* muscle is good for all considered strain rates and that the main characteristics of the behaviour are reproduced.

## 6. Concluding remarks

Experimental data on the rate-dependent response of an unstimulated muscle of the tobacco hornworm caterpillar, *Manduca sexta*, have been reported for engineering strain rates of =0.0144, 0.072, 0.36 and 1.8 s^{−1}. For each value of , the muscle was subjected to a total of five loading–unloading cycles of preconditioning up to a pre-selected maximum extension of *λ*_{max}=1.18. The last loading–unloading cycle was used to evaluate and compare the mechanical behaviour of the muscle. The data show dependence on the loading rate, large nonlinear elastic deformations, a hysteretic response during loading–unloading and stress softening during the first few cycles of repeated loading.

These strain rates were chosen to represent the range of muscle stretching speeds experienced during normal *Manduca* behaviour. When the caterpillar is exploring the substrate VIL strain rates are less than 0.1 s^{−1}, but during the rapid defensive strike reflex the rate can be as high as 2 s^{−1} for a brief period (Walters *et al*. 2001; and Michael Simon 2007, unpublished observations). Several features of the passive muscle properties described here are expected to be important for these behaviours. At all strain rates, the muscle was found to stress-soften with repeated cycles of loading and unloading. This softening recovered spontaneously when the muscle was left at its resting length for several minutes. This property could be important in hydrostatic animals such as *Manduca*. During inactivity, the increased passive stiffness of muscle will help to maintain the internal pressure and provide static stability. Once the muscle is stretched and released (presumably as it initiates movements), the decrease in passive stiffness will decrease the force required to re-stretch the muscle after shortening. Another feature was the increasing stiffness that occurred throughout periods of rapid loading. This is expected to be particularly important during high-velocity movements such as the strike reflex. In this behaviour, a strong mechanical stimulus (such as a pinch) applied to the posterior abdomen causes the caterpillar to contract major longitudinal muscles (such as VIL) on the stimulated side to bring the head quickly around to the contact site. Muscles on the opposite side are stretched rapidly, as much as eight times faster than during crawling. In the tests reported here, the initial (0–5% strain) stiffness was similar for all strain rates, but as stretching continued the stiffness increased much more rapidly for high-velocity and large magnitude strains. During a strike behaviour, the initial low stiffness of muscles on the unstimulated side would help the muscles on the stimulated side to accelerate lateral movement of the head. However, the relatively long distance from the abdominal muscles to the head produces a large moment and a potentially damaging fast stretch of the muscles and body wall tissues. The increased stiffness at high strain rates will passively resist such stretching. Together, these two properties illustrate how nonlinear passive properties of tissues such as muscles can automatically produce adaptive and behaviourally relevant responses without direct commands from the central nervous system.

One important aspect of muscle functions is that different molecular mechanisms are responsible for the time-dependent processes. In this paper we evaluate the response of the muscle for different, constant loading rates and associate the different macroscopic stress–deformation responses with changes in the molecular architecture. In particular, we consider unfolding of coiled domains during initial loading and partial refolding when the muscle returns to the reference configuration upon unloading. A systematic evaluation on the effect of recovery time during which a preconditioned muscle returns to an unconditioned state is currently being performed and is therefore not part of this paper. The associated experimental data showing the time-independent equilibrium configuration of a passive muscle and the time-dependent viscoelastic recovery at constant elongation from the preconditioned to the unconditioned state will be given elsewhere.

The experimental data reported in this paper have then been used to formulate a constitutive model for the mechanical response of the *Manduca* muscle at finite strains. Stress–strain relations for loading and unloading were developed using the general theory of a hyperelastic, transversely isotropic material. The theory has been modified to account for the hysteretic response of a preconditioned muscle during loading–unloading. Phenomenological relations were developed and included in the model to account for the molecular mechanisms responsible for the rate-dependent material response. Good agreement is demonstrated between the experimental data and the results of the numerical simulation.

## Acknowledgments

This research was supported by the National Science Foundation grants IBN 0342330 and IBN 0117135 to B.A.T. and a W.M. Keck Foundation Science and Engineering Program grant Biomimetic Technologies for Soft-bodied Robots.

## Footnotes

- Received April 13, 2007.
- Accepted June 1, 2007.

- © 2007 The Royal Society