Investigating Passive Muscle Mechanics With Biaxial Stretch

Introduction: The passive stiffness of skeletal muscle can drastically affect muscle function in vivo, such as the case for fibrotic tissue or patients with cerebral palsy. The two constituents of skeletal muscle that dominate passive stiffness are the intracellular protein titin and the collagenous extracellular matrix (ECM). However, efforts to correlate stiffness and measurements of specific muscle constituents have been mixed, and thus the complete mechanisms for changes to muscle stiffness remain unknown. We hypothesize that biaxial stretch can provide an improved approach to evaluating passive muscle stiffness.

Methods: We performed planar biaxial materials testing of passively stretched skeletal muscle and identified three previously published datasets of uniaxial materials testing. We developed and employed a constitutive model of passive skeletal muscle that includes aligned muscle fibers and dispersed ECM collagen fibers with a bimodal von Mises distribution. Parametric modeling studies and fits to experimental data (both biaxial and previously published) were completed.

Results: Biaxial data exhibited differences in time dependent behavior based on orientation (p < 0.0001), suggesting different mechanisms supporting load in the direction of muscle fibers (longitudinal) and in the perpendicular (transverse) directions. Model parametric studies and fits to experimental data exhibited the robustness of the model (<20% error) and how differences in tissue stiffness may not be observed in uniaxial longitudinal stretch, but are apparent in biaxial stretch.

Conclusion: This work presents novel materials testing data of passively stretched skeletal muscle and use of constitutive modeling and finite element analysis to explore the interaction between stiffness, constituent variability, and applied deformation in passive skeletal muscle. The results highlight the importance of biaxial stretch in evaluating muscle stiffness and in further considering the role of ECM collagen in modulating passive muscle stiffness.

Introduction

The human body is comprised of roughly 40% skeletal muscle – the tissue that drives locomotion, enables fine movements, and provides the capability to breathe in humans and animals alike. This is due to the innate ability of skeletal muscle to generate contractile force and thus drive movement of our musculoskeletal system. While skeletal muscle is a highly adaptable and regenerative tissue (Lieber, 2010; Lieber et al., 2017), neuromuscular conditions such as cerebral palsy, sarcopenia, and damage from acute injury can severely limit the ability of skeletal muscle to function properly (Lieber, 2010). Reductions in contractile capabilities can greatly impair muscle, however more recent work has highlighted the effects of passive muscle stiffness on form and function (Lieber and Fridén, 2019).

Dramatic increases in passive muscle stiffness, for example, can be detrimental for patients with cerebral palsy in comparison to healthy persons (Chapman et al., 2016; Lieber and Fridén, 2019). It follows then that understanding what mechanism(s) and/or constituent(s) in skeletal muscle dictate stiffness is necessary to treat these conditions and prevent extreme impairment. The two constituents that are recognized as the major contributors to the tensile stiffness of passive skeletal muscle are (1) muscle fibers (cells), and (2) the collagenous extracellular matrix (ECM) that provides the hierarchical organization of skeletal muscle (Huijing, 1999; Gillies and Lieber, 2011; Brynnel et al., 2018; Meyer and Lieber, 2018). Passive muscle stiffness has a non-linear and anisotropic nature that has been shown to vary between species and different muscles (Mohammadkhah et al., 2016). It should be noted here that throughout the manuscript we use the term “stiffness” to represent the intricate non-linear, anisotropic, and variable tensile material properties of passive skeletal muscle, and not the structural property k often used in Hooke’s Law that characterizes the structural stiffness of a physical object with specific dimensions and material properties.

Uniaxial tensile testing of longitudinal (along-fiber) muscle samples are the most common approach for evaluating tensile stiffness (Calvo et al., 2010; Sato et al., 2014; Lieber and Fridén, 2019). Other efforts to characterize the anisotropy of passive muscle have employed uniaxial stretch in both the longitudinal and transverse (cross-fiber) directions (Morrow et al., 2010; Takaza et al., 2012; Mohammadkhah et al., 2016; Wheatley et al., 2016b). However, during contraction and passive stretch, force is transmitted laterally both within skeletal muscle and between muscles (Huijing, 1999; Ramaswamy et al., 2011; Maas, 2019; Csapo et al., 2020), suggesting that muscle tissue is subject to a multi-axial stress state in vivo. This is further supported by the structure of the ECM, which consists of collagen fibrils that are dispersed around the transverse plane (Purslow, 1989; Purslow and Trotter, 1994; Gillies and Lieber, 2011). These observations raise the question as to whether uniaxial stretch is thus the most appropriate in vitro experimental technique to evaluate the stiffness of passively stretched muscle, or if multi-axial materials testing may provide certain benefits.

We propose the use of a biaxial tensile deformation as a method to elucidate the passive stiffness of skeletal muscle and have developed and employed both experimental and computational efforts to this end. This method tensions both the longitudinal (along-fiber) and transverse (cross-fiber) orientations simultaneously, which may enact mechanisms that are not observable with uniaxial stretch. Finally, we have previously shown the importance of stress relaxation in modeling passive muscle stiffness (Wheatley et al., 2016a, b), thus time dependence may also provide further insight into muscle stiffness and load sharing between muscle fibers and the ECM.

We also propose the use of computational modeling – in particular finite element analysis (FEA) – to study the passive response of skeletal muscle under both uniaxial and biaxial stretch. We aim to use a continuum-level constitutive model that accounts for stiffness of muscle fibers and the ECM and can capture the variability of stress-stretch behavior that has been observed experimentally (Mohammadkhah et al., 2016). FEA provides a scalable, robust computational tool to simulate skeletal muscle behavior (Jenkyn et al., 2002; Oomens et al., 2003; Blemker et al., 2005; Böl and Reese, 2008). Previous studies include models of muscle at the tissue level (Takaza et al., 2013; Wheatley et al., 2017b), whole muscle level (Blemker et al., 2005; Böl, 2010; Wheatley et al., 2018), with idealized geometries (Jenkyn et al., 2002; Lemos et al., 2008; Chi et al., 2010), and may incorporate the observed anisotropic (Van Loocke et al., 2006; Böl et al., 2014; Pietsch et al., 2014; Mohammadkhah et al., 2016; Wheatley et al., 2016b), hyperelastic (Meyer and Lieber, 2011; Gras et al., 2012; Simms et al., 2012; Wheatley et al., 2016a), and time dependent (Van Loocke et al., 2008; Gras et al., 2013; Wheatley et al., 2016a, c) characteristics of passive skeletal muscle. Thus, FEA is well-suited as a method to explore the tissue-level mechanics of muscle tissue across experimental data sets and loading conditions.

Comprehensively, we aim to explore if experimental and computational efforts to characterize passive muscle stiffness may be enhanced by biaxial stretch by (1) performing planar biaxial materials testing on passive skeletal muscle, (2) developing and employing a robust continuum-level constitutive model of muscle that captures uniaxial and biaxial stress-stretch behavior, and (3) using such a model to explore the similarities and differences between uniaxially and biaxially stretched muscle.

Materials and Methods

Experimental Planar Biaxial Testing

Porcine hind limbs were acquired from a local abattoir on the day of sacrifice for testing. Tissue was cooled and stored at 0°C prior to testing. No live animal handling was performed by any participants in this study. A total of four animals, seven muscles, and n = 16 total samples were used for testing. The biceps femoris muscle was harvested using standard dissection scalpels. Muscles were sliced along the orientation of fibers with a custom tool that provides 10 mm spacing between the dissection top and a high-profile histology blade (Labus and Puttlitz, 2016). Each ∼10 mm thick sample was then cut into a cruciform shape with a custom cruciform press, aligning the muscle fibers with one cruciform arm (Figure 1). Sample thickness was measured with a caliper mounted on a test stand that was zeroed to the stand platform. Thickness values were recorded in five locations on each sample – in the center of the sample and toward each cruciform arm – and averaged. Mean sample thickness was 8.90 mm with a standard error of 0.29 mm across the five measurements.

FIGURE 1

Figure 1. Planar biaxial materials testing overview, with (A) cruciform geometry, (B) representative planar biaxial sample, where the white arrow denotes the longitudinal direction and the white dashed square denotes approximate DIC region of interest (ROI), and (C) experimental stress-relaxation loading protocol schematic (note that axes are not to scale).

