Bone strain magnitude is correlated with bone strain rate in tetrapods: implications for models of mechanotransduction

1. Introduction

The mechanical loading of bones is a potent stimulus affecting adaptive bone modelling [1]. Because high strain magnitudes can increase the probability of a detrimental bone fracture [2], it has been hypothesized that bone adaptation via modelling should serve to decrease strain magnitudes in areas exposed to high loads [3–6]. Available data suggest that bone formation in a range of vertebrate taxa and bones is correlated with several stimuli, including strain gradients [7], strain magnitudes [8–14] and strain rates [15–18]. However, many of these studies were conducted using artificial loading regimes, and do not rely on the bone strain profiles experienced by animals during natural, unrestrained terrestrial locomotion and feeding. Because multiple tissue-level strain stimuli are capable of driving bone modelling, there is ambiguity as to which of these stimuli (or their combinations) leads to the initiation of mechanotransduction at the cellular level. In sum, it is not clear how osteocytes detect local differences in tissue-level strain magnitude, making it unclear how local changes in bone stress and strain magnitudes can elicit bone modelling responses [1,19].

Mechanistic modelling studies suggest that strain magnitudes might be high enough to directly excite osteocyte cell bodies within their lacunae [20,21]; cell-level strains greater than 5000 µε are required to excite osteocytes in vitro [22]. However, tissue-level principal strain magnitudes do not regularly surpass 4000 µε either for limb bones during vigorous locomotion [23] or for the jaws during feeding [24]. Furthermore, it has recently been shown that the osteocyte cell processes (which are housed within canaliculi), rather than the osteocyte cell bodies, are probably responsible for initiating the mechanosensory response [25]. Therefore, it is unlikely that coupling of osteocyte deformation to local bone deformation directly induces the intracellular response of osteocytes to tissue-level bone strain.

In contrast to the lack of mechanistic links between bone adaptive responses and bone strain magnitudes, dynamic loading provides a mechanism for strain amplification at the cellular level by inducing fluid flow through the bone's lacunar–canalicular porosity [26]. The cell membranes of the osteocyte processes are anchored to the canalicular wall by multiple tethering proteins that span the pericellular matrix [27]. Theoretical mechanotransduction models suggest that fluid flow in this matrix generates high drag on the tethering elements [28,29], which generates radial (hoop) strains in the osteocyte processes' cell membranes [29]. These strains have further been proposed to relate directly to the tissue-level modelling response [30]. In this way, models predict that interstitial fluid flow through a bone's lacunar–canalicular porosity is able to generate more than five times greater strain at the cell membrane than is required to elicit an intracellular response [29].

Mechanotransduction models also propose that the magnitude of drag experienced by the tethering elements increases in direct proportion to the velocity of fluid flow through the canaliculi [28]. Turner et al. [17] analogized compact bone with a water-soaked sponge and argued that the velocity of a bone's interstitial fluid will be related to the rate at which the bone experiences tissue-level mechanical strain. Turner et al.'s analogy has been supported experimentally in studies that have found a bone's streaming potential, which is caused by interstitial fluid flow, to increase with loading frequency [17,31–33]. Previous studies have also suggested that the velocity of interstitial fluid flow through the canaliculi is related to strain gradient [7], which is ultimately influenced by strain magnitude. However, numerous studies have found that new bone formation is directly proportional to strain rate when strain magnitude is held constant [15–17]. Drawing from these studies, it appears that large strain magnitudes are not capable of initiating bone modelling unless they are coupled to large strain rates [18].

Consideration of the data summarized above reveals a paradox. Although organisms would benefit from a mechanism linking adaptive bone modelling to changes in strain magnitude, the mechanism for detecting strain magnitude at the tissue level is unclear. Moreover, if the dynamic loading model provides a cell-based mechanism for modelling bone form, what is the selective advantage of linking adaptation in bone form to strain rate? Finally, it is unknown how the two potent tissue-level bone strain stimuli capable of inducing bone modelling—strain magnitude and strain rate—are related to each other during natural, unrestrained terrestrial locomotion. Based on our previous findings that bone strain magnitude is significantly correlated with strain rate in the tetrapod feeding system [34,35] and in the derived locomotor system of the river cooter turtle [23], we hypothesized that limb bone strain magnitude is highly correlated with strain rate, providing a link between tissue-level strain magnitudes and cellular-level fluid-flow rates. To test this hypothesis, we analysed bone strain data from a variety of tetrapod species, limb bones and locomotor behaviours to determine whether variation in strain magnitude is highly correlated with variation in strain rate.

