Three-Dimensional Muscle Architecture and Comprehensive Dynamic Properties of Rabbit Gastrocnemius, Plantaris and Soleus: Input for Simulation Studies

The vastly increasing number of neuro-muscular simulation studies (with increasing numbers of muscles used per simulation) is in sharp contrast to a narrow database of necessary muscle parameters. Simulation results depend heavily on rough parameter estimates often obtained by scaling of one muscle parameter set. However, in vivo muscles differ in their individual properties and architecture. Here we provide a comprehensive dataset of dynamic (n = 6 per muscle) and geometric (three-dimensional architecture, n = 3 per muscle) muscle properties of the rabbit calf muscles gastrocnemius, plantaris, and soleus. For completeness we provide the dynamic muscle properties for further important shank muscles (flexor digitorum longus, extensor digitorum longus, and tibialis anterior; n = 1 per muscle). Maximum shortening velocity (normalized to optimal fiber length) of the gastrocnemius is about twice that of soleus, while plantaris showed an intermediate value. The force-velocity relation is similar for gastrocnemius and plantaris but is much more bent for the soleus. Although the muscles vary greatly in their three-dimensional architecture their mean pennation angle and normalized force-length relationships are almost similar. Forces of the muscles were enhanced in the isometric phase following stretching and were depressed following shortening compared to the corresponding isometric forces. While the enhancement was independent of the ramp velocity, the depression was inversely related to the ramp velocity. The lowest effect strength for soleus supports the idea that these effects adapt to muscle function. The careful acquisition of typical dynamical parameters (e.g. force-length and force-velocity relations, force elongation relations of passive components), enhancement and depression effects, and 3D muscle architecture of calf muscles provides valuable comprehensive datasets for e.g. simulations with neuro-muscular models, development of more realistic muscle models, or simulation of muscle packages.