All materials testing was performed on a planar biaxial material testing system with 50 lb (∼220 N) load cells. Samples were gripped with 25 mm pyramid grips with an initial spacing of 30 mm between grips. Samples were subject to ten equibiaxial preconditioning cycles of 10% grip-to-grip strain (3 mm) and back to zero at 0.5 Hz prior to testing (Van Ee et al., 2000). A 0.02 N equibiaxial pre-load was then applied immediately prior to testing. The testing protocol included an equibiaxial ramp of 20% nominal (grip-to-grip) strain (6 mm) at 10%/s followed by a hold until 400 s to allow for tissue stress relaxation. Samples were then subject to equibiaxial constant rate stretch at 0.1%/s nominal (grip-to-grip) strain (0.03 mm/s) until failure. Failure was manually identified post hoc in stress-time curves where significant (>∼10%) decreases in stress were observed.

Digital image correlation software (Correlated Solutions, Inc.) was used to track strain during the constant rate ramp pull in a ∼10 × 10 mm region of interest (ROI) in the center of the sample. A solid in a reference configuration X that undergoes a deformation under an external load is placed into a deformed configuration x, which is described by the deformation gradient F (Eq. 1). For a 2D problem such as a single camera digital image correlation system, F is a 2 × 2 matrix of the deformations relative to orthogonal axes (Eq. 1). From F, 2D muscle ROI stretch λ (Eq. 2) can be calculated for the longitudinal and transverse orientations (Szczesny et al., 2012). Nominal (grip-to-grip) stretch was measured directly from grip displacement. Nominal stress S was determined by dividing load cell force by the product of sample arm length (30 mm) and mean sample thickness. A linearized modulus $E=\frac{\mathrm{\Delta }S}{\mathrm{\Delta }\mathrm{\lambda }}$ was calculated from the initial and final points of the constant ramp pull data. For comparative purposes between orientations, we used nominal stress and implemented a finite element model to determine ROI Cauchy (true) stress and material properties. All stress-stretch data were averaged either over time (for stress-relaxation data) or over stretch (for constant ramp pull data) for model fitting.

$$\mathbf{\text{F}}=\frac{\partial \mathbf{\text{x}}}{\partial \mathbf{\text{X}}}=\left[\begin{array}{cc}\hfill F_{11}\hfill & \hfill F_{12}\hfill \\ \hfill F_{21}\hfill & \hfill F_{22}\hfill \end{array}\right](1)$$

$${\mathrm{\lambda }}_1=\frac{\left|\partial {\mathbf{\text{x}}}_1\right|}{\left|\partial {\mathbf{\text{X}}}_1\right|}(2)$$

Constitutive Modeling

For a 3D solid subject to an external load, the governing linear momentum balance (Newton’s second law) can be written as Eq. 3 and the governing angular momentum balance can be written as Eq. 4, where σ is the Cauchy (true) stress, ρ is the density of the solid, b is the body force vector, and a is the acceleration vector. Assuming equilibrium with negligible body force, Eq. 3 reduces to ∇⋅σ = 0. For full derivations and further reading, the reader is directed toward Holzapfel’s Non-linear Solid Mechanics (Holzapfel, 2000).

$$\nabla \cdot \mathrm{\sigma }+\mathrm{\rho }\mathbf{\text{b}}=\mathrm{\rho }\mathbf{\text{a}}(3)$$

$${\mathrm{\sigma }}_{ji}={\mathrm{\sigma }}_{ij}(4)$$

The mechanical properties of musculoskeletal soft tissues such as skeletal muscle are often modeled with a strain energy density function, including any specified intricacies such as non-linearity, anisotropy, and nearly incompressibility. To characterize the variable nature of passive muscle anisotropy and non-linearity (Mohammadkhah et al., 2016), we have employed a continuum model Ψ_tot that includes contributions of an isotropic ground matrix ${\overline{\text{Ψ}}}_{\text{iso}}$, muscle fibers ${\overline{\text{Ψ}}}_{\text{fibers}}$, the collagenous ECM ${\overline{\text{Ψ}}}_{\text{ECM}}$, and a volumetric response Ψ_vol (Eq. 5). As muscle exhibits nearly incompressible behavior (Van Loocke et al., 2006; Takaza et al., 2012), this formulation features a decoupled deviatoric response $\overline{\text{Ψ}}={\overline{\text{Ψ}}}_{\text{iso}}+{\overline{\text{Ψ}}}_{\text{fibers}}+{\overline{\text{Ψ}}}_{\text{ECM}}$ and a dilatational response Ψ_vol. Here deformation is characterized by the volume ratio J, the deviatoric right Cauchy–Green deformation tensor $\overline{\mathbf{\text{C}}}=J^{-\frac{2}{3}}{\mathbf{\text{F}}}^{\mathbf{T}}\mathbf{\text{F}}$, the first deviatoric invariant of $\overline{\mathbf{\text{C}}}$ denoted by ${\overline{I}}_1$, and the deviatoric pseudo-invariant ${\overline{I}}_4=\mathbf{\text{m}}\cdot \overline{\mathbf{\text{C}}}\cdot \mathbf{\text{m}}$ that measures the square of muscle fiber stretch whose direction is defined by the unit vector m (Holzapfel, 2000). The Cauchy (true) stress σ can then be defined as a function of the constitutive model (Eq. 6, where dev(−) is the deviatoric operator, p is hydrostatic pressure and I is the identity matrix) (Holzapfel, 2000; Maas et al., 2012). While further detail is provided below regarding specific constitutive relations, the general model employed here is an uncoupled, fiber-reinforced material with two families of fibers – aligned muscle fibers and bimodal, continually distributed ECM collagen fibers (Figure 2; Yousefi et al., 2018; Bleiler et al., 2019). This formulation attributes muscle fibers and the ECM as the main load-bearing constituents in passively stretched muscle (Smith et al., 2019).

FIGURE 2

Figure 2. Schematic of passive muscle model for extracellular matrix and muscle fibers. A representative 2D square of muscle tissue shows the longitudinal or muscle fiber direction (red) and two families of collagen fiber dispersion (green) offset from the muscle fiber direction by an angle γ.

$${\text{Ψ}}_{\text{tot}}(\overline{\mathbf{\text{C}}},J)={\overline{\text{Ψ}}}_{\text{iso}}\left({\overline{I}}_1\right)+{\overline{\text{Ψ}}}_{\text{fibers}}\left({\overline{I}}_4\right)+{\overline{\text{Ψ}}}_{\text{ECM}}\left(\overline{\mathbf{\text{C}}}\right)+{\text{Ψ}}_{\text{vol}}\left(J\right)(5)$$

$$\mathrm{\sigma }=\text{dev}\left(2J^{-1}\overline{\mathbf{\text{F}}}\frac{\partial \overline{\text{Ψ}}}{\partial \overline{C}}{\overline{\mathbf{\text{F}}}}^T\right)+p\mathbf{\text{I}}(6)$$

The isotropic ground matrix model was modeled with an uncoupled neo-Hookean strain energy density formulation, and the volumetric term with a logarithmic function (Eqs 7 and 8, where c₁ is a shear-like modulus and k is a bulk-like modulus). Due to the highly anisotropic, non-symmetric, and nearly incompressible nature of skeletal muscle (Van Loocke et al., 2006; Takaza et al., 2012; Mohammadkhah et al., 2016), c₁ was selected as a low (but non-zero) constant value for this study and k was selected to ensure near-incompressibility, as provided in Table 1.

TABLE 1

Table 1. All model material parameters and units, with fixed values provided and omitted values for parameters that were varied in this study.

$${\overline{\text{Ψ}}}_{\text{iso}}\left({\overline{I}}_1\right)=c_1\left({\overline{I}}_1-3\right)(7)$$

$${\text{Ψ}}_{\text{vol}}\left(J\right)=\frac{k}{2}{\left(\ln J\right)}^2(8)$$