Strain magnitude could be increased by increasing strain rate while load duration is held constant, increasing load duration while strain rate is held constant, or some combination of the two [34]. Evidence of modulation of strain magnitude via loading duration comes from correlations between muscle activity duration and increasing functional demands on the musculoskeletal system. For example, the lizard Chamaeleo calyptratus adapts to increased incline, which is associated with an increased power requirement [36], by increasing EMG amplitude and activity duration of the gastrocnemius and tibialis anterior muscles with no associated change in kinematics [37]. Furthermore, EMG burst duration, relative to undulatory cycle time, increases in anguillid eels during locomotion across land, which probably requires increased force production compared with aquatic locomotion [38]. Muscle contraction places loads on associated skeletal elements [39–41] and, in the hindlimb, muscle contractile activity is correlated with peak bone strain magnitude [42]. Therefore, increases in the duration of muscle activity (load duration) may lead to increased bone strain magnitudes.

However, evidence that strain magnitude is likely to be modulated via changes in load rate across a wide range of vertebrates comes from observations that the locomotor and feeding systems of many tetrapods operate highly rhythmically [35,43]. Low variation in locomotor cycle duration implies that variation in forces and strains must be accommodated through variation in load rate rather than in load duration. Moreover, muscle force modulation during locomotion occurs through the orderly recruitment of motor units [44–46], whereby increased muscle force generation is achieved by recruitment of progressively larger and faster motor units. This provides a motor-control mechanism to explain how rate modulation of limb bone loading occurs [34].

This study is the first to test whether the magnitude of tissue-level bone strain is encoded in the rate at which the bones are loaded across a wide range of vertebrate taxa, bones and locomotor behaviours. This would indicate that the way that force is modulated in musculoskeletal systems during organism-level behaviours (rate modulation) encodes the magnitude of tissue-level strains in a manner that cellular-level processes can and do detect. The strain rate encoding of strain magnitude across a phylogenetically broad group of tetrapod species and two different functional systems (feeding and locomotion) would suggest that this is a fundamental mechanism of functional skeletal adaptation in tetrapods.

2. Results

(a) Bivariate analyses

Regardless of strain type (ɛ₁ principal strain, ɛ₂ principal strain and shear strain), the absolute magnitude of strain was always significantly positively correlated with strain rate (slopes ranging from 0.017 to 0.757; see M–R rows of electronic supplementary material, tables S1–S3) and only sometimes correlated with load duration across all species and limb bone elements (figure 1; electronic supplementary material, tables S1–S3). By contrast, ɛ₁, ɛ₂ and shear strain magnitude were only correlated (p < 0.05) with load duration in three of nine, four of nine and four of six species, respectively. On average, strain rate explained 49.29 ± 18.75% (mean ± s.d.) more of the variation in strain magnitude than load duration. When strain magnitude was correlated with both strain rate and load duration, strain rate always explained more of the variation in strain magnitude than load duration (figure 1). Individual effects were evident in the bivariate analyses both qualitatively and quantitatively (figure 2; electronic supplementary material, table S4). However, individual trends mirrored collective trends for the species, showing more (and stronger) correlations between strain magnitude and strain rate, compared with load duration (figure 2).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

(b) Multivariate analyses

Multivariate analyses were conducted to determine which independent variable explained more strain magnitude variation when individual effects and interaction effects between the independent variables (strain rate and load duration) were taken into account (figure 3; electronic supplementary material, tables S1–S3). The summary of fit for each multiple regression always yielded r² > 0.75, with 21 of 23 fits having r² > 0.90. Strain (ε₁, ε₂ and shear) magnitude was always significantly correlated with both strain rate and load duration. Strain rate β coefficients (standardized partial slope) were 1.05–2.06 (ɛ₁), 1.09–2.03 (ɛ₂) and 1.06–1.86 (shear) times higher than load duration β coefficients (figure 3). The presence of a larger β coefficient for strain rate suggests that variation in this factor always explains more strain magnitude variation than does load duration.