Introduction
About 600 skeletal muscles enable complex movements like locomotion or laughing in humans. Musculoskeletal models with an increasing number of muscles have been used to investigate human and animal movements [1][2][3][4]. Although many phenomenological Hill-type [5][6][7][8] and biophysical [9,10] muscle models have been developed to represent the dynamics of isolated muscles, the phenomenological modeling approach dominates in musculoskeletal modeling for simplicity and low computational cost. Because individual muscle parameters are lacking, a common approach in such studies is to scale one set of muscle parameters, e.g. by maximum isometric force and fiber length, to fit all muscles [2][3][4]. However, muscles differ in their individual muscle properties not only by scale, and these differences impact simulation results [11].
The active contractile properties of muscle fibers can be characterized by the hyperbolic force-velocity relation [5] and the force-length relation [12]. It has generally been found that slow twitch fibers exhibit a more bent force-velocity relation and lower maximum shortening velocity than fast ones [13]. Moreover, it is known that muscles exhibit a complex architecture consisting of fascicles with different lengths and pennation angles [14][15][16]. The specific architecture influences the shape of the muscle's force-length relation [17]. Thus, simple scaling of these relations introduces errors of unknown magnitude. A broader database of lumped muscle model parameters for specific muscles would improve the situation; in parallel, development of more realistic three-dimensional (3D) muscle models requires individual muscle architecture as input [18][19][20].
Passive muscle properties vary considerably between muscles too [21]. For instance, passive forces of frog gastrocnemius arise already on the ascending limb of the force-length relation and reach about 30% maximum isometric force (F im ) at optimal muscle length [22]. In contrast, passive forces of the frog semitendinosus muscle arise on the descending limb of the force-length relation [23]. Depending on muscle function tendons exhibit different mechanical properties [24], e.g. the tendon of the human wrist mover M. extensor carpi radialis longus used for positioning tasks is very stiff (1.8% strain at F im , [25]) compared with the more compliant human gastrocnemius tendon exploited for elastic recoil (5% strain at F im , [26]).
In Hill-type muscle models, the active contractile properties of the fibers are represented by a contractile component (CC). Its length is usually taken to be the mean optimal fiber length, and sometimes the mean pennation angle is considered [6]. Passive tissues in parallel to the fibers like connective tissue and titin (though titin may be considered to be a semi-active element, [27]) can be represented by a parallel elastic component (PEC). Tendon and aponeurosis thought to act in series with the fibers [22] are represented by a serial elastic component (SEC). In the structurally more convincing arrangement of these three components [28] the SEC is in series to both the CC and the PEC (Fig 1). Typical constitutive functions describing the components are depicted in Fig 1. Since more than 60 years it is known that muscle force further depends on contraction history [30]. For example, force is enhanced in the isometric phase following active stretching (force enhancement, FE) and depressed following active shortening (force depression, FD) compared with the corresponding isometric muscle force. Force enhancement effects can be much larger (2 F im at lengths with no filament overlap, [31]) than force depression effects (0.05-0.2 F im , [30,32]). Currently, the causes of these history effects remain a matter of scientific debate [8,[33][34][35]. Discussed mechanisms are e.g. the contribution of half-sarcomere chain dynamics [34,36,37] or non-cross-bridge contributions to muscle tension [27,31,38]. Muscle specific differences in contraction history are rarely examined, particularly as experimental protocols and conditions differ, which hampers comparison.
The calf musculature is frequently used as a research object in muscle experiments and simulations (gastrocnemius (GAS), soleus (SOL) and/or plantaris (PLA); e.g. [21,[39][40][41]). These distal muscles enable comparably easy surgical access, and their long distal tendons simplify fixation to force measurement equipment in animal models. Some of these muscles were observed in recent in vivo studies with respect to their 3D behavior [42,43]. Understanding e. g. 3D muscle deformations and related force effects in muscle packages by neuromuscular simulations requires consistent specific muscle properties and muscle architectures. Gained insights into fundamental 3D muscle functions may then be generalized which makes the collection of comprehensive data sets important from a clinical perspective. To our knowledge, such comprehensive consistent muscle properties are not available for the calf musculature.
The aim of this study is to provide comprehensive data sets of specific muscle properties (force-length relation, force-velocity relation, force-strain relations of SEC and PEC, activation time constant) and the 3D architecture of the superficial rabbit calf muscles for future research. To achieve this, we determine muscle properties of GAS, PLA, and SOL in in situ experiments (n = 6 per muscle) and measure the 3D architecture by manual digitization ( [16]; n = 3 per muscle). Because there is no generally accepted model of history effects, we provide standardized data for force enhancement and force depression (n = 3 per muscle) which can be used to adapt parameters of custom models describing these effects. Striving for completeness with respect to the shank musculature (Fig 2), we provide the muscle properties for further shank muscles (flexor digitorum longus (FDL), extensor digitorum longus (EDL) and tibialis anterior (TA), n = 1).

Experimental setup
Experiments on female New Zealand white rabbits (Oryctolagus cuniculus, n = 21, age: about 16 weeks) with an average weight of 3.37 ± 0.51 kg (mean ± SD, Table 1) were carried out in strict accordance with the recommendations of the German animal welfare law (Tierschutzgesetz, BGBl. I 1972, 1277, section 8). The protocol of this study was approved by the competent authority for animal welfare in Thuringia, Germany (Landesamt für Verbraucherschutz (Abteilung Gesundheitlicher und technischer Verbraucherschutz); Permit Number: 02-022/11 and 02-027/14). All experiments were performed under anesthesia with natrium pentobarbital (Nembutal, 80 mg/kg body weight) and Bupivacain (1 ml, 0.5%, epidural), and all efforts were made to minimize suffering. Experimental setup, anesthesia and preparation of rabbit SOL have been described earlier [43]. Procedures were similar for all muscles analyzed in this study ( Table 1). In short, the specific muscle was freed from its surrounding tissues and the rabbit was fixed by clamping hip, knee and ankle with three pairs of bone pins. The distal tendon was attached horizontally to a muscle lever system (Aurora scientific 310B-LR). The sciatic nerve was stimulated (Aurora Scientific 701C) with 100 μs square wave pulses at 100-140 Hz (supramaximal tetanic muscle stimulation). Body temperature was maintained at 39°C using a heating pad. The muscle was sprinkled with heated (39°C) physiological saline solution during the entire experiment.

Experiments for determination of muscle properties
First, the muscle-tendon complex length (L MTC_0 ) was measured in situ with a micrometer at an ankle and knee joint angle of 90° (Table 1). To determine the specific muscle properties (force-length relation, force-velocity relation, force-strain relation of the SEC and PEC, and activation time constant), isometric, isotonic and isokinetic contractions were performed. We identified the force-velocity relation by a series of about 10 isotonic contractions against forces in the range of 0.1 F im to 0.9 F im . Similar to [44], the force-strain relation of the SEC was obtained from the length tension data of a quick (exceeding maximum contraction velocity v CCmax ) isokinetic contraction accounting for fiber shortening. The active force-length relation of the CC and the force-strain relation of the PEC were determined from a series of 15-20 isometric contractions (with length increments of 1-2 mm) considering the interaction of the passive elastic structures [28]. The activation parameter τ, describing calcium concentration and thus muscle activation, was determined from an isometric contraction at optimal muscle length using a Hill-type model [29]. For details of these parameter determinations, see [29,45]. Muscle and animal mass as well as the muscle-tendon complex length L MTC_0 measured at ankle and knee joint angles of 90°(cf. Fig 1). n: number of muscles. *Architecture of GAS, PLA, and SOL was determined from the left legs of three rabbits. doi:10.1371/journal.pone.0130985.t001 History effects were identified for GAS, PLA, and SOL (n = 3 per muscle) using isokinetic ramps [32] with three different shortening / lengthening velocities (0.35, 0.7, and 1.4 l fm /s, where l fm is the mean fascicle length determined by manual digitization; see section 2.4), and a ramp length of 0.3 l fm . The isokinetic ramp started after an isometric pre-contraction (200 to 400 ms) which was sufficient to reach maximum isometric force. After the end of the ramp, stimulation continued (SOL: 1300 ms; other muscles: 500 ms) to allow sufficient force recovery during the subsequent isometric phase. Longer continued stimulation was used for SOL because this slow twitch fibered muscle [46] is more fatigue resistant than the other muscles. Shortening and lengthening ramps for determination of force depression and force enhancement started at optimum muscle length plus 0.15 l fm and minus 0.15 l fm , respectively. Force depression and force enhancement were identified as the difference in force between the force at the end of the ramp experiment and the force at the end of an isometric reference contraction of same duration at the same target length. For FDL, EDL, and TA no muscle architecture and thus no mean fascicle length (l fm ) was determined. Thus, for these muscles isokinetic ramps had velocities of 5, 10 and 20 mm/s starting at optimum muscle length plus 2 mm (force depression) and minus 2 mm (force enhancement) with a ramp length of 4 mm.

Hill-type muscle model
Parameters were determined for a typical Hill-type muscle model [29] consisting of a PEC in parallel to the CC, and a SEC connected to both of them. Using this model, the active force generated by the contractile component F cc is the difference between the single forces in the SEC and the PEC and can be described by a typical product approach [6]: In Eq (1), A is the muscle activation and f l as well as f v are factors that describe the forcelength and force-velocity relations normalized to F im , respectively.
The CC force-length relationship was described with a piecewise linear equation where f c is the force at which the ascending limb changes slope, l CC is the CC length, l CCopt the optimal CC length and l 1 , l 2 , l 3 , l 4 are specific lengths that are crucial for the sarcomere forcelength relationship (Fig 1). The force-velocity relationship follows the Hill hyperbola [5] for concentric contractions, with v CCmax being the maximal CC shortening velocity, and curv = a/F im (damping increases with decreasing curv, see Fig 1; a describes the force asymptote [5]) is an inverse measure of the relation's curvature.
The latency between the supramaximal stimulation and the muscle activation was modelled as a first order linear differential equation [47] dA dt where the activation parameter τ lumps the time constants of calcium influx from the sarcoplasmic reticulum into the sarcoplasm, and A(t = 0) = 0. The SEC force-elongation relationship F SEC (Δl SEC ) was taken from [6] where Δl SEC1 and F 1 are the elongation and the force at which the force-elongation relation changes from exponential to linear. k was calculated from Δl SEC1 , F 1 , the dimensionless shape parameter k sh and the constraint that stiffness at Δl SEC = Δl SEC1 is the same for each equation.
A PEC force-elongation relation depending on k 1 and k 2 was taken from [48].

Determination of muscle architecture
Muscle architecture of GAS, PLA, and SOL was determined from the left legs of three rabbits (R1 to R3, m = 3.28 ± 0.17 kg) by manual digitization [16,43]. After the rabbit was killed with an overdose of pentobarbital, the leg was amputated above the knee. The skin was removed and the preparation fixed in Bouin solution (an aqueous solution of picric acid, acetic acid and formaldehyde minimizing tissue shrinkage [49,50]) for 48 h [14] at knee and ankle angles of about 79°and 93°, respectively. Subsequently, the bone-muscle preparation was cast in wax to provide additional mechanical stability during the digitizing process. For the digitization of the whole muscle architecture of each muscle, small fascicle bundles were successively dissected with a micro forceps. Their original position was then recorded using a manual 3D digitizer (MicroScribe MLX) with a sampling frequency of 5 Hz and an accuracy of 0.07 mm. This process was repeated until all fascicles of the individual muscle were recorded. Each fascicle was described by 20 points. During the dissection and the digitizing, the palm of the hand holding the digitizer-handpiece was placed on the preparation-table to minimize movement of the digitizer tip (< 0.1 mm). Fascicle length and pennation angle were calculated as described in [16]. In addition to recording the fascicle bundle positions, the insertions and origins of the muscles were recorded for each animal (S1, S5, and S9 Datasets).

Statistics
Prior to analysis, muscle parameters were normalized using pertinent units. Parameters were tested for normal distribution using the Kolmogorov-Smirnov-Test with Lilliefors correction. All data were normally distributed. The Levené test was used to check variance homogeneity.
To test whether muscle properties differ between the muscles (GAS, PLA, SOL) an analysis of variance (ANOVA) was calculated. In case that the ANOVA demonstrated significant main effects, post hoc analyses were performed using the Tukey HSD test if variances were homogenous. Otherwise the Tamhane test was used. The significance level was set at P < 0.05. All analyses were performed using SPSS 22 (IBM Corp, Armonk, NY, USA). The effect sizes f were calculated as where σ m is the standard deviation of the population means, and σ the within-population standard deviation [51]. The effect sizes were classified as low (f = 0.1), medium (f = 0.25) and large (f > 0.40) [51].

Muscle properties
Active and passive muscle properties vary considerably between GAS, PLA, and SOL. We found significant differences for the parameters l CCopt , v CCmax , curv, Δl SEC1 /l SEC0 , k, l SEC0 , and l PEC0 ( Table 2). All experimental force-velocity relationships feature the typical hyperbolic shape (Fig 3, second row) observed by [5]. Maximum shortening velocity of GAS (13.5 ± 1.7 l CCopt /s) is about twice the value of SOL (6.4 ± 1.0 l CCopt /s). Maximum shortening velocity of PLA (10.1 ± 3.3 l CCopt /s) is in between these values. The curv values of the force-velocity relation are similar for GAS and PLA (0.47 ± 0.09 and 0.41 ± 0.16, respectively) but about three times the value for SOL (0.15 ± 0.05). The calf muscles exhibited a characteristic force-length dependency (Fig 3, upper row) which is attributable to the muscle fiber force-length relationship [12]. Maximum isometric forces produced at optimum muscle lengths by GAS, PLA, and SOL are 161.3 ± 18.2 N, 86.4 ± 21.3 N, and 24.1 ± 5.8 N, respectively. Considering a muscle tissue density of 1.056 g/cm 3 [52], as well as a mean muscle mass (Table 1) and mean optimal fiber length (  3). The standard deviations of the determined muscle properties are small, with the exception of the force-strain relation of the parallel elastic component (Fig 3, bottom row) which is about three times the standard deviation observed for the SEC. Forces of GAS, PLA, and SOL were enhanced following stretching and were depressed following shortening compared with the corresponding isometric forces ( Table 2, Fig 4). For all muscles, force depression was inversely related to the ramp velocity (Table 2). This effect was pronounced for the SOL. In contrast, force enhancement was independent of ramp velocity. The magnitude of force enhancement increased from GAS (% 8% F im ), to SOL (% 11% F im ) up to PLA (% 17% F im ). At the slowest ramp velocity (0.35 l fm /s), force depression of GAS and PLA (% 17.5% F im ) was about two fold that of SOL (9.7% F im ).
We report additionally obtained muscle properties of further shank muscles (FDL, EDL, TA; n = 1) for completeness in the Supporting Information (S1 Text).

Architecture
General architectural properties of GAS, PLA, and SOL are listed in Table 3. Spatial coordinates of the fascicles of rabbit R1 are presented in Fig 5. Three dimensional fascicle data including origin and insertion of the GAS, PLA, and SOL of the three animals (R1, R2, R3) are provided in txt-format in the Supporting Information (S1-S12 Datasets). SOL (Fig 5, green fascicles) exhibits simple unipennate muscle architecture while GAS (medialis: light red fascicles; lateralis: yellow fascicles) and PLA (dark red fascicles) show more complex bipennate muscle architectures. For each animal, the mean pennation angles of all three muscles (GAS, PLA, and SOL) were almost similar (Table 3). Differences appear to be about 1-2°, only. In between the different animals, variations were slightly larger. Mean pennation angles in R1 are about 5°l arger than in R3. These variations in pennation angles have only small (< 0.02 F im ) impact on the calculation of muscle force. Mean fascicle lengths (l fm ) are larger for the heavier animals (R2 and R3). For each specific animal, GAS and SOL exhibit about the same mean fascicle lengths but in general PLA is about 30% shorter.

Discussion
Experiments performed within this study provide comprehensive data sets for the rabbit calf muscles GAS, PLA, and SOL consisting of manually digitized 3D muscle architectures and specific muscle properties including the quantification of history effects.
The active force-length relation was described by the theoretical sarcomere force-length relationship [12]. As demonstrated in recent studies [28,59] this relation enabled the accurate prediction of experimental rabbit and cat muscle forces. Winters [59] reproduced the active force-length relations of rabbit TA, EDL, and extensor digitorum II based on myofilament lengths using a scaled sarcomere model. In contrast, we fitted the experimental force-length data using a piecewise linear equation. However, the results are consistent for the most part. Starting from optimum muscle length, the rabbit muscles were able to shorten by 50% in both studies. For lengthening muscle, force production is limited to 1.6 l CCopt in the scaled sarcomere model, but reaches about 2 l CCopt for GAS, PLA, and SOL (Fig 3). Differences might be due to more complex muscle architectures, especially of GAS and PLA (Fig 5, see Sect. 4.2), influencing the width of the force-length relation [60]. Also, to avoid muscle damage induced by high passive forces, muscles in our study were lengthened only up to passive forces of about 0.2 F im (Fig 3, upper row, marked by a white circle). Thus, the slope of the descending limb of the force length relation was determined using limited experimental data and should be considered with caution. The change in slope at the ascending limb of the force-length relation is fixed at 0.7 F im in the scaled sarcomere model, but appears at higher forces (f c = 0.85 F im ) in our measurements. However, the change in slope of experimental TA force-length relations appears at (black) and isometric reference contraction (grey) determined 500ms (GAS, PLA) and 1300ms (SOL) after the end of the ramp, shown exemplarily for the slowest (0.35 l fm /s) ramp.
doi:10.1371/journal.pone.0130985.g004 Mean pennation angle and fascicle length determined from the left legs of three rabbits (R1, R2, R3) by manual digitization. 1838, 1773, and 1523 fascicles have been digitized for R1, R2, and R3, respectively. Their lengths and pennation angles are normally distributed. The muscle architecture is exemplarily shown for R1 in Fig 5. Note that the complete 3D data are provided in the Supporting Information (S1-S12 Datasets). Rabbit Calf Muscle Properties higher forces in the study of [59] (their Fig 3A) which is in agreement with our observations on TA (S1 Table). Series elastic structures exhibit a typical [6] nonlinear force-strain relationship (Fig 3, third  row). The mean maximum strain at F im was 5.5% and 3.6% L SEC0 for PLA and SOL (Fig 3), respectively, which is expected for tendinous tissue [61,62]. The SEC of the GAS was more compliant (7.7% strain at F im ). This may be due to a higher proportion of aponeuroses in the muscle tendon complex which may be more compliant than tendons [63].
The standard deviations of the determined muscle properties are small, with the exception of the force-strain relation of the parallel elastic component (Fig 3, lower row). This has also been reported for other muscles [43,64,65] and may be related to variations in connective tissues (fascia, epimysium, perimysium, endomysium) or titin-isoforms (e.g. [66]).
The behavior of the rabbit muscles observed in this study is mostly consistent with history effects observed in other muscles. As found in our study, force enhancement is independent of stretch velocity [67], and force depression decreases with increasing ramp velocity [30]. The influence of these history effects on the determination of the muscle properties is discussed in the Supporting Information (S2 Text).
Using the same experimental setup and conditions we found muscle specific differences in the amount and in the ratio of force enhancement and force depression ( Table 3, S1 Table). These muscle specific properties might be explained by synthesis of different titin isoforms in different rabbit muscles as reported by [66]. Interestingly, these differences in titin isoforms are not related to the fiber type. Interaction of titin with the actin myofilament is assumed to be responsible for the history dependence of muscle contractions [8,27]. However, experimental and modeling evidence is necessary to demonstrate the conceivable relation between muscle specific titin isoforms and muscle specific history dependence of muscle contraction.
The mechanism and function of contraction history effects are a matter of scientific debate [8,[33][34][35]. Force enhancement enables the muscle to withstand high forces during eccentric contractions. Rode [27] suggested that force depression is an unwanted by-product of desired force enhancement, and does not occur in stretch-shortening cycles associated with bouncing gaits. The biarticular PLA has comparably short muscle fibers (Table 2) and long tendons, and is appropriate to work as a spring during hopping [58,68]. Indeed, PLA exhibits comparably high force enhancement (Table 2) which enables generating high forces during eccentric contractions. However, the primary function of muscles working as motors during locomotion, e.g. the monoarticular SOL, is to produce positive work [41,69,70]. For these muscles force depression seems to be counterproductive because it reduces positive work. Indeed, SOL exhibited much lower force depression than GAS and PLA (Table 2). These findings support the idea that history effects represent an adaptation to the specific muscle function [27].

Muscle architecture
Studies providing 3D architectural data of muscle packages consisting of several synergistic muscles are rare. Using diffusion tensor imaging 3D muscle architecture of e.g. the human calf [71], human thigh [72], human forearm [73] and mouse hindlimb muscles [74] have been examined. The architecture of the human back [75] and cavy forelimb muscles [14] was determined by manual digitization. However, to the author's best knowledge, there is no consistent 3D data set of the complete rabbit calf muscle architecture.
Measurements comparable to our experiments were performed on rabbit SOL [43]. The authors examined ten SOL muscles using manual digitization and reported a mean fascicle length of 16.6 ± 2.6 mm which is similar to our value obtained from three muscles (17.1 ± 2.6 mm). A slightly lower mean pennation angle (9.9 ± 2.8°vs. 13.0 ± 2.0°) might be due to smaller ankle joint angles reported in their study. In agreement with our results, Hiepe [76] reported mean fascicle lengths of 16.2 ± 9.1 mm for rabbit GAS medialis using DTI. The higher standard deviation in their study might be due to limitations of the diffusion tensor imaging method, e.g. generated fiber tracts may cross muscle borders due to equally oriented adjacent structures, resulting in artificial fascicle traces that are too long [16,74].
There are few studies dealing with architectural measurements of the PLA [77][78][79][80]. None of them utilizes the rabbit as animal model or provides 3D architectural data what makes a direct comparison of fiber lengths and pennation angle impossible. Savelberg [78] mentioned that the non-parallel arrangement of the aponeuroses inside the rat's PLA features different pennation angles and various fiber lengths. This is in agreement with our observations of the complex bipennate rabbit PLA muscle architecture.
Models of muscle packages are required to understand transverse interaction of muscles with each other [81][82][83] and with the skeleton, internal muscle forces, and the influence of muscle architecture on contraction dynamics and muscle deformation. In addition to deviating muscle architectures, observed significant differences in normalized dynamical muscle parameters suggest that it is important to use specific muscle parameters in neuromuscular models aiming at understanding fundamental 3D muscle functions. The revelation of such fundamental 3D muscle functions may be relevant from a clinical perspective when assessing the effects of muscle malfunction e.g. on the stabilization of joints during movement. Moreover, 3D muscle models may contribute to a deeper understanding of widespread diseases as chronic low back pain which are accompanied by several changes in the muscle structure as atrophy [84], steatosis [85] or altered fiber type distribution [86].