The muscle fiber contribution term was defined as a power law to model non-linear stress-stretch behavior of passive muscle when stretched in the direction of muscle fibers (Eq. 9, where ξ is a modulus-like parameter and β is the power parameter) (Takaza et al., 2012; Wheatley et al., 2016b). The ECM strain energy density function defines the behavior of a continually dispersed, 3D bimodal von Mises distribution of tension-only fibers (Ateshian et al., 2009). The formulation presented here is modified from an ellipsoidal bivariate von Mises distribution to describe the anisotropic and inhomogeneous collagen fiber distribution in articular cartilage (Zimmerman and Ateshian, 2019). Due to the continually dispersed nature of the fibers, the strain energy density function is an integration over a unit sphere of volume V of the product of the distribution R(n) (where n is the orientation of the ECM collagen fibers) and the fiber constitutive law ${\overline{\text{Ψ}}}_n$ (Eq. 10). The distribution R(n) is further broken into two functions using spherical angle functions P(θ) (where θ is the azimuth angle) and Q(φ) (where φ is the declination angle) (Eq. 11). If one assumes that the ECM fibers have no directional preference in the transverse plane (perpendicular to muscle fibers), then P(θ) decomposes to the circle equation (Eq. 12). The remaining dispersion term, Q(φ) describes the ECM collagen fiber dispersion in the along-fiber plane (Figure 2) with a bimodal von Mises function that includes the primary ECM collagen fiber orientation angle γ (the angle of offset from muscle fibers to ECM collagen fibers) and a dispersion term d that characterizes the degree of alignment of ECM collagen fibers (Eq. 13). This equation also includes an integration term q(d,γ) that enforces ∫R(n)dV = 1. By varying the ECM orientation angle γ and the dispersion term d, the relative density of ECM collagen fibers can be continuously defined throughout the solid. Finally, a neo-Hookean type fiber constitutive law was used for the ECM collagen fibers (Eq. 14), with a modulus μ that is a function of the square of the collagen fiber stretch ${\overline{I}}_n$ (FEBio User Manual 2. 8, 2018). These formulations also assume that fibers (both muscle fibers and ECM collagen fibers) can only sustain tension, not compression (Weiss et al., 1996).

$${\overline{\text{Ψ}}}_{\text{fibers}}\left({\overline{I}}_4\right)=\frac{\mathit{\text{ξ}}}{\mathrm{\beta }}{\left({\overline{I}}_4-1\right)}^{\mathrm{\beta }}(9)$$

$${\overline{\text{Ψ}}}_{\text{ECM}}\left(\mathbf{C}\right)=\int R\left(\mathbf{n}\right){\overline{\text{Ψ}}}_n\left({\overline{I}}_n\right)dV(10)$$

$$\int R\left(\mathbf{n}\right)dV=1=\int P\left(\mathrm{\theta }\right)Q\left(\mathrm{\varphi }\right)dV(11)$$

$$P\left(\mathrm{\theta }\right)={\left[\left({\cos }^2\mathrm{\theta }+{\sin }^2\mathrm{\theta }\right)\right]}^{-1/2}(12)$$

$$Q\left(\mathrm{\varphi }\right)=\frac{1}{q(d,\mathrm{\gamma })}\left\{\exp \left[2d{\cos }^2\left(\mathrm{\varphi }+\mathrm{\gamma }\right)\right]+\exp \left[2d{\cos }^2\left(\mathrm{\varphi }-\mathrm{\gamma }\right)\right]\right\}(13)$$

$${\overline{\text{Ψ}}}_n\left({\overline{I}}_n\right)=\frac{\mathrm{\mu }}{4}{\left({\overline{I}}_n-1\right)}^2(14)$$

A quasi-linear Prony series viscoelastic formulation was used to model stress-relaxation of passively stretched skeletal muscle (Wheatley et al., 2016a). Briefly, the deviatoric stress $\overline{\mathrm{\sigma }}$ can be defined as a function of a convolution integral (Eq. 15, where G(t) is the relaxation function, t is time, and ζ is an integration variable). A Prony series relaxation function (Eq. 16) enables the use of viscoelastic coefficients g_i and associated time constants τ_i that characterize the amount and rate of relaxation, respectively. For this study, we fixed τ_i terms as spaced parameters to ensure a broad range of relaxation rates (Table 1) and varied g_i terms (Vaidya and Wheatley, 2019).

$$\overline{\mathrm{\sigma }}\left(t\right)=\underset{-\mathrm{\infty }}{\overset{t}{\int }}G\left(t-\mathrm{\zeta }\right)\frac{d\overline{\mathrm{\sigma }}}{d\mathrm{\zeta }}d\mathrm{\zeta }(15)$$

$$G\left(t\right)=1+\sum _{i=1}^4{g}_i\text{exp}\left(-\frac{t}{{\mathrm{\tau }}_i}\right)(16)$$

Finite Element Modeling

All finite element modeling results presented here were conducted using the open source finite element package FEBio (Maas et al., 2012). A custom plugin was written to apply the bivariate von Mises distribution of the ECM collagen fibers R(n). To model biaxial stretch, a symmetric, eighth cruciform finite element model consisting of 2,184 linear hexahedral elements was developed (Figure 3B). This model was chosen to represent a cruciform 30 mm × 30 mm in width and height and a thickness of 8.9 mm. In addition to symmetric boundary conditions (Figure 3B), the face of each cruciform arm was fixed to a rigid body and subject to displacements to mimic the experimental protocol (Figure 3C). Reaction force divided by initial arm cross sectional area was calculated as nominal stress and total model length was divided by initial model length to determine nominal stretch. By using the solid mixture capabilities in FEBio, separate viscoelastic parameters were assigned to the stress contributions from the muscle fibers and ECM.

FIGURE 3

Figure 3. (A) Representative color contour plot of longitudinal (horizontal) stretch from digital image correlation. (B) Deformed symmetric cruciform finite element model with model ROI (black dashed line) symmetric boundary conditions (note that Z-symm is on the bottom, hidden face), and color contour plot of longitudinal stretch. (C) Deformed symmetric cruciform finite element model, showing initial cross-sectional area and reaction force at the model boundary (grip locations) used to determine model nominal stress. (D) Undeformed and deformed single element finite element model. Model Cauchy stress was directly output from the single element.

Constitutive parameters were optimized to experimental data by fitting model nominal stress to experimental nominal stress. Model ROI stretch was calculated based on the position of model ROI surface nodes similar to experimental stretch (Eqs 1 and 2), then used as a validation to experimental DIC stretch (Figures 3A,B). Parameter optimization was completed in two steps – first the viscoelastic Prony series parameters were fit to normalized stress-relaxation data for both the longitudinal and transverse directions, then hyperelastic parameters were fit to the full set of longitudinal and transverse stress data. This approach has the advantage of reducing the overall number of parameters needed to be optimized at any given step in the process by first determining stress relaxation behavior and then hyperelastic stiffness (Vaidya and Wheatley, 2019). All optimization was performed in MATLAB using constrained non-linear optimization (lsqnonlin) by varying model parameters and minimizing the sum of squared residuals between model (σ^m) and experimental (σ^e) stresses as an objective function obj across all experimental data points (total number npts) (Eq. 17). Nominal stress was used for fitting of our cruciform finite element model (Figure 3C) to planar biaxial experimental data as Cauchy (true) stress cannot be estimated from experimental planar biaxial tests without a correction factor, which is typically determined from FEA.

$$obj=\sum _{i=1}^{npts}{\left({\mathrm{\sigma }}_i^e-{\mathrm{\sigma }}_i^m\right)}^2(17)$$

For comparisons across previously published experimental studies of uniaxial stretch of passive skeletal muscle, a simplified approach of a single linear hexahedral finite element model was implemented (Figure 3D). Three previously published studies of materials testing of skeletal muscle under uniaxial tension were identified (Table 2). These studies provide a range of data across species, muscles, and orientations for model comparison and fitting. The finite element model was fit to the experimental studies by comparing literature Cauchy (true) stress to model Cauchy stress as a function of directional stretch. Data from Wheatley et al. (2016b) were zeroed following an initial stress-relaxation phase for consistency with other data. Following all model fitting, optimized parameters for each data set were used to simulate stress-stretch behavior under uniaxial and equibiaxial stretch. Finally, a simple parameter study was conducted to further highlight the differences in stress-stretch behavior between uniaxial and biaxial stretch. Cauchy stress was used for fitting to previously published uniaxial data and for parametric studies as it requires fewer assumptions to estimate than under biaxial conditions, Cauchy stress is reported in the literature cited here, and the use of Cauchy stress allows for a simplified single element finite element model. No stress conversions were calculated directly from a push forward or pull back operation in this work, as all planar biaxial model fitting used nominal stress only, and all uniaxial and parametric studies used Cauchy stress only.

TABLE 2

Table 2. Summary of data of passively stretched skeletal muscle used in this study.

To summarize, we fit the cruciform finite element model to planar biaxial data, then used the optimized parameters to simulate uniaxial stretch. Conversely, we fit the simplified finite element model to uniaxial data from three different previously published studies, then used the optimized parameters to simulate biaxial stretch.

Statistics

Stress relaxation data were normalized to sample peak stress and fit to a power law model (Eq. 18, where σ_n is normalized stress, t is relaxation time, and a and b are power law coefficients) to characterize the rate of relaxation between orientations. The power law b coefficients (rate of relaxation), stress at three time points – peak stress, end of relaxation, and end of ramp pull, and linearized modulus from the pull phase were compared between orientations using a paired t-test. A linear regression was performed to investigate the potential effect of post-mortem time on modulus for both directions. For all statistical tests, significance was set at p < 0.05.