• Download figure
• Open in new tab
• Download PowerPoint

Interaction effects between strain rate and load duration were always significant for each species's regression of ɛ₁, ɛ₂ and shear strain. However, the variance inflation factors (VIFs) for strain rate and load duration were usually low (below 3.0), suggesting minimal effects from multicollinearity, and independent effects of strain rate and load duration on strain magnitude. With one exception (ɛ₂ strain data from the radius of the goat Capra hircus), the weak relationship between strain rate and load duration was negative.

(c) Multivariate analyses accounting for speed

In order to test whether the relationship between strain magnitude and strain rate is maintained across changes in locomotor speed, multivariate analyses were conducted on strain data from the tibiotarsus of the emu (Dromaius novaehollandiae) at duty factors of 0.65, 0.55 and 0.40 (figure 4; electronic supplementary material, table S5–S7). Regardless of duty factor, ɛ₁ and ɛ₂ strain magnitudes were always significantly correlated with both strain rate and load duration. However, as in our previous analyses, the β coefficient for strain rate was always larger than that for load duration (figure 4). The strain rate β coefficients from regressions of the ɛ₁ data were 1.50, 1.78 and 1.28 times larger than the load duration β coefficients at duty factors of 0.65, 0.55 and 0.40, respectively. For the ɛ₂ data regressions, the strain rate coefficients were 1.72, 1.32 and 2.03 times larger than those for load duration at these same duty factors.

• Download figure
• Open in new tab
• Download PowerPoint

Emu tibiotarsus data were further pooled into two sets (one for ɛ₁ and one for ɛ₂) that each included all three duty factors. In each case, the strain rate β coefficient was greater than the load duration β coefficient (figure 4) by a factor of 1.23 for ɛ₁ and 2.06 for ɛ₂. In order to fully account for the effects of duty factor, a second multivariate analysis was conducted that included cycle duration (as a proxy for speed) as an independent variable. The results remained consistent with the previous analysis. In the latter analysis, strain rate and load duration β coefficients changed from 0.83 to 0.81 and from 0.68 to 0.66 for the ɛ₁ dataset, and from 1.01 to 0.99 and from 0.49 to 0.47 for the ɛ₂ dataset (figure 4). Thus, inclusion of cycle duration in the regression model only increased the ratio of the strain rate β coefficient to load duration β coefficient by a factor of 0.007 for the ɛ₁ dataset and 0.031 for the ɛ₂ dataset.

3. Discussion

In the tetrapod feeding system, bone strain magnitude is significantly correlated with the rate at which that strain develops [34,35]. Such a correlation also has been identified in a single study of a highly derived locomotor system (turtles) operating at low speeds [23]. This study evaluated whether this relationship holds more generally across the locomotor systems of a broader range of taxa and locomotor styles, and whether strain rate or load duration explained more of the variation in strain magnitude across this sample. Our analysis of locomotor strain data across diverse tetrapod species, limb bones and locomotor styles shows that limb bone strain magnitude is always significantly correlated with strain rate, but not always with load duration, and that strain rate (rather than load duration) explains more variation in strain magnitude. Thus, high correlations between strain magnitude and strain rate are a general feature of tetrapod bone loading in locomotor and feeding systems, whether in cyclic loading events (e.g. mammal chewing; mammal, bird, turtle and alligator walking) or in discrete loading events (e.g. frog jumping).