$${\mathrm{\sigma }}_n=a{t}^b(18)$$

Model fits to experimental data were evaluated with an average percent error for each experimental data point, normalized root mean square error (NRMSE, where 1 is a perfect fit and −∞ is the worst possible fit), and root mean square error (RMSE, in kPa) (Vaidya and Wheatley, 2019).

Results

Experimental Planar Biaxial Data

Biaxial data showed that longitudinal direction nominal stress was greater and decreased at a faster rate during stress relaxation than transverse direction stress. This was supported by both visual analysis of normalized relaxation (Figure 4A) and statistical analysis (Figure 4B). Specifically, the paired t-tests suggest that the power law b coefficient was greater in the longitudinal orientation (p < 0.0001). Stress was greater at the peak (p = 0.021), end of relaxation phase (p = 0.037), and end of constant rate pull phase (p = 0.0063), and the linearized modulus was greater in the longitudinal direction versus the transverse direction (p = 0.028). The power law fits provided excellent agreement to experimental data visually (Figure 4A) and with mean R² values of 0.985 and 0.974 for longitudinal and transverse data, respectively. Linear regression results showed that modulus was not correlated with post-mortem time (p > 0.6, R² < 0.02 for both directions).

FIGURE 4

Figure 4. (A) Normalized stress relaxation data (shown as dashed mean curves and standard error bars) and a power law fit to the mean data (shown as solid curves). Note that power law fits and experimental data are visually overlapping and thus nearly indistinguishable. (B) Bar graphs for mean power law b coefficient, stress data, and linearized modulus with standard error bars.

Model Fitting

The use of constrained non-linear optimization produced a strong fit of the biaxial finite element model to experimental data, both for the stress relaxation phase as well as the constant ramp pull phase. This is shown both visually (Figures 5A,B,D) and through statistical analysis (NRMSE > 0.9, Table 3). Measured experimental stretch in the sample ROI from digital image correlation and predicted model ROI stretch are provided for model validation (Figures 5C,D). The model showed strong agreement to transverse stretch data, and overpredicted longitudinal stretch somewhat.

FIGURE 5

Figure 5. Experimental and model data of passive muscle subject to planar biaxial stretch. (A) Nominal stress relaxation step data (open circles with standard error bars) and model fits (solid curves), with the first 10 s of these data for clarity shown at right, (B) longitudinal constant rate nominal stress-stretch curves for all experimental samples (thin curves) and model (thick curve), (C) longitudinal constant rate ROI stretch-nominal stretch curves, (D) transverse constant rate nominal stress-stretch curves, and (E) transverse constant rate ROI stretch-nominal stretch curves.

TABLE 3

Table 3. Statistical fitting results between model and experimental stress-stretch data.

The single element uniaxial model exhibited strong fitting capabilities across all uniaxial stretch data sets as observed visually (Figure 6) and by evaluating the statistical differences between model outputs and experimental data (NRMSE > 0.85, Table 3). Specifically, the model was able to match a wide range of anisotropy and non-linearity between data sets, including directions of greatest stiffness of the transverse direction (Takaza et al., 2012; Wheatley et al., 2016b) and 45° (Mohammadkhah et al., 2016).

FIGURE 6

Figure 6. Modeling fits to uniaxial tensile experimental data from the previously published works of (A) Wheatley et al. (2016b), (B) Takaza et al. (2012), and (C) Mohammadkhah et al. (2016). Experimental data are shown as open circles and model data are solid curves.

Optimized hyperelastic parameter values (Table 4) demonstrate the variability of passive muscle material properties. For the constitutive model used here, there was a particularly wide range of muscle fiber stiffness ξ (4–110 kPa), ECM modulus μ (28–1,700 kPa), and ECM orientation angles γ (32–90°). It should be noted here that the Mohammadkhah et al. (2016) chicken data best fit produced a negligible muscle fiber modulus, hence the reported value of 0 kPa. Finally, to ensure a unique set of parameters for each optimal fit, the ECM fiber dispersion parameter d was fixed for some of the optimizations, as shown in Table 4. Optimized viscoelastic parameter values (Table 5) further highlight the differences in viscoelastic behavior between orientations, as muscle fiber g_i terms were larger than those applied to the ECM term. This shows greater relaxation for the muscle fiber term in comparison to the ECM term.

TABLE 4

Table 4. Optimized parameter values from fits to various experimental data sets.

TABLE 5

Table 5. Optimized viscoelastic parameters g_i and associated time constants τ_i from fits to planar biaxial experimental data.

Modeling Biaxial and Uniaxial Stretch

Simulating both uniaxial stretch and biaxial stretch with each optimized parameter set showed different effects of biaxial stretch on model response (Figure 7 and Table 6). Specifically, stress-stretch curves from Wheatley et al. (2016b). Parameters were largely unaffected by biaxial versus uniaxial stretch, while biaxial stretch greatly increased stiffness for the Mohammadkhah et al. (2016). Parameter set. A parametric study of uniaxial and biaxial stretch for two different parameter sets – one with Aligned fibers and one with Dispersed fibers – shows the models exhibit nearly identical stiffness behavior under uniaxial stretch (Figure 8A) but distinctly different behavior under biaxial stretch (up to 119% difference, Figure 8B). This was observed for both the longitudinal and transverse directions, highlighting the role of the dispersed ECM fibers and assumptions of anisotropy in altering model behavior. Parameter values used (Table 7) fall within those optimized to experimental data (Table 4).

FIGURE 7

Figure 7. Stress-stretch curves for simulated uniaxial stretch (dashed) and biaxial stretch (solid) for optimized parameters from (A) biaxial data presented in this study, (B) Wheatley et al. (2016b), and (C) Mohammadkhah et al. (2016). The increase in stress-stretch curve stiffness with biaxial stretch versus uniaxial stretch is denoted with arrows. Note that some models predict a negligible increase in stiffness (Wheatley) and others a major increase in stiffness (Mohammadkhah).

TABLE 6

Table 6. Differences in Cauchy stress at λ = 1.3 between uniaxial and biaxial stretch conditions for all optimized parameter sets.

FIGURE 8

Figure 8. Parametric study stress-stretch curves for (A) uniaxial stretch and (B) biaxial stretch. Note that for both models (Aligned – solid curves and Dispersed – dashed curves) they exhibit nearly identical uniaxial behavior for both longitudinal (back curves) and transverse (blue curves), but distinctly different behavior when subject to biaxial stretch, with percentage differences between Aligned and Dispersed shown at λ = 1.3.

TABLE 7

Table 7. Parameter values for the models shown in Figure 8.

Discussion

Planar Biaxial Testing

We have presented here, to the best of our knowledge, the first set of experimental data of passive skeletal muscle subject to planar biaxial stretch. From these data, we determined that the porcine hind limb tissue tested exhibited the following characteristics: (1) faster relaxation for longitudinal samples (Figures 4A,B) and (2) greater stiffness for the longitudinal direction (Figure 4B). For the viscoelastic response, we previously measured greater relaxation rate with greater longitudinal stiffness of muscle tissue (Wheatley et al., 2016b), which agrees with the results seen here. The differences in relaxation rate between orientations suggest two mechanisms that support load in passively stretched skeletal muscle. This observation is supported by ongoing efforts that have shown that passive muscle stiffness in mammals is dictated by both the collagenous ECM (Meyer and Lieber, 2011, 2018) and muscle fibers themselves (Brynnel et al., 2018). Here we suggest that both constituents may contribute to longitudinal stiffness, and that measuring the anisotropic viscoelastic response may further elucidate load sharing between constituents.

It is known that muscle exhibits stress relaxation at both the muscle fiber level (Meyer et al., 2011; Rehorn et al., 2014) and the whole muscle or tissue level (Best et al., 1994; Gras et al., 2013; Wheatley et al., 2016a). Meyer et al. (2011) measured ∼95% stress relaxation when stretching muscle fibers at 2,000%/s and ∼80% stress relaxation at 20%/s, which exceeds stress relaxation observed in highly collagenous tissues such as tendon (Atkinson et al., 1999). Based on these findings and the observation of less relaxation in the transverse direction from our data, the ECM may exhibit less stress-relaxation than muscle fibers. This requires further experimental efforts to confirm or deny, however. To the best of our knowledge, there have been no studies that have compared viscoelastic behavior between single fiber and tissue level samples or have tried to measure the viscoelastic properties of the ECM directly or indirectly. Such a study would help contextualize tissue-level measurements of longitudinal and transverse viscoelastic behavior in regards to the contribution of muscle fibers and the ECM to tissue stiffness.

In comparison to previously published data, most studies of passively stretched muscle have observed a greater stiffness in the transverse direction in comparison to the longitudinal direction, albeit to varying degrees (Takaza et al., 2012; Mohammadkhah et al., 2016; Wheatley et al., 2016b). One previous study observed greater stiffness in the longitudinal direction (Morrow et al., 2010). These comprehensively suggest that anisotropy may be variable in skeletal muscle and may depend on a range of physiological factors. Exploring the link between anisotropy and in vivo function was outside the scope of this work but would be appropriate for future studies. We have previously hypothesized that a greater longitudinal stiffness was the result of rigor mortis (Wheatley et al., 2016b), however in this work all testing in this study was completed within seven hours to reduce this risk (Van Ee et al., 2000; Van Loocke et al., 2006) and tissue stiffness was not correlated with post-mortem testing time (p > 0.1, R² < 0.2). It is thus unlikely that our data are driven by rigor mortis alone. Use of a relaxing agent (Meyer and Lieber, 2018) could be used to further prevent the effects of post-mortem stiffening. Nonetheless, the data presented here should not be viewed as a comprehensive set of muscle material properties, but as a validation of an experimental and computational technique to investigate muscle stiffness.

Constitutive Modeling of Experimental Data

Fitting results show the capability of our constitutive model to accurately simulate the range of experimentally observed anisotropic and non-linear stress-stretch behavior. This is shown both visually (Figures 5, 6) and through statistical evaluation (Table 3). We also used the experimental biaxial ROI stretch data for model validation (Figure 5). These results comprehensively suggest that our model is well-suited for studying the tissue-level mechanics of passively stretched skeletal muscle. To encourage a unique solution for each data set, the ECM fiber dispersion parameter d was fixed based on a qualitative comparison to muscle ECM fiber dispersion (Purslow and Trotter, 1994). Additionally, dataset characteristics such as viscoelasticity (biaxial data), data at 45° (Mohammadkhah and Takaza), and non-linear longitudinal data (Wheatley) enforced a unique solution for each optimization.

The constitutive model used in this study includes a non-linear muscle fiber term ${\overline{\text{Ψ}}}_{\text{fibers}}\left({\overline{I}}_4\right)$ which is a function of ${\overline{I}}_4$, the square of muscle fiber stretch. We chose to employ a power law function for this term as it models the non-linear stress-stretch response of muscle in the longitudinal orientation (Mohammadkhah et al., 2016; Wheatley et al., 2016a) with only two parameters. The optimized values for the modulus-like parameter ξ (0–108 kPa) and for the power coefficient β (2.35–3.78) are reasonable, although ξ = 0 for the Mohammadkhah data is questionable. However, Mohammadkhah et al., 2016 data was obtained from chicken pectoralis muscle tissue, which as they note has a higher collagen content (Nishimura, 2010). This may partially explain why our model optimization approach identified a negligible muscle fiber modulus for this data set if collagen is dominating the stress-stretch response. Our remaining muscle fiber modulus-like values of 4.2–108 kPa compare reasonably to experimental observations of ∼40 kPa in mice (Meyer and Lieber, 2018).

Our model of muscle ECM ${\overline{\text{Ψ}}}_{\text{ECM}}\left(\overline{\mathbf{\text{C}}}\right)$ describes the collagen fibers with a neo-Hookean hyperelastic model (with a shear modulus μ) and a bimodal von Mises distribution (with angle γ and dispersion d). Our use of a single modulus term is a simplification of a highly complex combination of ECM collagen amount, type, crosslinking, and crimp. While these each have been studied in regards to tissue stiffness through either experimentation (Smith and Barton, 2014; Chapman et al., 2015; Mohammadkhah et al., 2018; Lieber and Fridén, 2019; Smith et al., 2019) or modeling (Gindre et al., 2013; Bleiler et al., 2019; Spyrou et al., 2019; Valentin and Simms, 2020), developing a unique set of parameters from healthy tissue-level data that incorporates each of these was outside of the scope of this work. This also does not address the different layers of ECM structure such as perimysium and endomysium. We instead chose to use the approach of minimizing the number of model parameters while ensuring a strong fit to experimental data.

Purslow and Trotter (1994) measured muscle ECM collagen fiber orientations under a range of physiological conditions and found that the primary fiber alignment angle was dependent on stretch, but ranged from ∼20–80°. Qualitatively, recent mammalian ECM scanning electron microscopy by Sleboda et al. (2020) found that multilayered, collagen-rich ECM was common between a range of species but that microstructure was less consistent. These studies suggest that ECM fiber angle and dispersion may vary with a range of mechanical, anatomical, and physiological factors such as animal size and muscle fiber type distributions. The optimized ECM fiber angles we determined (32–90°) thus seem reasonable.

In considering specific modeling studies relevant to this work, Yucesoy et al. (2002) modeled the muscle fibers and ECM as distinct but linked constituents. Gindre et al. (2013) developed a microstructural model of a muscle fiber wrapped with a single family of dispersed ECM fibers to explore titin and ECM contributions. Yousefi et al. (2018) showed how a similar model of the ECM as the only load bearing mechanism with two perfectly reinforcing fiber directions can describe the observed anisotropy in passively stretched bovine, porcine, and chicken muscle. Bleiler et al. (2019) designed and formulated a passive constitutive model with dispersed collagen fibers surrounding muscle fibers that could be integrated into a finite element simulation. Teklemariam et al. (2019) used a similar micromechanical approach with distinct muscle fiber and ECM domains. Spyrou et al. (2019) developed a multiscale model that employed homogenization from a microstructurally derived model to a continuum-level response.

While each of these studies present advantages for modeling the passive response of skeletal muscle, we have chosen to use a similar approach to Yousefi et al. (2018) with the extension of the model to include ECM collagen dispersion and muscle fiber stiffness. After applying assumptions for the low-stiffness isotropic ground matrix (Wheatley et al., 2017a) and near-incompressibility (Takaza et al., 2012), this model required five parameters to describe the hyperelastic response – two for the muscle fibers (stiffness and non-linearity) and three for the ECM (stiffness, direction, and dispersion). The advantage of this approach is a relatively low number of parameters while still enabling model robustness. The use of a Prony series viscoelastic model may increase the overall number of parameters of the model, but as we have shown in this and previous works (Vaidya and Wheatley, 2019), those parameters can be optimized with a two-step fitting procedure. Based on stress-stretch data alone, it would be unclear how load is shared between the ECM and muscle fibers. However, the stress-relaxation data shows distinct time dependent differences between longitudinal and transverse stress relaxation rate (Figure 4). This suggests load may be supported by both muscle fibers and ECM, and perhaps more so the muscle fibers in the longitudinal direction.

It should be noted that the model chosen here enables a wide range of stress-stretch behavior and is generally informed by muscle physiology, but is not derived from microstructure and does not account for effects of interaction between the ECM and muscle fibers. The parameters (such as ECM fiber angle and dispersion) may be generally related to tissue microstructure, but are not direct analogs. One must be careful not to conclude concrete microstructural findings based on the fitting results presented here.

Modeling Uniaxial Versus Biaxial Stretch

Expanding our modeling from fitting to simulations of uniaxial versus biaxial stretch showed variability between data sets (Figure 7 and Table 7). Generally speaking, materials exhibit greater stiffness when stretched biaxially versus uniaxially. However, for highly anisotropic materials with multiple families of fibers, the effect may not be as dramatic as expected, as shown in the uniaxial versus biaxial comparisons of the Wheatley et al. (2016b) parameter set (Figure 7B and Table 6). In this case, the model ECM fibers are aligned perpendicular to muscle fibers and have low dispersion and during biaxial stretch each set of fibers are recruited independently. Conversely, the Mohammadkhah et al. (2016) parameter set increased in excess of 100% in both the longitudinal and transverse directions (Figure 7C and Table 6). Here the ECM fibers are oriented between directions and highly dispersed, which recruits these fibers during both longitudinal and transverse stretch. Thus, the biaxial deformation will stretch the ECM fibers to a greater amount.