Our finding that tissue-level strain magnitudes are rate modulated has particular salience in the context of the fluid-flow mechanotransduction model for bone adaptation. The fluid-flow model proposes that the magnitude of drag experienced by osteocyte tethering elements is directly proportional to the velocity of interstitial fluid flow through the canaliculi [28], which is in turn related to tissue-level strain rate [17]. Consequently, the magnitude of the intracellular response of osteocytes to drag imposed on the cell process tethering elements and axial strains of integrins [29] is directly influenced by the rate at which tissue-level strains develop. However, until now, variation in tissue-level strain rates has not been linked to variation in tissue-level strain magnitudes, so it has not been clear how the modelling response elicited by interstitial fluid flow could be linked to the strain magnitudes that determine bones' risk of failure. Strong correlations between strain magnitude and strain rate in both feeding [34,35] and a diversity of locomotor systems (figure 5) provide a basis for hypothesizing such a link. Because strain magnitude also partly determines the velocity of interstitial fluid flow (via strain gradient), we further hypothesize that a correlated increase in strain magnitude and strain rate would help to prevent ambiguity in the encoding of tissue-level strain magnitude at the cellular level through fluid-flow velocity. In summary, our data are consistent with a hypothesis that links tissue-level, in vivo bone loading regimes during natural, unrestrained behaviours with theoretical, cellular-level models of mechanotransduction [28,29], in which tissue-level strain magnitude is encoded at the cellular level by tissue-level bone strain rates.

• Download figure
• Open in new tab
• Download PowerPoint

Bone is a viscoelastic material, with material properties such as Young's modulus increasing with increasing strain rate [47]. Nonetheless, there is little indication that such properties would impact our conclusions. In this study, we found that species-average strain rate ranged from 439 µε s⁻¹ in the tegu (Tupinambis merianae) to 18 966 µε s⁻¹ in the emu (D. novaehollandiae). However, available data suggest that Young's modulus only increases slightly over the three orders of magnitude that bracket our range of in vivo strain rates [48]. For example, the Young's modulus of the bovine femur is 15.2, 17.2 and 17.9 GPa at strain rates of 1000, 10 000 and 100 000 µε s⁻¹, respectively [48]. Thus, at a load of 50 MPa, bovine femoral bone strain decreases from 0.00329 to 0.00279 at strain rates of 1000 and 100 000 µε s⁻¹, respectively. These strain-rate-induced differences in strain magnitude would introduce an encoding error at the cellular level of approximately 15% (similar to the sensory noise level observed in the visual system [49]) when strain rate spans three orders of magnitude. This bracket of experimental strain rates encompasses a range of values approximately 80 000 µε s⁻¹ greater than the largest species-average rate calculated in this study. The extremely high loads and strain rates that occur during the infrequent performance of extreme behaviours might be associated with an encoding error in strain magnitude at the cellular level owing to the viscoelastic properties of bone, but such error would cause strain magnitude to be overestimated, rather than underestimated.

It is possible that rate modulation of bone strain magnitudes could encode information about loading frequency in situations where strain magnitude is relatively invariant because when strain magnitude is held constant, increases in loading frequency are associated with increases in strain rate [50]. Encoding of loading frequency might be advantageous because fatigue loading is known to produce micro-cracks and weaken bone [2,51,52]. However, we do not find this argument compelling because, although the initial loading cycles of a series strongly influence bone formation [53], bone cells rapidly accommodate to ‘routine loading signals’ [50,53], suggesting that bone cells do not record a ‘memory’ of loading frequency. During normal rhythmic locomotion, increases in loading frequency within gaits are also associated with increases in both strain magnitude [54,55] and (as shown here) strain rate. Thus, although resistance to fatigue damage is an important aspect of bone function, the evidence does not support the suggestion that strain rate encoding of load frequency is an important mechanism of fatigue resistance.

During rhythmic locomotion, stance phase and bone loading duration decrease as locomotor speed increases, necessarily resulting in increased ground reaction forces and bone strain magnitudes [55]. Locomotion dynamics therefore almost require that limb bone strain magnitudes are correlated with strain rates and not loading durations. However, it is important to point out that there is no necessary relationship between strain magnitude and loading rate in those musculoskeletal systems that do not support body mass and resist ground reaction forces during cyclic activity. For example, the mammalian feeding system displays high correlations between strain magnitude and strain rate during rhythmic mastication, despite the fact that body-mass dynamics are not an important determinant of feeding-system dynamics [34].