The potential physiological relevance of a case where biaxial stretch and uniaxial stretch exhibit similar stress-stretch behavior can be seen in Figure 8 and Table 7. Here we have identified two sets of parameters that fall within the previously optimized values that have nearly indistinguishable uniaxial stress-stretch behavior in both the longitudinal and transverse orientations (Figure 8A). When subject to biaxial stretch however, Dispersed shows drastic changes in stiffness while Aligned is largely unaffected (113% difference in the longitudinal direction between models). This presents a simplified case where two muscles that may seem to have the same mechanical properties when stretched uniaxially would in fact have quite different mechanics when subject to a more complex deformation. In effect, these differences are “hidden” by uniaxial stretch. This could partially explain that differences in longitudinal stiffness between cerebral palsy and healthy muscle cannot be explained by collagen content, quantity, and cross-linking alone (Chapman et al., 2015; Lieber and Fridén, 2019; Smith et al., 2019).

Smith et al. (2019) discuss these collagen content-passive stiffness correlations and a relatively minor contribution of collagen crosslinking that this observation “… suggests the intriguingly possibility that higher-order structures may determine tissue stiffness to a greater extent than molecular components.” We suggest that ECM collagen fiber orientation and dispersion may be these “higher-order structures” and show with our model how differences in tissue stiffness could be hidden by uniaxial stretch (Figure 8A). As noted above, the technique employed here is a continuum-level hyperelastic constitutive model. We do not imply that this model is a direct prediction of tissue microstructure, only that our model has shown robust and accurate stress-stretch behavior and that similar mechanisms may be present.

Another consideration for uniaxial versus biaxial stretch is the observation of transverse load transmission in contracting muscle as well as laterally between individual muscles (Huijing, 1999; Yucesoy et al., 2008). If load generated longitudinally by muscle fibers is transmitted transversally through the ECM, then muscle tissue must be subject to a multi-axial stress state in vivo. Our parametric study suggests that biaxial stretch could enact a stiffening effect to longitudinal stress-stretch response in comparison to uniaxial only (Figures 8A,B). For tissues that exhibit higher load sharing of the ECM, this effect could be exaggerated, and biaxial stretch would in effect increase the perceived tissue stiffness, and thus perhaps increase the efficiency of load transfer in vivo during contraction. However, further experimental research is needed to confirm if this is the case for biaxially stretched skeletal muscle. Nonetheless, we have highlighted the importance of a biaxial deformation in passively stretched skeletal muscle, and hope that this consideration can drive future work to better understand load transmission in vitro and in vivo.

Limitations and Future Directions

This work is not without limitations. Firstly, the geometry selection of a simplified, idealized cruciform or single element cuboid is clearly not a representation of the geometric/structural complexities of whole, in vivo skeletal muscle. However, the experimental data used in this study are generated from tissue samples isolated from whole muscle, and thus do not represent a full in vivo environment either. This isolation is necessary to accurately determine tensile material properties. While a more detailed set of geometric finite element models could be developed to match average specimen geometry from each uniaxial experiment, this may not necessarily yield improvements in fit or different study conclusions. The advantage of our geometric approach is in computational efficiency and simplicity – as the optimization protocol that fit the model to experimental stress-stretch curves does not require significant computation time and is highly stable. Nonetheless, as experimental and finite element models of in vivo muscle deformation have shown complex strains (Blemker and Delp, 2005; Böl et al., 2015), use and validation of this model in such cases would be a significant benefit to the field.

While our constitutive model exhibits robustness in simulating tensile stress-stretch behavior (Figure 6), it does not model microstructural and physiological characteristics such as collagen crosslinking, multiple collagen types, or muscle fiber-ECM interactions. Including such components would likely yield increased robustness and physiological accuracy of such a model. However, our model has exhibited efficacy in simulating a wide range of passive muscle stretch. We have shown that this model can inform future studies of ECM structure – such as collagen fiber orientation and dispersion – while fitting tissue-level data and maintaining experimental observations such as near-incompressibility.

It should be noted that model validation across uniaxial and biaxial stretch would strengthen future applications of this model. Additionally, experiments such as biaxial materials testing coupled with decellularization or muscle fiber isolation would provide necessary insight into the extent to which this model or future improved models can accurately characterize load sharing between muscle fibers and the ECM. This would greatly strengthen this work, and provide strong efficacy for application of this model to in vivo conditions of muscle impairment such as cerebral palsy. We have also not explored the model response under compression or during active contraction as those are outside the scope of this work. Thus, this model should be viewed not as a comprehensive model of passive skeletal muscle, but as an effective tool in better understanding passive muscle stiffness.

Conclusion

Based on section “Results and Discussion” of this work, we have made the following observations, recommendations, and conclusions:

(1) We performed biaxial stress-relaxation testing on passive skeletal muscle and suggest that this approach can be used to effectively characterize passive muscle mechanics.

(2) Our model of a dispersed ECM contribution and aligned muscle fibers was able to exhibit broad variability in simulating and fitting tensile stiffness, non-linearity, and anisotropy of passive skeletal muscle.

(3) This model, in conjunction with experimental data, exhibited the role of biaxial stretch in measuring passive muscle stiffness and suggesting future work to explore inconsistent correlations between muscle ECM collagen measurements and passive stiffness.

Future validation, development, and employment of modeling and biaxial experimentation would elucidate the role of the ECM in in vivo muscle function, and help explain how detrimental changes to muscle stiffness – such as those observed in cerebral palsy – may be explained by ECM structure.

Data Availability Statement

The datasets generated for this study are available on request to the corresponding author.

Author Contributions

BW contributed to the study design, experimentation and modeling, data processing, and manuscript development.

Funding

This material is based upon work supported by the National Science Foundation under Grant No. 1828082.

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

The author would like to acknowledge Brandon K. Zimmerman for his assistance adapting his collagen fiber dispersion model for this work. The author also would like to acknowledge Kristen Fu for her work developing biaxial protocols and assistance processing digital image correlation images.

References

Ateshian, G. A., Rajan, V., Chahine, N. O., Canal, C. E., and Hung, C. T. (2009). Modeling the matrix of articular cartilage using a continuous fiber angular distribution predicts many observed phenomena. J. Biomech. Eng. 131:061003. doi: 10.1115/1.3118773

PubMed Abstract | CrossRef Full Text | Google Scholar

Atkinson, T. S., Ewers, B. J., and Haut, R. C. (1999). The tensile and stress relaxation responses of human patellar tendon varies with specimen cross-sectional area. J. Biomech. 32, 907–914. doi: 10.1016/S0021-9290(99)00089-5

CrossRef Full Text | Google Scholar

Best, T. M., McElhaney, J., Garrett, W. E., and Myers, B. S. (1994). Characterization of the passive responses of live skeletal muscle using the quasi-linear theory of viscoelasticity. J. Biomech. 27, 413–419. doi: 10.1016/0021-9290(94)90017-5

CrossRef Full Text | Google Scholar

Bleiler, C., Ponte Castañeda, P., and Röhrle, O. (2019). A microstructurally-based, multi-scale, continuum-mechanical model for the passive behaviour of skeletal muscle tissue. J. Mech. Behav. Biomed. Mater. 97, 171–186. doi: 10.1016/J.JMBBM.2019.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Blemker, S. S., and Delp, S. L. (2005). Three-dimensional representation of complex muscle architectures and geometries. Ann. Biomed. Eng. 33, 661–673. doi: 10.1007/s10439-005-1433-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Blemker, S. S., Pinsky, P. M., and Delp, S. L. (2005). A 3D model of muscle reveals the causes of nonuniform strains in the biceps brachii. J. Biomech. 38, 657–665. doi: 10.1016/j.jbiomech.2004.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Böl, M. (2010). Micromechanical modelling of skeletal muscles: from the single fibre to the whole muscle. Arch. Appl. Mech. 80, 557–567. doi: 10.1007/s00419-009-0378-y

CrossRef Full Text | Google Scholar

Böl, M., Ehret, A. E., Leichsenring, K., Weichert, C., and Kruse, R. (2014). On the anisotropy of skeletal muscle tissue under compression. Acta Biomater. 10, 3225–3234. doi: 10.1016/j.actbio.2014.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Böl, M., Leichsenring, K., Ernst, M., Wick, C., Blickhan, R., and Siebert, T. (2015). Novel microstructural findings in M. plantaris and their impact during active and passive loading at the macro level. J. Mech. Behav. Biomed. Mater. 51, 25–39. doi: 10.1016/j.jmbbm.2015.06.026

PubMed Abstract | CrossRef Full Text | Google Scholar

Böl, M., and Reese, S. (2008). Micromechanical modelling of skeletal muscles based on the finite element method. Comput. Methods Biomech. Biomed. Engin. 11, 489–504. doi: 10.1080/10255840701771750

PubMed Abstract | CrossRef Full Text | Google Scholar

Brynnel, A., Hernandez, Y., Kiss, B., Lindqvist, J., Adler, M., Kolb, J., et al. (2018). Downsizing the molecular spring of the giant protein titin reveals that skeletal muscle titin determines passive stiffness and drives longitudinal hypertrophy. eLife 7:e40532. doi: 10.7554/eLife.40532

PubMed Abstract | CrossRef Full Text | Google Scholar

Calvo, B., Ramírez, A., Alonso, A., Grasa, J., Soteras, F., Osta, R., et al. (2010). Passive nonlinear elastic behaviour of skeletal muscle: experimental results and model formulation. J. Biomech. 43, 318–325. doi: 10.1016/j.jbiomech.2009.08.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapman, M. A., Meza, R., and Lieber, R. L. (2016). Skeletal muscle fibroblasts in health and disease. Differentiation 92, 108–115. doi: 10.1016/j.diff.2016.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Chapman, M. A., Pichika, R., and Lieber, R. L. (2015). Collagen crosslinking does not dictate stiffness in a transgenic mouse model of skeletal muscle fibrosis. J. Biomech. 48, 375–378. doi: 10.1016/j.jbiomech.2014.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Chi, S., Hodgson, J., Chen, J., Reggie Edgerton, V., Shin, D. D., Roiz, R. A., et al. (2010). Finite element modeling reveals complex strain mechanics in the aponeuroses of contracting skeletal muscle. J. Biomech. 43, 1243–1250. doi: 10.1016/j.jbiomech.2010.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Csapo, R., Gumpenberger, M., and Wessner, B. (2020). Skeletal Muscle Extracellular Matrix – What Do We Know About Its Composition, Regulation, and Physiological Roles? A Narrative Review. Front. Physiol. 11:253. doi: 10.3389/fphys.2020.00253

PubMed Abstract | CrossRef Full Text | Google Scholar

FEBio User Manual, 2. 8 (2018). Available online at: https://help.febio.org/FEBio/FEBio_um_2_8/ (accessed January 16, 2020).

Google Scholar

Gillies, A. R., and Lieber, R. L. (2011). Structure and function of the skeletal muscle extracellular matrix. Muscle Nerve 44, 318–331. doi: 10.1002/mus.22094

PubMed Abstract | CrossRef Full Text | Google Scholar

Gindre, J., Takaza, M., Moerman, K. M., and Simms, C. K. (2013). A structural model of passive skeletal muscle shows two reinforcement processes in resisting deformation. J. Mech. Behav. Biomed. Mater. 22, 84–94. doi: 10.1016/j.jmbbm.2013.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Gras, L. L., Mitton, D., Viot, P., and Laporte, S. (2012). Hyper-elastic properties of the human sternocleidomastoideus muscle in tension. J. Mech. Behav. Biomed. Mater. 15, 131–140. doi: 10.1016/j.jmbbm.2012.06.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Gras, L. L., Mitton, D., Viot, P., and Laporte, S. (2013). Viscoelastic properties of the human sternocleidomastoideus muscle of aged women in relaxation. J. Mech. Behav. Biomed. Mater. 27, 77–83. doi: 10.1016/j.jmbbm.2013.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Holzapfel, G. A. (2000). Nonlinear Solid Mechanics. Chichester: Wiley.

Google Scholar

Huijing, P. A. (1999). Muscle as a collagen fiber reinforced composite: a review of force transmission in muscle and whole limb. J. Biomech. 32, 329–345. doi: 10.1016/S0021-9290(98)00186-9

CrossRef Full Text | Google Scholar

Jenkyn, T., Koopman, B., Huijing, P. A., Lieber, R. L., and Kaufman, K. R. (2002). Finite element model of intramuscular pressure during isometric contraction of skeletal muscle. Phys. Med. Biol. 47, 4043–4061. doi: 10.1088/0031-9155/47/22/309

CrossRef Full Text | Google Scholar

Labus, K. M., and Puttlitz, C. M. (2016). Viscoelasticity of brain corpus callosum in biaxial tension. J. Mech. Phys. Solids 96, 591–604. doi: 10.1016/j.jmps.2016.08.010

CrossRef Full Text | Google Scholar

Lemos, R. R., Epstein, M., and Herzog, W. (2008). Modeling of skeletal muscle: the influence of tendon and aponeuroses compliance on the force-length relationship. Med. Biol. Eng. Comput. 46, 23–32. doi: 10.1007/s11517-007-0259-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Lieber, R. L. (2010). Skeletal Muscle Structure, Function, and Plasticity. Philadelphia, PA: Lippincott Williams and Wilkins.

Google Scholar

Lieber, R. L., and Fridén, J. (2019). Muscle contracture and passive mechanics in cerebral palsy. J. Appl. Physiol. 126, 1492–1501. doi: 10.1152/japplphysiol.00278.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Lieber, R. L., Roberts, T. J., Blemker, S. S., Lee, S. S. M., and Herzog, W. (2017). Skeletal muscle mechanics, energetics and plasticity. J. Neuroeng. Rehabil. 14:108. doi: 10.1186/s12984-017-0318-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Maas, H. (2019). Significance of epimuscular myofascial force transmission under passive muscle conditions. J. Appl. Physiol. 126, 1465–1473. doi: 10.1152/japplphysiol.00631.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Maas, S. A., Ellis, B. J., Ateshian, G. A., and Weiss, J. A. (2012). FEBio: finite elements for biomechanics. J. Biomech. Eng. 134:011005. doi: 10.1115/1.4005694

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, G., and Lieber, R. L. (2018). Muscle fibers bear a larger fraction of passive muscle tension in frogs compared with mice. J. Exp. Biol. 221:jeb.182089. doi: 10.1242/jeb.182089

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, G. A., and Lieber, R. L. (2011). Elucidation of extracellular matrix mechanics from muscle fibers and fiber bundles. J. Biomech. 44, 771–773. doi: 10.1016/j.jbiomech.2010.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Meyer, G. A., McCulloch, A. D., and Lieber, R. L. (2011). A Nonlinear Model of Passive Muscle Viscosity. J. Biomech. Eng. 133:091007. doi: 10.1115/1.4004993

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohammadkhah, M., Murphy, P., and Simms, C. K. (2016). The in vitro passive elastic response of chicken pectoralis muscle to applied tensile and compressive deformation. J. Mech. Behav. Biomed. Mater. 62, 468–480. doi: 10.1016/j.jmbbm.2016.05.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Mohammadkhah, M., Murphy, P., and Simms, C. K. (2018). Collagen fibril organization in chicken and porcine skeletal muscle perimysium under applied tension and compression. J. Mech. Behav. Biomed. Mater. 77, 734–744. doi: 10.1016/j.jmbbm.2017.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrow, D. A., Haut Donahue, T. L., Odegard, G. M., and Kaufman, K. R. (2010). Transversely isotropic tensile material properties of skeletal muscle tissue. J. Mech. Behav. Biomed. Mater. 3, 124–129. doi: 10.1016/j.jmbbm.2009.03.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishimura, T. (2010). The role of intramuscular connective tissue in meat texture. Anim. Sci. J. 81, 21–27. doi: 10.1111/j.1740-0929.2009.00696.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Oomens, C. W. J., Maenhout, M., van Oijen, C. H., Drost, M. R., and Baaijens, F. P. (2003). Finite element modelling of contracting skeletal muscle. Philos. Trans. R. Soc. Lond. B. Biol. Sci. 358, 1453–1460. doi: 10.1098/rstb.2003.1345

PubMed Abstract | CrossRef Full Text | Google Scholar

Pietsch, R., Wheatley, B. B. B., Haut Donahue, T. L., Gilbrech, R., Prabhu, R., Liao, J., et al. (2014). Anisotropic compressive properties of passive porcine muscle tissue. J. Biomech. Eng. 136:111003. doi: 10.1115/1.4028088

PubMed Abstract | CrossRef Full Text | Google Scholar

Purslow, P. P. (1989). Strain-induced reorientation of an intramuscular connective tissue network: implications for passive muscle elasticity. J. Biomech. 22, 21–31. doi: 10.1016/0021-9290(89)90181-4

CrossRef Full Text | Google Scholar

Purslow, P. P., and Trotter, J. A. (1994). The morphology and mechanical properties of endomysium in series-fibred muscles: variations with muscle length. J. Muscle Res. Cell Motil. 15, 299–308.

Google Scholar

Ramaswamy, K. S., Palmer, M. L., van der Meulen, J. H., Renoux, A., Kostrominova, T. Y., Michele, D. E., et al. (2011). Lateral transmission of force is impaired in skeletal muscles of dystrophic mice and very old rats. J. Physiol. 589, 1195–1208. doi: 10.1113/jphysiol.2010.201921