Instead, common mechanisms of bone adaptation can be linked to the orderly recruitment of motor units, which is an important principle of motor control in both locomotor and feeding systems [34,35]. In vertebrate locomotor and feeding systems, as more force is required, small motor neurons, which innervate small motor units consisting of slow twitch fibres, are recruited first, followed by motor units with progressively larger motor neurons and faster fibre types [56–68]. Force modulation through orderly recruitment of motor units has significant advantages for motor control, but it also makes rate modulation of locomotor and feeding forces, and their associated strains, almost inevitable. The demonstration of orderly recruitment of motor units in bony fishes [69] suggests that rate modulation of force may be a widespread and ancient feature of vertebrate musculoskeletal systems. Thus, the evolution of bone adaptation mechanisms that took advantage of this organization may not be surprising.

4. Material and methods

All bone strain data used in this study were recorded as part of unrelated prior studies. Bone strain data were analysed from a variety of tetrapod species, different limb bones and different locomotor behaviours: turtle (Pseudemys concinna) femur [23]; opossum (Didelphis virginiana) femur [70]; chicken (Gallus gallus) femur [71]; goat (C. hircus) radius [72]; emu (Dromaius novaehollandiae) tibiotarsus (TBT) and femur [73]; tegu (T. merianae) femur [74]; alligator (Alligator mississipiensis) humerus [75]; frog (Lithobates catesbeiana) femur [75]. The data were appropriate for the current study because they represented uninterrupted sequences of locomotion consisting of at least five cycles in all species except the frog (L. catesbeiana), which exhibits discontinuous locomotion in which each jump is a discrete loading event. A subset of the data (the TBT of D. novaehollandiae) was analysed across different speeds.

All data analysed in this study were collected from rosette strain gauges attached to the surface of each animal's bone, in accordance with standard methods [76]. Data were sometimes sampled at different rates across the original source studies, due to interspecific differences in typical locomotor speeds. Slight variations were also present in the precise locations of each gauge relative to the neutral axis among individuals within the same study. However, there were no major differences in strain gauge application or data collection methods among the original studies that were the sources of these data, facilitating our comparisons in this analysis. A summary of the strain data from each study, including the sampling rate, number of individuals and step cycles analysed per species, can be found in electronic supplementary material, table S8.

Strain data were analysed using a custom MATLAB (MathWorks, Natick, MA, USA) routine by J.I.-D., who extracted the following variables from each locomotor (step) cycle: peak strain magnitude, duration of strain development (load duration), the rate of strain development (strain rate) and in a subset of the data, step cycle (stance + swing) duration. Calculations and equations for bone strain variables followed previously published methods for mandibular strain [34]. It is rare that bone strain returns to a magnitude of zero during the swing portion of a step cycle because of inertial and muscular forces during limb movement. Therefore, variables were evaluated relative to 25% of peak strain in order to focus analyses on the major loading event that occurs during the step cycle (stance). Load duration was calculated as the time between 25% of peak strain (the time at which 25% of peak strain magnitude was reached during loading) and peak strain. Strain rate was calculated as the average rate of strain development between 25% of peak strain and peak strain magnitude. Step cycle duration was estimated as the average time from the preceding strain peak to the current peak, and from the current peak to the following peak.

Ethics

All experimental procedures were carried out under various Institutional Animal Care and Use Committee guidelines. See original papers for more information; citations can be found in the electronic supplementary material, table S8.

Data accessibility

All original bone strain data are available on request from each contributing author. Publication details for each set of bone strain data can be found in the electronic supplementary material, table S8. The data extracted from the previously published bone strain traces (load duration, strain rate and, in some instances, step cycle) and used in our analysis are available in the electronic supplementary material, table S9.

Authors' contributions

B.R.A., J.I.-D. and C.F.R. designed research; R.W.B, M.T.B., M.T.C., N.R.E. and R.P.M. contributed data; B.R.A., J.I-D. and C.F.R. analysed data; B.R.A., R.W.B., R.P.M., N.R.E. and C.F.R. wrote the paper.

Competing interests

We declare we have no competing interests.

Funding

This material is based on work supported by the National Science Foundation under grants DGE-0903637 (an Integrative Graduate Education and Research Traineeship supporting B.R.A), IOB-0517340 (to R.W.B.) and BCS-010913 (to C.F.R.).