PubMed Abstract | CrossRef Full Text | Google Scholar

Rehorn, M. R., Schroer, A. K., and Blemker, S. S. (2014). The passive properties of muscle fibers are velocity dependent. J. Biomech. 47, 687–693. doi: 10.1016/j.jbiomech.2013.11.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Sato, E. J., Killian, M. L., Choi, A. J., Lin, E., Esparza, M. C., Galatz, L. M., et al. (2014). Skeletal muscle fibrosis and stiffness increase after rotator cuff tendon injury and neuromuscular compromise in a rat model. J. Orthop. Res. 32, 1111–1116. doi: 10.1002/jor.22646

PubMed Abstract | CrossRef Full Text | Google Scholar

Simms, C. K., Van Loocke, M., and Lyons, C. G. (2012). Skeletal Muscle in Compression: Modeling Approaches for the Passive Muscle Bulk. Int. J. Multiscale Comput. Eng. 10, 143–154. doi: 10.1615/IntJMultCompEng.2011002419

PubMed Abstract | CrossRef Full Text | Google Scholar

Sleboda, D. A., Stover, K. K., and Roberts, T. J. (2020). Diversity of extracellular matrix morphology in vertebrate skeletal muscle. J. Morphol. 281, 160–169. doi: 10.1002/jmor.21088

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, L. R., and Barton, E. R. (2014). Collagen content does not alter the passive mechanical properties of fibrotic skeletal muscle in mdx mice. Am. J. Physiol. Physiol. 306, C889–C898. doi: 10.1152/ajpcell.00383.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, L. R., Pichika, R., Meza, R. C., Gillies, A. R., Baliki, M. N., Chambers, H. G., et al. (2019). Contribution of extracellular matrix components to the stiffness of skeletal muscle contractures in patients with cerebral palsy. Connect. Tissue Res. doi: 10.1080/03008207.2019.1694011 [Online ahead of print]

PubMed Abstract | CrossRef Full Text | Google Scholar

Spyrou, L. A., Brisard, S., and Danas, K. (2019). Multiscale modeling of skeletal muscle tissues based on analytical and numerical homogenization. J. Mech. Behav. Biomed. Mater. 92, 97–117. doi: 10.1016/j.jmbbm.2018.12.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Szczesny, S. E., Peloquin, J. M., Cortes, D. H., Kadlowec, J. A., Soslowsky, L. J., and Elliott, D. M. (2012). Biaxial tensile testing and constitutive modeling of human supraspinatus tendon. J. Biomech. Eng. 134:021004. doi: 10.1115/1.4005852

PubMed Abstract | CrossRef Full Text | Google Scholar

Takaza, M., Moerman, K. M., Gindre, J., Lyons, G., and Simms, C. K. (2012). The anisotropic mechanical behaviour of passive skeletal muscle tissue subjected to large tensile strain. J. Mech. Behav. Biomed. Mater. 17, 209–220. doi: 10.1016/j.jmbbm.2012.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Takaza, M., Moerman, K. M., and Simms, C. K. (2013). Passive skeletal muscle response to impact loading: experimental testing and inverse modelling. J. Mech. Behav. Biomed. Mater. 27, 214–225. doi: 10.1016/j.jmbbm.2013.04.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Teklemariam, A., Hodson-Tole, E., Reeves, N. D., and Cooper, G. (2019). A micromechanical muscle model for determining the impact of motor unit fiber clustering on force transmission in aging skeletal muscle. Biomech. Model. Mechanobiol. 18, 1401–1413. doi: 10.1007/s10237-019-01152-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Vaidya, A. J., and Wheatley, B. B. (2019). An experimental and computational investigation of the effects of volumetric boundary conditions on the compressive mechanics of passive skeletal muscle. J. Mech. Behav. Biomed. Mater. 2019:103526. doi: 10.1016/j.jmbbm.2019.103526

PubMed Abstract | CrossRef Full Text | Google Scholar

Valentin, T., and Simms, C. (2020). An inverse model of the mechanical response of passive skeletal muscle: implications for microstructure. J. Biomech. 99:109483. doi: 10.1016/j.jbiomech.2019.109483

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Ee, C. A., Chasse, A. L., and Myers, B. S. (2000). Quantifying skeletal muscle properties in cadaveric test specimens: effects of mechanical loading, postmortem time, and freezer storage. J. Biomech. Eng. 122, 9–14. doi: 10.1115/1.429621

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Loocke, M., Lyons, C. G., and Simms, C. K. (2006). A validated model of passive muscle in compression. J. Biomech. 39, 2999–3009. doi: 10.1016/j.jbiomech.2005.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Loocke, M., Lyons, C. G., and Simms, C. K. (2008). Viscoelastic properties of passive skeletal muscle in compression: stress-relaxation behaviour and constitutive modelling. J. Biomech. 41, 1555–1566. doi: 10.1016/j.jbiomech.2008.02.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Weiss, J. A., Maker, B. N., and Govindjee, S. (1996). Finite element implementation of incompressible, transversely isotropic hyperelasticity. Comput. Methods Appl. Mech. Eng. 135, 107–128. doi: 10.1016/0045-7825(96)01035-1033

CrossRef Full Text | Google Scholar

Wheatley, B. B., Odegard, G. M., Kaufman, K. R., and Haut Donahue, T. L. (2018). Modeling skeletal muscle stress and intramuscular pressure: a whole muscle active-passive approach. J. Biomech. Eng. 140:081006. doi: 10.1115/1.4040318

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheatley, B. B., Morrow, D. A., Odegard, G. M. G. M., Kaufman, K. R. K. R., and Haut Donahue, T. L. T. L. (2016a). Skeletal muscle tensile strain dependence: hyperviscoelastic nonlinearity. J. Mech. Behav. Biomed. Mater. 53, 445–454. doi: 10.1016/j.jmbbm.2015.08.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheatley, B. B., Odegard, G. M., Kaufman, K. R., and Donahue, T. L. H. (2016b). How does tissue preparation affect skeletal muscle transverse isotropy? J. Biomech. 49, 3056–3060. doi: 10.1016/j.jbiomech.2016.06.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheatley, B. B., Pietsch, R. B., Haut Donahue, T. L., and Williams, L. N. (2016c). Fully non-linear hyper-viscoelastic modeling of skeletal muscle in compression. Comput. Methods Biomech. Biomed. Engin. 19, 1181–1189. doi: 10.1080/10255842.2015.1118468

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheatley, B. B., Odegard, G. M., Kaufman, K. R., and Haut Donahue, T. L. (2017a). A validated model of passive skeletal muscle to predict force and intramuscular pressure. Biomech. Model. Mechanobiol. 16, 1011–1022. doi: 10.1007/s10237-016-0869-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wheatley, B. B., Odegard, G. M., Kaufman, K. R., and Haut Donahue, T. L. (2017b). A case for poroelasticity in skeletal muscle finite element analysis: experiment and modeling. Comput. Methods Biomech. Biomed. Engin. 20, 598–601. doi: 10.1080/10255842.2016.1268132

PubMed Abstract | CrossRef Full Text | Google Scholar

Yousefi, A.-A. A. K., Nazari, M. A., Perrier, P., Panahi, M. S., and Payan, Y. (2018). A new model of passive muscle tissue integrating Collagen Fibers: consequences for muscle behavior analysis. J. Mech. Behav. Biomed. Mater. 88, 29–40. doi: 10.1016/j.jmbbm.2018.07.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Yucesoy, C. A., Koopman, B. H. F. J. M., Grootenboer, H. J., and Huijing, P. A. (2008). Extramuscular myofascial force transmission alters substantially the acute effects of surgical aponeurotomy: assessment by finite element modeling. Biomech. Model. Mechanobiol. 7, 175–189. doi: 10.1007/s10237-007-0084-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yucesoy, C. A., Koopman, B. H. F. J. M. F. J. M., Huijing, P. A., and Grootenboer, H. J. (2002). Three-dimensional finite element modeling of skeletal muscle using a two-domain approach: linked fiber-matrix mesh model. J. Biomech. 35, 1253–1262. doi: 10.1016/S0021-9290(02)00069-6

CrossRef Full Text | Google Scholar

Zimmerman, B. K., and Ateshian, G. A. (2019). “A three dimensional rotationally nonsymmetric continuous fiber distribution for articular cartilage,” in Proceedings of the 16th International Symposium on Computer Methods in Biomechanics and Biomedical Engineering, New York, NY.

Google Scholar