Modeling of active skeletal muscles: a 3D continuum approach incorporating multiple muscle interactions

1 Introduction

As a major type of human muscle, skeletal muscle provides locomotion, maintains posture, and stabilizes bones and joints (Fung, 1993). A single skeletal muscle is an organ with a hierarchical structure, composed of complex fibrous actuators organized in different levels of aggregation (Gotti et al., 2020). Each skeletal muscle is wrapped in the epimysium (or deep fascia), an irregular fibrous connective tissue sheath which separates individual muscle and allows the enveloped muscle to move independently. As such, skeletal muscle behaviors are highly heterogeneous, and the muscle properties are extremely complicated (Blemker, 2017). Due to the structural, architectural, and mechanical complexity of the skeletal muscles, in vivo characterization of some physical quantities cannot be made without highly invasive methods (Dao and Tho, 2018). Computational modeling of skeletal muscles has become an indispensable modern method to determine these quantities, understand the complex heterogeneous behaviors of skeletal muscles, and the role of muscle pathology in movement biomechanics.

Phenomenological and biophysical approaches have been developed to model active muscle, including the most widely used Hill-type (Hill, 1938) and Huxley-type (Huxley, 1957) muscle models. The biophysical Huxley-type model for muscle contraction considered the cross-bridges and dynamics of myosin filaments within muscle. Usually, it needs a detailed set of experimental data to calibrate the model, and it is computationally expensive (Khodaei et al., 2013). Recent advances in multiscale biophysical muscle modeling have enabled simulation from the microscopic half-sarcomere level to the whole-muscle organ level (Röhrle et al., 2012; Heidlauf et al., 2016; Spyrou et al., 2019). As the most widely used landmark model to study mechanical behaviors of muscles, the Hill-type model (Hill, 1938) uses a macroscopic phenomenological approach to describe active muscle properties by introducing a contractile component with force–length and force–velocity relations for the entire muscle in one dimension (Zajac, 1989). To enhance the one-dimensional (1D) Hill-type model’s biofidelity and to represent the shape of muscles more realistically, it has been extended into three-dimensional (3D) Hill-type models for musculotendon dynamics. Most of the developed Hill-type phenomenological models for active muscles included passive and active components (Dao and Tho, 2018). The passive component has been commonly characterized by hyperelastic materials using the Neo-Hookean model (Göktepe et al., 2014), the Mooney-Rivlin models (Wu et al., 2014), or the Ogden model (Hedenstierna et al., 2008; Mo et al., 2019). The sophisticated Ogden model can accurately describe the nonlinear response of highly deformable tissues but is computationally expensive and requires more material parameters than the other two models. The Neo-Hookean type model, which is a popular choice, was adopted in the current study for its computational efficiency and ease of implementation but it can be limited in its applicability to tissues with highly deformable and compressible states. The definition of the active component of muscles requires the utilization of hierarchical fibers and their activation mechanism, which can be modeled by various approaches.

Over the past 20 years, there have been many exciting developments in constitutive formulations for active muscle behaviors. Based on the model of the cardiac muscle proposed by Humphrey and Yin (1987), Martins et al. (1998 and 2006) modeled skeletal muscles as a fiber-reinforced incompressible hyperelastic composite. Fernandez et al. (2005) presented a finite element (FE) framework for modelling the mechanical response of the rectus femoris muscle to a flexion loading. However, the simulations in that study simplified several aspects including the spatial distribution of material properties and excitation-contraction coupling. A number of studies have modeled the skeletal muscles by a combination of Hill-type 1D truss elements (or beam elements) for the active part and 3D solid elements for the passive part, i.e., discrete approaches for modeling skeletal muscles. For instance, Hedenstierna et al. (2008) modeled active muscles using a combination of passive non-linear, hyperelastic viscoelastic solid elements and active Hill-type truss elements, which were validated by rabbit leg muscles with simple geometry; Feller et al. (2016) incorporated a strain-dependent muscle activation scheme using 1D beam elements into the neck of a full human body model for a low-loading condition through a head fall test. Aiming to predict injury responses of lower limb-pelvis during the whole joint movement; Mo et al. (2018) developed a lower extremity model with active muscles in a recent study. The muscles were modeled as a 3D passive unit with an isotropic hyperelastic material and Hill-type 1D truss elements, which were created surrounding the 3D mesh of passive unit. They also implemented a muscle control strategy for combining multibody dynamics and finite element analysis (FEA) (Mo et al., 2019). In these studies, the distribution of the multiple truss or beam elements needs to be defined based on the mesh of solid elements which was a tedious process and can cause numerical instability issues.

There have been several advancements in the constitutive modeling for simulating complex active skeletal muscle behaviors at the macroscopic scale, including tissue strain, intramuscular pressure, muscle force, fatigue, and parameter identification (Blemker et al., 2005; Grasa et al., 2011; Khodaei et al., 2013; Lu et al., 2010 and 2011a; Tang et al., 2007 and 2009; Wheatley et al., 2018; Klotz et al., 2021). For example, Blemker et al. (2005) created a 3D FE model of the biceps brachii and compared the predicted tissue strains with experimental data. Lu et al. (2010) developed a transversely visco-hyperelastic model of skeletal muscle tissue under high strain rates. The model was validated using experimental studies of rabbit tibialis anterior (TA) muscle with simplified cylindrical geometry and a uniform muscle fiber direction. A similar study has been presented by Khodaei et al. (2013) using a transversely isotropic visco-hyperelastic model. Wheatley et al. (2018) presented an FE model of rabbit TA with a hyperporo-viscoelastic constitutive law to predict intramuscular pressure and muscle force under passive and active conditions.

In addition to the complexity of constitutive representations of skeletal muscles, the definition of fiber trajectories or architecture is a critical step to provide important input to the muscle model. There are several different approaches to characterize the arrangement of fibers, such as parallel-fibered architecture using parallel fiber distribution along a uniform direction or at some specific angle (Blemker et al., 2005; Clemen et al., 2017; Lu et al., 2010 and 2011a; Tang et al., 2007 and 2009), bipennate fiber orientation (Fernandez et al., 2005), and fusiform-shaped fiber distribution (Grasa et al., 2011). Other approaches include defining a fiber map template which is morphable to the target muscle geometry (Blemker and Delp, 2005; Röhrle and Pullan, 2007), using a method with non-uniform rational B-splines (NURBS) to characterize the fiber orientation (Lu et al., 2011b). However, these approaches have their own limitations. For instance, the standard fiber map templates may not be adequate to characterize various muscle architectures (Blemker, 2017), and alternative approaches either rely on oversimplified fiber orientations (Blemker et al., 2005; Tang et al., 2007 and 2009; Fernandez et al., 2005; Grasa et al., 2011) or encounter challenges with characterizing muscles that have multiple branches (Lu et al., 2011b). This study employed a physics-based method that utilizes computational fluid dynamics (CFD) to determine the muscle fiber architecture of muscles with realistic geometry. Inouye et al. (2015) recently proposed this approach and Handsfield et al. (2017) demonstrated it to be effective in predicting subject-specific muscle fiber patterns.

While many computational approaches have been developed to model the complex mechanical behaviors of individual active muscles, few studies have incorporated multiple muscles into a system-level model (Dao and Tho, 2018). For example, Elyasi et al. (2022) combined FE and multibody models into a lower limb model and to predict surgery’s impact on patellar position, but the model has limitations due to simplified anatomical details and material properties (e.g., all muscles were modelled as a passive state). Similarly, Spyrou and Aravas (2012) developed an FE model of the lower leg with active muscles to simulate plantar flexion but did not account for muscle interactions or detailed anatomical relationships. To date, very few studies have been able to model the complex anatomical relationships between organs, such as muscle-to-muscle and muscle-to-bone, in a biofidelic manner for a body model with multiple muscles.

This work presents the development of an efficient framework for modeling 3D active skeletal muscles at a system level to simulate realistic joint motion (see Supplementary Material for the flowchart of implementation process). The muscles were modeled as active, quasi-incompressible, transversely isotropic and hyperelastic solids, which were implemented into the commercial FE software ABAQUS/Explicit through user-defined material subroutines (VUMAT). A parametric study was conducted to validate the muscle model with ideal 3D geometry at the tissue level, followed by the modeling of the quadriceps muscle group combined with a previously described detailed FE representation of the human knee (Ali et al., 2016; Hume et al., 2018). To demonstrate the effectiveness of the active muscle model in the lower limb, the system level model performance was demonstrated by simulating knee flexion-extension with active quadriceps in a lower limb model. Computational fluid dynamics (CFD) was employed to specify fiber directions for the quadriceps muscles. The interactions between adjacent muscles, and between muscles and femur were innovatively simulated using the cohesive contact approach in ABAQUS. This allows sliding between parts under shearing but restricts separation between the involved parts under tension.

2 Methods and materials

2.1 Constitutive modeling of skeletal muscle

The constitutive model is the core of the physics of modeling active skeletal muscle. The muscle constitutive law used in this work was based on the strain energy approach, and the general framework of the Hill-type hyperelastic active muscle model was derived from previous studies (Blemker et al., 2005; Tang et al., 2009; Lu et al., 2011a; Khodaei et al., 2013; Kleinbach et al., 2017). The muscle was modeled as a fiber-reinforced composite (Figure 1A) in which the strain energy density function includes the contributions from the composite and the contributions from the fibers (Blemker, 2017). The strain energy density function of active muscle can be expressed by an additive form:

$W=W_{\mathrm{I}}+W_{\mathrm{f}}+W_{\mathrm{V}}(1)$

where $W_{\mathrm{I}}$, $W_{\mathrm{f}}$, and $W_{\mathrm{V}}$ represent the strain energy densities for the muscle matrix, the muscle fiber and related volume change, which can be given by the three equations as follows (Tang et al., 2009)

$W_{\mathrm{I}}=c\left\{\exp \left[b\left({\overline{I}}_1^C-3\right)-1\right]\right\}(2)$

in which the ${\overline{I}}_1^C$ is the first invariant of the right Cauchy-Green strain tensor ($\mathbf{C}$), b and c are material parameters associated with the matrix (Chen and Zeltzer, 1992)

$W_{\mathrm{f}}={\int }_1^{{\overline{\lambda }}_{\mathrm{f}}}\left[{\sigma }_{\mathrm{P}\mathrm{E}}\left(\lambda \right)+{\sigma }_{\mathrm{S}\mathrm{E}\mathrm{E}}\left(\lambda ,{\lambda }_{\mathrm{s}}\right)\right]d\lambda (3)$

where ${\sigma }_{\mathrm{P}\mathrm{E}}\left(\lambda \right)$ signifies the stress contributed by the parallel element (PE), ${\sigma }_{\mathrm{S}\mathrm{E}\mathrm{E}}\left(\lambda ,{\lambda }_{\mathrm{s}}\right)$ is the stress contributed by the series elastic element (SEE), λ refers to the fiber stretch ratio, ${\overline{\lambda }}_{\mathrm{f}}$ represents the fiber stretch ratio with the volume change eliminated, i.e., modified fiber stretch ratio, and ${\lambda }_{\mathrm{s}}$ represents the stretch ratio in the SEE.

$W_{\mathrm{V}}=\frac{1}{D}{\left(J-1\right)}^2(4)$

where $D$ is the constant of compressibility, and $J$ denotes the Jacobian of the deformation gradient.

FIGURE 1

FIGURE 1. The three-element Hill’s muscle model: (A) illustration of the model as a fiber-reinforced composite, and (B) the schematic diagram of the model consisting of three elements [a parallel element (PE), a series elastic element (SEE) and a contractile element (CE)].

The stress ${}^{t+\Delta t}{\sigma }_{\mathrm{P}\mathrm{E}}$ produced by the parallel element (PE) of the muscle at time $t+\Delta t$ can be expressed by

${}^{t+\Delta t}{\sigma }_{\mathrm{P}\mathrm{E}}={\sigma }_0{f}_{\mathrm{P}\mathrm{E}}\left({}^{t+\Delta t}{\overline{\lambda }}_{\mathrm{f}}\right)(5)$

where ${\sigma }_0$ is the maximum isometric stress, and the function $f_{\mathrm{P}\mathrm{E}}\left({}^{t+\Delta t}{\overline{\lambda }}_{\mathrm{f}}\right)$ can be defined by (Chen and Zeltzer, 1992)

$f_{\mathrm{P}\mathrm{E}}\left({}^{t+\Delta t}{\overline{\lambda }}_{\mathrm{f}}\right)=\left\{\begin{array}{c}A{\left({}^{t+\Delta t}{\overline{\lambda }}_{\mathrm{f}}-1\right)}^2,\\ 0,\end{array}\right.\begin{array}{l}\begin{array}{cc}\mathrm{f}\mathrm{o}\mathrm{r}& {}^{t+\Delta t}{\overline{\mathrm{\lambda }}}_{\mathrm{f}}\end{array}>1\\ \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}(6)$

in which $A$ is the parameter related to the stress in the PE.

Accordingly, the stress produced by the SEE at time $t+\Delta t$ is specified as (Pinto and Fung, 1973; Fung, 1993; Tang et al., 2009)

${}^{t+\Delta t}{\mathrm{\sigma }}_{\mathrm{S}\mathrm{E}\mathrm{E}}=\beta \left[{\mathrm{e}}^{\alpha \left({}^{t+\Delta t}{\lambda }_{\mathrm{s}}-1\right)}-1\right]={\mathrm{e}}^{\alpha \Delta {\lambda }_{\mathrm{s}}}\beta {\mathrm{e}}^{\alpha \left({}^t{\lambda }_{\mathrm{s}}-1\right)}-\beta ={\mathrm{e}}^{\alpha \Delta {\lambda }_{\mathrm{s}}}\left({}^t{\sigma }_{\mathrm{S}\mathrm{E}\mathrm{E}}+\beta \right)-\beta (7)$

with

${}^t{\mathrm{\sigma }}_{\mathrm{S}\mathrm{E}\mathrm{E}}=\beta \left[{\mathrm{e}}^{\alpha \left({}^t{\lambda }_{\mathrm{s}}-1\right)}-1\right](8)$

where $\alpha$ and $\beta$ are the material constants associated with the stress in the PE, and the relation ${}^{t+\Delta t}{\lambda }_{\mathrm{s}}={}^t{\lambda }_{\mathrm{s}}+\Delta {\lambda }_{\mathrm{s}}$ was used in the derivation process for Eq. 7.

The stress produced in the contractile element (CE) is defined as

${}^{t+\Delta t}{\sigma }_{\mathrm{C}\mathrm{E}}={\sigma }_0\cdot f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)\cdot f_{\mathrm{v}}\left({\dot{\lambda }}_{\mathrm{m}}\right)\cdot f_{\mathrm{t}}\left(t+\Delta t\right)(9)$

where the $f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)$, $f_{\mathrm{v}}\left({\dot{\lambda }}_{\mathrm{m}}\right)$, and $f_{\mathrm{t}}\left(t+\Delta t\right)$ denote the function of muscle stretch, the function of muscle velocity, and the muscle activation function, respectively. Their mathematical expressions are given by the equations below (Gordon et al., 1966; Böl and Reese, 2008):

$f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)=\left\{\begin{array}{l}0,\\ 9{\left(r_{\mathrm{f}}^{\mathrm{\lambda }}-0.4\right)}^2,\\ 1-4{\left(1-r_{\mathrm{f}}^{\mathrm{\lambda }}\right)}^2,\\ 9{\left(r_{\mathrm{f}}^{\mathrm{\lambda }}-1.6\right)}^2,\\ 0,\end{array}\right.\begin{array}{l}\begin{array}{cc}\mathrm{i}\mathrm{f}& r_{\mathrm{f}}^{\mathrm{\lambda }}<0.4\end{array}\\ \begin{array}{cc}\mathrm{i}\mathrm{f}& 0.4\le r_{\mathrm{f}}^{\mathrm{\lambda }}<0.6\end{array}\\ \begin{array}{cc}\mathrm{i}\mathrm{f}& 0.6\le r_{\mathrm{f}}^{\mathrm{\lambda }}<1.4\end{array}\\ \begin{array}{cc}\mathrm{i}\mathrm{f}& 1.4\le r_{\mathrm{f}}^{\mathrm{\lambda }}<1.6\end{array}\\ \begin{array}{cc}\mathrm{i}\mathrm{f}& r_{\mathrm{f}}^{\mathrm{\lambda }}\ge 1.6\end{array}\end{array}(10)$

where the dimensionless constant $r_{\mathrm{f}}^{\mathrm{\lambda }}={}^t{\overline{\lambda }}_{\mathrm{f}}/{\lambda }_{\mathrm{o}\mathrm{p}\mathrm{t}}$, and ${\lambda }_{\mathrm{o}\mathrm{p}\mathrm{t}}$ is the optimal fiber stretch at which the sarcomeres reach optimal length (Blemker et al., 2005).

It is generally known that the force decreases with increasing velocity during concentric muscle contraction (Hill, 1970; Fung, 1993). The hyperbolic function of muscle velocity $f_{\mathrm{v}}\left({\dot{\lambda }}_{\mathrm{m}}\right)$ can be expressed in terms of the ratio of stretch rates $r_{\mathrm{m}}^{\mathrm{\lambda }}={\dot{\lambda }}_{\mathrm{m}}/{\dot{\lambda }}_{\mathrm{m}}^{\min }$ as (Böl and Reese, 2008)

$f_{\mathrm{v}}\left({\dot{\lambda }}_{\mathrm{m}}\right)=\left\{\begin{array}{l}\frac{1-r_{\mathrm{m}}^{\mathrm{\lambda }}}{1+k_{\mathrm{c}}{r}_{\mathrm{m}}^{\mathrm{\lambda }}},\\ d-\left(d-1\right)\frac{1+r_{\mathrm{m}}^{\mathrm{\lambda }}}{1-k_{\mathrm{c}}{k}_{\mathrm{e}}{r}_{\mathrm{m}}^{\mathrm{\lambda }}},\end{array}\right.\begin{array}{l}\begin{array}{cc}\mathrm{i}\mathrm{f}& {\dot{\lambda }}_{\mathrm{m}}\le 0\end{array}\\ \\ \begin{array}{cc}\mathrm{i}\mathrm{f}& {\dot{\lambda }}_{\mathrm{m}}>0\end{array}\end{array}(11)$

where ${\dot{\lambda }}_{\mathrm{m}}$ is stretch rate in the CE, ${\dot{\lambda }}_{\mathrm{m}}^{\min }$ stands for the minimum stretch rate, $k_{\mathrm{c}}$ and $k_{\mathrm{e}}$ are dimensionless parameters controlling the curvature of the concentric phase (shortening) (Hill, 1938) and eccentric phase (lengthening) of the curve, and $d$ is a dimensionless constant to describe the offset of the eccentric function.

Muscle activation behavior is complicated and further investigation of activation patterns is needed (Ting et al., 2012; Roux et al., 2021). In this work, an exponential function has been employed to characterize muscle activation (Meier and Blickhan, 2000):

$f_{\mathrm{t}}\left(t\right)=\left\{\begin{array}{l}{n}_1,\\ n_1+\left(n_2-n_1\right)\cdot h_{\mathrm{t}}\left(t,t_0\right),\\ n_1+\left(n_2-n_1\right)\cdot h_{\mathrm{t}}\left(t_1,t_0\right)-\left[\left(n_2-n_1\right)\cdot h_{\mathrm{t}}\left(t_1,t_0\right)\right]\cdot h_{\mathrm{t}}\left(t,t_1\right),\end{array}\right.\begin{array}{l}\begin{array}{cc}\mathrm{i}\mathrm{f}& t\le t_0\end{array}\\ \begin{array}{cc}\mathrm{i}\mathrm{f}& t_0<t\le t_1\end{array}\\ \begin{array}{cc}\mathrm{i}\mathrm{f}& t>t_1\end{array}\end{array}(12)$

with

$h_{\mathrm{t}}\left(t_{\mathrm{j}},t_{\mathrm{k}}\right)=\left\{1-\exp \left[-S\cdot \left(t_{\mathrm{j}}-t_{\mathrm{k}}\right)\right]\right\}(13)$

where $t_0$ and $t_1$ denote the activation time and the deactivation time of the muscle, separately, $n_1$ represents the activation level before and after the activation, $n_2$ denotes the activation level during the activation process, and $S$ is the exponential factor (Lu et al., 2010).

Since the stresses produced in the SEE and CE are identical (Figure 1B), the unknown $\Delta {\lambda }_{\mathrm{s}}$ can be obtained by Equations 7 and (9), which leads to the governing equation to solve the stress increment of SEE (Tang et al., 2009):

$f\left(\Delta {\lambda }_{\mathrm{s}}\right)=\left(w_2+w_3\Delta {\lambda }_{\mathrm{s}}\right)e^{\alpha \Delta {\lambda }_{\mathrm{s}}}-w_4\Delta {\lambda }_{\mathrm{s}}-w_5=0(14)$

where for muscle shortening (concentric phase)

$w_2=\left(\beta +{}^t{\sigma }_{\mathrm{s}}\right)\left(1+\frac{k_{\mathrm{c}}\cdot w_1}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}\right)(15)$

$w_3=-\left(\beta +{}^t{\sigma }_{\mathrm{s}}\right)\frac{k\cdot k_{\mathrm{c}}}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}(16)$

$w_4=\frac{-\beta \cdot k_{\mathrm{c}}+f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)\cdot f_{\mathrm{t}}\left(t+\Delta t\right)}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}k(17)$

$w_5=\beta +f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)\cdot f_{\mathrm{t}}\left(t+\Delta t\right)-\frac{f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)\cdot f_{\mathrm{t}}\left(t+\Delta t\right)-\beta \cdot k_{\mathrm{c}}}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}{w}_1(18)$

with

$w_1=\left(1+k\right){}^{t+\Delta t}{\lambda }_{\mathrm{f}}-{}^t{\lambda }_{\mathrm{m}}-k^t{\lambda }_{\mathrm{s}}(19)$

and for muscle lengthening (eccentric phase)

$w_2=\left(\beta +{}^t{\sigma }_{\mathrm{s}}\right)\left(1-\frac{k_{\mathrm{c}}\cdot k_{\mathrm{e}}\cdot {\alpha }_1}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}\right)(20)$

$w_3=\left(\beta +{}^t{\sigma }_{\mathrm{s}}\right)\frac{k\cdot k_{\mathrm{c}}\cdot k_{\mathrm{e}}}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}(21)$

$w_4=\frac{\beta \cdot k_{\mathrm{c}}\cdot k_{\mathrm{e}}+f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)\cdot f_{\mathrm{t}}\left(t+\Delta t\right)\cdot \left(d\cdot k_{\mathrm{c}}\cdot k_{\mathrm{e}}+d-1\right)}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}k(22)$

$w_5=\beta +{}^t{f}_{\mathrm{\lambda }}\left({\overline{\lambda }}_{\mathrm{f}}\right)\cdot f_{\mathrm{t}}\left(t+\Delta t\right)-\frac{f_{\mathrm{\lambda }}\left({}^t{\overline{\lambda }}_{\mathrm{f}}\right)\cdot f_{\mathrm{t}}\left(t+\Delta t\right)\cdot \left[1-d\left(1+k_{\mathrm{c}}\cdot k_{\mathrm{e}}\right)\right]-\beta \cdot k_{\mathrm{c}}\cdot k_{\mathrm{e}}}{{\dot{\lambda }}_{\mathrm{m}}^{\min }\cdot \Delta t}{w}_1(23)$

where $k$ denotes the ratio of the length of CE to that of SEE, which is set as 0.3 (Lu et al., 2011a).

2.2 Parametric study and validation by muscle elongation simulation

The above constitutive model was implemented in user-defined subroutines (VUMAT) using ABAQUS/Explicit, which contains more than 14 material parameters (Table 1). The identification and calibration of these material parameters is a significant aspect of developing the credible muscle model. Based on elongation simulations using both passive and active muscle constitutive models with idealized 3D geometry, the parametric study including investigation of the parameter sensitivity was performed to calibrate and validate the muscle model.

TABLE 1

TABLE 1. Material parameters used for the muscle model.

The FE model of muscle with fusiform-shaped geometry (Figure 2A), which was adopted in previous studies (Böl and Reese, 2008; Lu et al., 2011a) was created in ABAQUS to simulate the muscle response under uniaxial extension. It has a length of 50 mm, and the diameters of the end section (minimum area) and mid-belly cross-section (maximum area) were set as 9 and 17.5 mm, respectively. The fiber orientation was assumed to be aligned along a uniform direction parallel to the axis of the model, since the fiber dispersion is minor in the fusiform skeletal muscles (Zajac, 1989). The simulations were carried out with two different scenarios: elongation of passive muscle (Figure 2B) and elongation of active muscle (Figure 2C). In the passive simulation, one of the ends of the muscle model was fixed and the other end was pulled with a constant velocity (5 mm/s) for 2 s. In the active simulation, the dimension and mesh were the same as the passive muscle elongation. While the muscle was also fixed on one end, the muscle was held constant in length for 0.5 s to ensure full activation, then the muscle was pulled with a constant velocity of 5 mm/s for 2 s (the muscle was activated during this process). In a similar fashion to the parameter study conducted in one of our authors’ previous work (Lu et al., 2011a), several parameters were determined through parameter studies. These parameters include D, b, c, d, $k_{\mathrm{c}}$ and $k_{\mathrm{e}}$, as demonstrated in Table 1. The simulated engineering stress response was collected and compared with the available experimental data. In the passive muscle elongation, the stress-strain response was compared with the experimental data reported by Davis et al. (2003), which was measured from the isolated rabbit TA muscle. The active muscle elongation response was compared with the stress-strain response of the active rabbit TA from experiments reported by Myers et al. (1998).

FIGURE 2

FIGURE 2. The FE model of muscle with fusiform-shaped geometry: (A) initial configuration (B) a deformed configuration of passive muscle (20% elongation), and (C) a deformed configuration of active muscle (20% elongation) and 100% activation.

2.3 Estimation of the muscle fiber architecture by CFD

In a CFD simulation aimed at estimating muscle fiber architecture (Inouye et al., 2015), the regions of muscle origin and insertion are designated as the “inlet” and “outlet” of the fluid flow, respectively. Furthermore, the surface of the muscle is assumed to be impenetrable, except for the inlet and outlet, as demonstrated by Handsfield et al. (2017). The resultant velocity vector field produced by the CFD model can then be used to represent muscle fiber architectural arrangements. In some recent studies, the flow vectors predicted from the CFD fiber simulation approach were compared with fiber fields determined by diffusion tensor imaging (DTI), which exhibited similar fiber maps (Varvik et al., 2021). Figure 3 shows an example of CFD fiber simulation for the rectus femoris (RF) using ABAQUS/CFD (CAE Version 6.13–4). In the simulations, the origin and insertion served as inlet and outlet boundary condition surfaces of the Newtonian flow with an inlet velocity. The velocity vector map of the simulated flow was then processed to estimate the map of fiber directions.

FIGURE 3

FIGURE 3. The illustration of CFD fiber simulation: (A) flow simulation setup for the rectus femoris (RF) with inlet and outlet boundary condition surfaces, and (B) velocity vectors were exported from simulation to characterize the fiber arrangements.

2.4 Modeling of multiple muscles and anatomical interactions: quadriceps response in knee flexion-extension test

Regarding 3D active skeletal muscle simulation, most current models have been developed with a single muscle configuration (Dao and Tho, 2018). Only a limited number of studies integrated multiple active skeletal muscles into a body region (Stelletta et al., 2017). In this study, we extended the 3D active muscle model presented above into a body system level with multiple lower extremity muscles (Figure 4A), aiming to examine the effectiveness and stability of simulating motion with multiple instances of the above presented pipeline for single active muscle modeling. Simulation of the quadriceps mechanism including the knee is important for the investigation of muscle pathologies and treatments that impact lower limb function, such as acute injury (Puladi et al., 2022), fibrosis (Hung et al., 2014), and sarcoma (Houdek et al., 2021). To simulate active quadriceps, we have a baseline model of the right lower limb which included subject-specific bony structures and cartilages of the knee segmented from computed tomography (CT) and magnetic resonance imaging (MRI) (Hume et al., 2019). It has a 3 DOF ball-joint representing the hip joint, and 6 DOF joints to characterize the patellofemoral (PF) and tibiofemoral (TF) joints (Harris et al., 2016). The patellar tendon and patellofemoral ligaments were modeled as hyperelastic quadrilateral meshes embedded with nonlinear springs (Hume et al., 2019). The 3D deformable active quadriceps muscles were created using material parameters calibrated in the previously presented simplified model, including vastus intermedius VI), vastus lateralis (VL), vastus medialis (VM), and rectus femoris (RF). The 3D geometries of the quadriceps muscles were segmented from images of the Visible Human Male (Ackerman, 1998; Andreassen et al., 2023). The fiber distributions of each individual muscle were modeled using the CFD simulation method described in section 2.3. A knee flexion-extension during a squat was simulated in ABAQUS/Explicit from full extension to 90 deg knee joint flexion, then extended back to full extension. To achieve the flexion-extension motion of the knee in a squat, the pelvis was assigned a kinematic profile in the vertical direction. The distal end of tibia was pinned, and tibial rotation was permitted in the sagittal plane. At approximately 35° of knee flexion, the quadriceps muscles were activated to 50% of their maximum level (Krishnan et al., 2011; O'Bryan et al., 2022) and this level was maintained until maximum flexion was achieved.

FIGURE 4

FIGURE 4. FE model with deformable quadriceps to simulate knee flexion-extension: (A) configuration of initial position (full extension), and (B) deformed configuration with displacement map for 90 deg knee joint flexion during the squat.

It remains a computational challenge to model interactions among muscles. To simulate realistic interactions (Stecco et al., 2008), a computational solution should restrict separations between muscles under tension, but allow relative sliding caused by shearing. A feasible approach is to employ cohesive contact (CC) approach in ABAQUS, which has been used to model bonded interfaces with the possibility of damage or failure in fracture mechanics (Kulkarni et al., 2021). After implementation of the approach using several simplified models, the cohesive contact approach was then utilized to simulate the muscle interactions in the current quadriceps model. The parameters of the cohesive properties used in this study include: the stiffness coefficients K_nn = K_ss = K_tt = 0.05 N/mm³, and the damage initiation (maximum separation) D_nn = D_ss = D_tt = 30 mm. The approach can indeed transmit loads normal to the boundary of part, which restricts separation of anatomical structures by tension, but permits relative sliding. In this study, the effective stress (Zeng et al., 2021a) was reported and can be compared with previous studies (Li et al., 2019). In this study, the 95th percentile of von Mises stress (95% vM stress) based on the 95th percentile peak element response was used to characterize the model outcomes (peak stress means 95% vM stress hereafter). Similar to the idea of using 95th percentile maximum principal strain (MPS) in brain tissue biomechanics (Panzer et al., 2012), 95th percentile von Mises stress may avoid or limit potential numerical artifacts on a single element or very few elements associated with using the 100th percentile-ranked maximum element value (Gabler et al., 2019).

3 Results

In this section, some key results are provided for both validation of single muscle elongation and response of a group of muscles (quadriceps) in knee flexion-extension.

3.1 Passive and active behavior validation by elongation tests

For uniaxial tensile tests, the deformed configurations (20% elongation) of passive muscle and active muscle were compared with the experimental testing results. Differences of deformation were found between passive and active elongation simulations (Figures 2B,C). The sensitivity analysis showed that the compressibility constant D and the maximum isometric stress ${\sigma }_0$ have considerable influence on the engineering stress response, which indicates that it is critical to choose appropriate values for D and ${\sigma }_0$ (Lu et al., 2011a). The simulated responses agreed well with the corresponding experimental curves (Myers et al., 1998; Davis et al., 2003), and only slight differences were observed (e.g., the curves with more than 15% engineering strains) (Figure 5).

FIGURE 5

FIGURE 5. The engineering stress-strain response for the single muscle under uniaxial elongation: (A) comparison of passive muscle response with experimental data from Davis et al. (2003) (B) comparison of active muscle response with experimental data from Myers et al. (1998).

3.2 Response in knee flexion-extension

The deformed configuration with displacement map for 90 deg knee joint flexion was displayed in Figure 4B, which can be compared with the initial configuration with 0 deg of knee extension. Figure 6A provided the configuration of the deformable active muscles along with the rigid bones, as well as the stress distribution (up to the peak stress) in the muscles of the quadriceps under the 90 deg knee joint flexion. For the peak stress, it occurred in the ends of muscles (e.g., insertions and origins), as well as some interaction regions between anatomical structures due to contacts (e.g., the contacted areas between distal parts of VI and VM). The comparisons of the peak stress values (i.e., 95% vM stress) between passive and active muscles were collected in Figure 6B. It can be noticed that the peak stress was mostly located at the distal ends of VI and VM, likely due to the large tension in this region and associated muscle-muscle interactions and muscle-bone interactions. It is apparent that the maximum magnitude of the peak stress in each muscle produced by the model with active quadriceps was higher to a certain extent than the peak stress produced by the passive muscle model. Compared to the musculoskeletal FE model of lower extremity developed previously (Li et al., 2019) for kinematics during walking gait, the peak stresses predicted by the models in the current study were lower, but of a similar magnitude (MPa). However, the maximum stresses in Li et al. (2019) were calculated without using 95th percentile maximum element response of the stress, which can potentially be affected by numerical artifacts with using the 100th percentile value (Carlsen, et al., 2021).

FIGURE 6

FIGURE 6. The predicted effective stress in the quadriceps: (A) the stress distribution in the quadriceps muscles at the 90° knee joint flexion, and (B) comparison of the peak stress values (95% percentile of von Mises stress) between passive and active muscles for different angles (30°, 45°, 60°, and 90°) of knee extension (VL: vastus lateralis, VM: vastus medialis, RF: rectus femoris, and VI: vastus intermedius).

4 Discussion

This study presented a framework for modeling active muscles by 3D continuum approach and incorporated multiple realistic muscles with interactions into a lower limb that is capable of simulating lower extremity motion. The general computational framework of Hill-type hyperelastic 3D active muscle model was presented, which integrated passive and active properties of muscle into a continuum material model based on a strain energy approach (Blemker et al., 2005; Tang et al., 2009; Lu et al., 2011a; Khodaei et al., 2013; Kleinbach et al., 2017). Compared to discrete approaches to model the muscle using 1D beam/truss elements combined with a 3D solid matrix, the presented framework can deal with complex geometries and is more convenient to implement without the tedious process to define muscle fibers restricted by the mesh of the matrix (Feller et al., 2016; Hedenstierna et al., 2008; Mo et al., 2018 and 2019). To estimate the muscle fiber arrangements, the CFD fiber simulation approach was used and applied for muscles with 3D realistic geometry. The model of muscle was then implemented at a system level into a model of the quadriceps with multiple active muscles and detailed representation of the knee joint under knee bending conditions, which demonstrated the effectiveness and numerical stability of the presented muscle model for simulating complex musculoskeletal systems. The cohesive contact approach was adopted to model interactions between muscles and between muscle and bone. To the authors’ knowledge, this work is the first study which utilized a cohesive contact approach for modeling interactions among anatomical structures in a dynamic simulation with faithful 3D muscle geometry.

Despite the accurate representation of the skeletal muscle geometry and fiber architecture, parameter identification or determination is another challenge of modeling realistic active muscles (Dao and Tho, 2018). The material parameters of different active skeletal muscle models exhibited a wide range of values and variations, likely due to different skeletal muscle specimens used and differences in the experimental methods and modeling approaches. This study describes the active muscle model by multiple constitutive parameters, which characterize the stress in the matrix, PE, SEE, and/or CE. The uniaxial elongation simulations of both passive and active muscles using 3D ideal geometry were conducted to identify the sensitivity of main parameters (also see Lu et al., 2011a) and determine the values of the adjustable material parameters. Compared with the simulation results in Lu et al. (2011a), the presented muscle model in this study with calibrated material parameters displayed a closer correlation with the experimental curves in Figure 5, particularly for large engineering strain (e.g., >10%). Regarding the definition of fiber orientations, the physics-based approach using CFD simulation was utilized to model the subject-specific muscle fiber arrangements. The method was compared against the fiber orientations observed with DTI (Inouye et al., 2015; Handsfield et al., 2017; Varvik et al., 2021), which has shown to be a physiologically reasonable and cost-effective simulation approach for predicting fiber trajectories. Compared to the other modeling approaches to create the fiber map such as mapping technique from different fiber templates (Blemker et al., 2005), the CFD fiber simulation technique is computationally more efficient and can handle subject-specific and muscle-specific differences without template-confined limitations such as the requirement of a structured mesh.

Few simulations have been reported using active 3D skeletal muscle models, and most of the loading cases used simplified loading conditions, such as single muscle elongation, shortening, or contraction (Dao and Tho, 2018). In the context of modeling multiple active muscles behaviors, it is still comparatively rare and almost all prior methods adopted 1D discrete elements to model the active properties (Feller et al., 2016; Mo et al., 2018) with a combination of continuum elements for passive properties (Khodaei, et al., 2013). If the active muscle model is only simulated with very simple geometry under basic loading scenarios (e.g., uniaxial tension, shearing, shortening, or lengthening), the numerical stability and effectiveness of the model might be overlooked for multiple realistic muscles with large deformation within a body region or even a whole human body model. In this study, the presented method for single active muscle modeling has been extended to the quadriceps under dynamic simulation of knee bending at a system level, which involved multiple active muscles and interactions among muscles. Modeling and simulation of interactions between anatomical structures (e.g., muscles) is still challenging in several ways, including mechanical and physiological aspects, as well as numerical issues. While a small number of prior studies have included multiple 3D muscles, they were limited by either not including contact between muscles or representing contact in non-physiological ways (Stelletta et al., 2017). For example, the quadriceps were modeled without considering contacts or other interactions between muscles (Fernandez and Hunter, 2005). The common approaches to model interactions between anatomical structures include sharing nodes between the mesh of structures (e.g., the Global Human Body Models Consortium (GHBMC) extremity models (Schwartz et al., 2015; Zeng et al., 2021b), fixed constraint(s) between organs (Fernandez et al., 2005), contact (Stelletta et al., 2017), and combination of contact and springs (Böl et al., 2011). Physiologically, the muscles of the lower extremity such as quadriceps are enveloped by their deep fascias, permitting relative sliding of the muscles because of their epimysium (Stecco et al., 2008). In this work, for the first time at the quadriceps, we modeled the interactions among anatomical structures using a cohesive contact approach. It can restrict separation of the structures (e.g., muscle to muscle, muscle to bone) under tension, but still permit relative sliding due to shearing. The four bundles of the quadriceps experienced large deformation during the knee bending process (Figure 6A). It should be emphasized that the specific aim of this example was to show the numerical stability and applicability of the presented 3D active skeletal muscle modeling framework for multiple muscle configurations with anatomical interactions at a system level. Because the responses of multiple muscles in a complex system can be affected by a wide variety of factors (e.g., geometry, locations of origins and insertions, attached tendons and ligaments, surrounding structures, interactions, boundary and loading conditions, etc.), the current study did not seek a definitive quantification of the stress or strain of the muscles. However, the trends of simulated results agreed with some published studies including a previously developed musculoskeletal FE model of the lower extremity (Li et al., 2019), which showed the stresses were concentrated at the muscle’s insertion and origin regions most likely because of the small muscle cross-sectional area of the two ends of the muscles. Compared to the model presented by Li et al. (2019) for kinematics during walking gait, the peak stresses predicted by the model in the current study were lower, but of a similar magnitude (MPa). However, the maximum stresses in Li et al. (2019) were calculated without using 95th percentile maximum element response of the stress, which can potentially be affected by numerical artifacts with using the 100th percentile value (Carlsen, et al., 2021). The current study also found that large stress levels appeared at the contacted regions between muscles, including the insertion ends of the VI and VM. Additionally, as expected, the peak stress in each muscle with active properties appeared to be higher than the corresponding muscle with passive properties only (Figure 6B). The stress increased in the localized regions of the active muscle since it was stiffened under the active state due to the new contributor of active component integrated with the passive stiffness (Winters, 2012).

Although the present work has yielded promising results, there are some limitations in this study, which can be considered and improved in future research. First, the presented model was validated by the testing data of only two simple loading cases with elongation from rabbit leg muscles, which were the most used muscle type in the literature for model validation and parametric identification. There was a paucity of available experimental data regarding active muscle model validations under different external mechanical stimuli, particularly experimental data (i.e., electrophysiological recordings) of human active skeletal muscles in pathophysiological conditions and for specific aspects of the numerical implementation of the model. Likewise, there is no model validation data developed or presented in the literature for the purposes of the active quadriceps response in knee bending. Additionally, it is possible that the knee angle-dependent differences could be influenced by variations in muscle activation levels (Kooistra et al., 2005). The current study employed the same activation pattern for all muscles, and did not account for the gradient recruitment of fibers along the muscle due to the presence of different neuromotor units within the muscle. Second, there were some model simplifications due to factors such as insufficient representation of anatomical details or material properties. For example, the patellar tendon and patellofemoral ligaments were not calibrated, although the knee joint in the adopted baseline model was validated by comparing kinematics of the FE model with experimental kinematic data (Baldwin et al., 2009; Ali et al., 2016). Third, regarding the characterization of multiple muscle interaction, the quantification of force transmission between muscles remains a challenge from both experimental and simulation perspectives. Although we believe that the presented cohesive contact approach is more efficient and realistic than other existing approaches (Fernandez and Hunter, 2005; Böl et al., 2011; Schwartz et al., 2015; Zeng et al., 2021b) to model organ interactions based on the performance compared to the limited available evidence, the interaction effects on the muscle forces have not been quantified since there is a lack of experimental testing data for validation thus far. Finally, the external fascia and skin were not considered in the current quadriceps model, so further research is needed to accurately represent the thigh to consider the muscle packing effect (Stecco et al., 2008).

5 Conclusion

This study proposed a computational framework of 3D muscle simulation by Hill-type hyperelastic active muscle model, which was validated by single muscle loading cases and extended for modeling multiple muscles with interactions during movement at a system level. The presented constitutive model and relevant material parameters were validated by the experimental data based on the rabbit TA muscle (Myers et al., 1998; Davis et al., 2003). The simulation results were in good agreement with the experimental measured elongation responses, which indicated that the current model with determined parameters can characterize the passive and active behaviors of muscles. After conducting single muscle level validation, the active muscle model was implemented into multiple muscle configurations, i.e., the quadriceps response in knee flexion-extension. For each muscle with complex shape, the muscle fiber orientations were defined using the CFD fiber simulation approach. For the interactions among anatomical structures, the effect of contact between muscles and between muscle and femur was qualitatively simulated by the cohesive contact approach, which can restrict separation of the components under tension, but still permit relative sliding due to shearing. The muscle model was stable and effective for modeling multiple realistic 3D muscles with interactions, and capable of predicting the trends of muscle responses through dynamics simulation of knee flexion and extension. Therefore, it is feasible for the presented framework to link the tissue-level active muscle behavior and dynamics of human movements. It is expected that this study can encourage additional computational and experimental studies of skeletal muscle mechanics, particularly muscle activation, quantification of muscle interactions, constitutive model development, parameters determination, and multiscale modeling for neuromuscular dynamic simulation. In a larger scope, this study provides a promising framework toward developing realistic multiscale active human body models for orthopedic treatment, rehabilitation, implant mechanics, and impact biomechanics.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This study was supported in part by the NIH U01 AR072989 project, funded by the NIH National Institute of Arthritis and Musculoskeletal and Skin Diseases, the National Institute of Biomedical Imaging and Bioengineering, and the National Institute of Child Health and Human Development.

Conflict of interest

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

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2023.1153692/full#supplementary-material

References

Ackerman, M. J. (1998). The visible human project. Proc. IEEE 86, 134–140. doi:10.3233/SHTI210988

CrossRef Full Text | Google Scholar

Ali, A. A., Shalhoub, S. S., Cyr, A. J., Fitzpatrick, C. K., Maletsky, L. P., Rullkoetter, P. J., et al. (2016). Validation of predicted patellofemoral mechanics in a finite element model of the healthy and cruciate-deficient knee. J. Biomech. 49 (2), 302–309. doi:10.1016/j.jbiomech.2015.12.020

PubMed Abstract | CrossRef Full Text | Google Scholar

Andreassen, T. E., Hume, D. R., Hamilton, L. D., Walker, K. E., Higinbotham, S. E., and Shelburne, K. B. (2023). Three dimensional lower extremity musculoskeletal geometry of the visible human female and male. Sci. Data 10 (1), 34. doi:10.1038/s41597-022-01905-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Baldwin, M. A., Clary, C., Maletsky, L. P., and Rullkoetter, P. J. (2009). Verification of predicted specimen-specific natural and implanted patellofemoral kinematics during simulated deep knee bend. J. Biomech. 42 (14), 2341–2348. doi:10.1016/j.jbiomech.2009.06.028

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 (5), 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 (4), 657–665. doi:10.1016/j.jbiomech.2004.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Blemker, S. S. (2017). “Chapter 17 Three-dimensional modeling of active muscle tissue: The why, the how, and the future,” in Biomechanics of living organs. Editors Y. Payan, and J. Ohayon (Oxford: Academic Press), 361–375.

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 (5), 489–504. doi:10.1080/10255840701771750

PubMed Abstract | CrossRef Full Text | Google Scholar

Böl, M., Sturmat, M., Weichert, C., and Kober, C. (2011). A new approach for the validation of skeletal muscle modelling using MRI data. Comput. Mech. 47, 591–601. doi:10.1007/s00466-010-0567-0

CrossRef Full Text | Google Scholar

Carlsen, R. W., Fawzi, A. L., Wan, Y., Kesari, H., and Franck, C. (2021). A quantitative relationship between rotational head kinematics and brain tissue strain from a 2-D parametric finite element analysis. Brain Multiphys 2, 100024. doi:10.1016/j.brain.2021.100024

CrossRef Full Text | Google Scholar

Chen, D. T., and Zeltzer, D. (1992). Pump it up: Computer animation of a biomechanically based model of muscle using the finite element method. Comput. Graph. 26, 89–98. doi:10.1145/142920.134016

CrossRef Full Text | Google Scholar

Clemen, C. B., Benderoth, G. E. K., Schmidt, A., Hübner, F., Vogl, T. J., and Silber, G. (2017). Human skeletal muscle behavior in vivo: Finite element implementation, experiment, and passive mechanical characterization. J. Mech. Behav. Biomed. Mat. 65, 679–687. doi:10.1016/j.jmbbm.2016.09.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Close, R. (1964). Dynamic properties of fast and slow skeletal muscles of the rat during development. J. Physiol. 173, 74–95. doi:10.1113/jphysiol.1964.sp007444

PubMed Abstract | CrossRef Full Text | Google Scholar

Dao, T. T., and Tho, M. H. B. (2018). A systematic review of continuum modeling of skeletal muscles: Current trends, limitations, and recommendations. Appl. Bionics Biomech. 2018, 1–17. doi:10.1155/2018/7631818

CrossRef Full Text | Google Scholar

Davis, J., Kaufman, K. R., and Lieber, R. L. (2003). Correlation between active and passive isometric force and intramuscular pressure in the isolated rabbit tibialis anterior muscle. J. Biomech. 36 (4), 505–512. doi:10.1016/s0021-9290(02)00430-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Elyasi, E., Perrier, A., Bailet, M., and Payan, Y. (2022). Biomechanical lower limb model to predict patellar position alteration after medial open wedge high tibial osteotomy. J. Biomech. 136, 111062. doi:10.1016/j.jbiomech.2022.111062

PubMed Abstract | CrossRef Full Text | Google Scholar

Feller, L., Kleinbach, C., Fehr, J., and Schmitt, S. (2016). “Incorporating muscle activation dynamics into the global human body model,” in International IRCOBI Conference, Malaga, Spain, September 14–16, 2016.

Google Scholar

Fernandez, J. W., and Hunter, P. J. (2005). An anatomically based patient-specific finite element model of patella articulation: Towards a diagnostic tool. Biomech. Model. Mechanobiol. 4 (1), 20–38. doi:10.1007/s10237-005-0072-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Fernandez, J. W., Buist, M. L., Nickerson, M. L., and Hunter, P. J. (2005). Modelling the passive and nerve activated response of the rectus femoris muscle to a flexion loading: A finite element framework. Med. Eng. Phys. 27 (10), 862–870. doi:10.1016/j.medengphy.2005.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Fung, Y. C. (1993). Biomechanics: Mechanical properties of living tissues. 2nd ed. New York, NY: Springer-Verlag.

Google Scholar

Gabler, L. F., Crandall, J. R., and Panzer, M. B. (2019). Development of a second-order system for rapid estimation of maximum brain strain. Ann. Biomed. Eng. 47 (9), 1971–1981. doi:10.1007/s10439-018-02179-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Göktepe, S., Menzel, A., and Kuhl, A. (2014). The generalized Hill model: A kinematic approach towards active muscle contraction. J. Mech. Phys. Solids 72, 20–39. doi:10.1016/j.jmps.2014.07.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, A. M., Huxley, A. F., and Julian, F. J. (1966). The variation in isometric tension with sarcomere length in vertebrate muscle fibres. J. Physiol. 184 (1), 170–192. doi:10.1113/jphysiol.1966.sp007909

PubMed Abstract | CrossRef Full Text | Google Scholar

Gotti, C., Sensini, A., Zucchelli, A., Carloni, R., and Focarete, M. L. (2020). Hierarchical fibrous structures for muscle-inspired soft-actuators: A review. Appl. Mat. Today 20, 100772. doi:10.1016/j.apmt.2020.100772

CrossRef Full Text | Google Scholar

Grasa, J., Ramírez, A., Osta, R., Muñoz, M. J., Soteras, F., and Calvo, B. (2011). A 3D active-passive numerical skeletal muscle model incorporating initial tissue strains. Validation with experimental results on rat tibialis anterior muscle. Biomech. Model. Mechanobiol. 10 (5), 779–787. doi:10.1007/s10237-010-0273-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Handsfield, G. G., Bolsterlee, B., Inouye, J. M., Herbert, R. D., Besier, T. F., and Fernandez, J. W. (2017). Determining skeletal muscle architecture with laplacian simulations: A comparison with diffusion tensor imaging. Biomech. Model. Mechanobiol. 16 (6), 1845–1855. doi:10.1007/s10237-017-0923-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, M. D., Cyr, A. J., Ali, A. A., Fitzpatrick, C. K., Rullkoetter, P. J., Maletsky, L. P., et al. (2016). A combined experimental and computational approach to subject-specific analysis of knee joint laxity. J. Biomech. Eng. 138 (8), 0810041–0810048. doi:10.1115/1.4033882

PubMed Abstract | CrossRef Full Text | Google Scholar

Hedenstierna, S., Halldin, P., and Brolin, K. (2008). Evaluation of a combination of continuum and truss finite elements in a model of passive and active muscle tissue. Comput. Methods Biomech. Biomed. Engin. 11 (6), 627–639. doi:10.1080/17474230802312516

PubMed Abstract | CrossRef Full Text | Google Scholar

Heidlauf, T., Klotz, T., Rode, C., Altan, E., Bleiler, C., Siebert, T., et al. (2016). A multi-scale continuum model of skeletal muscle mechanics predicting force enhancement based on actin-titin interaction. Biomech. Model. Mechanobiol. 15 (6), 1423–1437. doi:10.1007/s10237-016-0772-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hill, A. V. (1938). The heat of shortening and dynamics constants of muscles. Proc. R. Soc. B 126 (843), 136–195. doi:10.1098/rspb.1938.0050

CrossRef Full Text | Google Scholar

Hill, A. V. (1970). First and last experiments in muscle mechanics. Cambridge, UK: Cambridge University Press.

Google Scholar

Houdek, M. T., Wellings, E. P., Mallett, K. E., Honig, R. L., Rose, P. S., and Moran, S. L. (2021). Free functional latissimus dorsi reconstruction of the quadriceps and hamstrings following oncologic resection of soft tissue sarcomas of the thigh. Sarcoma 2021, 1–7. doi:10.1155/2021/8480737

CrossRef Full Text | Google Scholar

Hume, D. R., Navacchia, A., Ali, A. A., and Shelburne, K. B. (2018). The interaction of muscle moment arm, knee laxity, and torque in a multi-scale musculoskeletal model of the lower limb. J. Biomech. 76, 173–180. doi:10.1016/j.jbiomech.2018.05.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Hume, D. R., Navacchia, A., Rullkoetter, P. J., and Shelburne, K. B. (2019). A lower extremity model for muscle-driven simulation of activity using explicit finite element modeling. J. Biomech. 84, 153–160. doi:10.1016/j.jbiomech.2018.12.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, J. D., and Yin, F. C. (1987). On constitutive relations and finite deformations of passive cardiac tissue: I. A pseudostrain-energy function. J. Biomech. Eng. 109 (4), 298–304. doi:10.1115/1.3138684

PubMed Abstract | CrossRef Full Text | Google Scholar

Hung, N. N., Tan, D., and Do Ngoc Hien, N. (2014). Patellar dislocation due to iatrogenic quadriceps fibrosis: Results of operative treatment in 54 cases. J. Child. Orthop. 8 (1), 49–59. doi:10.1007/s11832-014-0564-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Huxley, A. F. (1957). Muscle structure and theories of contraction. Prog. Biophys. Biophys. Chem. 7, 255–318. doi:10.1016/s0096-4174(18)30128-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Inouye, J. M., Handsfield, G. G., and Blemker, S. S. (2015). “Fiber tractography for finite-element modeling of transversely isotropic biological tissues of arbitrary shape using computational fluid dynamics,” in Proceedings of the Conference on Summer Computer Simulation, Chicago, IL, July 26 - 29, 2015 (): Society for Computer Simulation International).

Google Scholar

Khodaei, H., Mostofizadeh, S., Brolin, K., Johansson, H., and Osth, J. (2013). Simulation of active skeletal muscle tissue with a transversely isotropic viscohyperelastic continuum material model. Proc. Inst. Mech. Eng. H. 227 (5), 571–580. doi:10.1177/0954411913476640

PubMed Abstract | CrossRef Full Text | Google Scholar

Kleinbach, C., Martynenko, O., Promies, J., Haeufle, D. F. B., Fehr, J., and Schmitt, S. (2017). Implementation and validation of the extended Hill-type muscle model with robust routing capabilities in LS-DYNA for active human body models. Biomed. Eng. Online 16 (1), 109. doi:10.1186/s12938-017-0399-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Klotz, T., Bleiler, C., and Röhrle, O. (2021). A physiology-guided classification of active-stress and active-strain approaches for continuum-mechanical modeling of skeletal muscle tissue. Front. Physiol. 12, 685531. doi:10.3389/fphys.2021.685531

PubMed Abstract | CrossRef Full Text | Google Scholar

Kooistra, R. D., de Ruiter, C. J., and de Haan, A. (2005). Muscle activation and blood flow do not explain the muscle length-dependent variation in quadriceps isometric endurance. J. Appl. Physiol. (1985) 98 (3), 810–816. doi:10.1152/japplphysiol.00712.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Krishnan, C., Allen, E. J., and Williams, G. N. (2011). Effect of knee position on quadriceps muscle force steadiness and activation strategies. Muscle Nerve 43 (4), 563–573. doi:10.1002/mus.21981

PubMed Abstract | CrossRef Full Text | Google Scholar

Kulkarni, S. S., Gupta, V., Ortiz, A., Das, H., Upadhyay, P., Barker, E., et al. (2021). Determining cohesive parameters for modeling interfacial fracture in dissimilar-metal friction stir welded joints. Int. J. Solids Struct. 216, 200–210. doi:10.1016/j.ijsolstr.2021.01.023

CrossRef Full Text | Google Scholar

Li, J., Lu, Y., Miller, S. C., Jin, Z., and Hua, X. (2019). Development of a finite element musculoskeletal model with the ability to predict contractions of three-dimensional muscles. J. Biomech. 94, 230–234. doi:10.1016/j.jbiomech.2019.07.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y. T., Zhu, H. X., Richmond, S., and Middleton, J. (2010). A visco-hyperelastic model for skeletal muscle tissue under high strain rates. J. Biomech. 43 (13), 2629–2632. doi:10.1016/j.jbiomech.2010.05.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y. T., Beldie, L., Walker, B., Richmond, S., and Middleton, J. (2011a). Parametric study of a Hill-type hyperelastic skeletal muscle model. Proc. Inst. Mech. Eng. H. 225 (5), 437–447. doi:10.1177/2041303310392632

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y. T., Zhu, H. X., Richmond, S., and Middleton, J. (2011b). Modelling skeletal muscle fibre orientation arrangement. Comput. Methods Biomech. Biomed. Engin. 14 (12), 1079–1088. doi:10.1080/10255842.2010.509100

PubMed Abstract | CrossRef Full Text | Google Scholar

Martins, J. A. C., Pires, E. B., Salvado, R., and Dinis, P. B. (1998). A numerical model of passive and active behavior of skeletal muscles. Comput. Methods Appl. Mech. Eng. 151 (3–4), 419–433. doi:10.1016/S0045-7825(97)00162-X

CrossRef Full Text | Google Scholar

Martins, J. A. C., Pato, M. P. M., and Pires, E. B. (2006). A finite element model of skeletal muscles. Virtual Phys. Prototyp. 1 (3), 159–170. doi:10.1080/17452750601040626

CrossRef Full Text | Google Scholar

Meier, P., and Blickhan, R. (2000). “FEM-Simulation of skeletal muscle: The influence of inertia during activation and deactivation,” in Skeletal muscle mechanics: From mechanisms to function. Editor W. Herzog (Chichester: John Wiley and Sons), 207–224.

Google Scholar

Mo, F., Li, F., Behr, M., Xiao, Z., Zhang, G., and Du, X. (2018). A lower limb-pelvis finite element model with 3D active muscles. Ann. Biomed. Eng. 46 (1), 86–96. doi:10.1007/s10439-017-1942-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Mo, F., Li, J., Dan, J., Liu, T., and Behr, T. (2019). Implementation of controlling strategy in a biomechanical lower limb model with active muscles for coupling multibody dynamics and finite element analysis. J. Biomech. 91, 51–60. doi:10.1016/j.jbiomech.2019.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Myers, B. S., Woolley, C. T., Slotter, T. L., Garrett, W. E., and Best, T. M. (1998). The influence of strain rate on the passive and stimulated engineering stress--large strain behavior of the rabbit tibialis anterior muscle. J. Biomech. Eng. 120 (1), 126–132. doi:10.1115/1.2834292

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Bryan, S. J., Taylor, J. L., D'Amico, J. M., and Rouffet, D. M. (2022). Quadriceps muscle fatigue reduces extension and flexion power during maximal cycling. Front. Sports Act. Living 3, 797288. doi:10.3389/fspor.2021.797288

PubMed Abstract | CrossRef Full Text | Google Scholar

Panzer, M. B., Myers, B. S., Capehart, B. P., and Bass, C. R. (2012). Development of a finite element model for blast brain injury and the effects of CSF cavitation. Ann. Biomed. Eng. 40 (7), 1530–1544. doi:10.1007/s10439-012-0519-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinto, J. G., and Fung, Y. C. (1973). Mechanical properties of the heart muscle in the passive state. J. Biomech. 6 (6), 597–616. doi:10.1016/0021-9290(73)90017-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Puladi, B., Ooms, M., Geijtenbeek, T., Trinler, U., Houschyar, K. S., Gruber, L. J., et al. (2022). Tolerable degree of muscle sacrifice when harvesting a vastus lateralis or myocutaneous anterolateral thigh flap. J. Plast. Reconstr. Aesthet. Surg. 77, 94–103. doi:10.1016/j.bjps.2022.10.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Röhrle, O., and Pullan, A. J. (2007). Three-dimensional finite element modelling of muscle forces during mastication. J. Biomech. 40 (15), 3363–3372. doi:10.1016/j.jbiomech.2007.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Röhrle, O., Davidson, J. B., and Pullan, A. J. (2012). A physiologically based, multi-scale model of skeletal muscle structure and function. Front. Physiol. 3, 358. doi:10.3389/fphys.2012.00358

PubMed Abstract | CrossRef Full Text | Google Scholar

Roux, A., Lecompte, J., Iordanoff, I., and Laporte, S. (2021). Modeling of muscular activation of the muscle-tendon complex using discrete element method. Comput. Methods Biomech. Biomed. Engin. 24 (11), 1184–1194. doi:10.1080/10255842.2020.1870039

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwartz, D., Guleyupoglu, B., Koya, B., Stitzel, J. D., and Gayzik, F. S. (2015). Development of a computationally efficient full human body finite element model. Traffic Inj. Prev. 16 (1), S49–S56. doi:10.1080/15389588.2015.1021418

PubMed Abstract | CrossRef Full Text | Google Scholar

Spyrou, L. A., and Aravas, N. (2012). Muscle-driven finite element simulation of human foot movements. Comput. Methods Biomech. Biomed. Engin. 15 (9), 925–934. doi:10.1080/10255842.2011.566564

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. Mat. 92, 97–117. doi:10.1016/j.jmbbm.2018.12.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Stecco, C., Porzionato, A., Lancerotto, L., Stecco, A., Macchi, V., Day, J. A., et al. (2008). Histological study of the deep fasciae of the limbs. J. Bodyw. Mov. Ther. 12 (3), 225–230. doi:10.1016/j.jbmt.2008.04.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Stelletta, J., Dumas, R., and Lafon, Y. (2017). “Chapter 23- modeling of the thigh: A 3D deformable approach considering muscle interactions,” in Biomechanics of living organs. Editors Y. Payan, and J. Ohayon (Oxford: Academic Press), 497–521.

CrossRef Full Text | Google Scholar

Tang, C. Y., Tsui, C. P., Stojanovic, B., and Kojic, M. (2007). Finite element modelling of skeletal muscles coupled with fatigue. Int. J. Mech. Sci. 49 (10), 1179–1191. doi:10.1016/j.ijmecsci.2007.02.002

CrossRef Full Text | Google Scholar

Tang, C. Y., Zhang, G., and Tsui, G. (2009). A 3D skeletal muscle model coupled with active contraction of muscle fibres and hyperelastic behaviour. J. Biomech. 42 (7), 865–872. doi:10.1016/j.jbiomech.2009.01.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Ting, L. H., Chvatal, S. A., Safavynia, S. A., and McKay, J. L. (2012). Review and perspective: Neuromechanical considerations for predicting muscle activation patterns for movement. Int. J. Numer. Method Biomed. Eng. 28 (10), 1003–1014. doi:10.1002/cnm.2485

PubMed Abstract | CrossRef Full Text | Google Scholar

Varvik, J., Besier, T. F., and Handsfield, G. G. (2021). Computational fluid dynamics simulations for 3D muscle fiber architecture in finite element analysis: Comparisons between computational fluid dynamics and diffusion tensor imaging. Int. J. Numer. Method Biomed. Eng. 37 (12), e3521. doi:10.1002/cnm.3521

PubMed Abstract | 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 (8), 0810061–0810068. doi:10.1115/1.4040318

PubMed Abstract | CrossRef Full Text | Google Scholar

Winters, T. M. (2012). Determinants of active and passive tension in skeletal muscle. PhD Dissertation. San Diego, CA: University of California.

Google Scholar

Wu, T., Hung, T., and Mithraratne, K. (2014). Generating facial expressions using an anatomically accurate biomechanical model. IEEE Trans. Vis. Comput. Graph 20 (11), 1519–1529. doi:10.1109/TVCG.2014.2339835

PubMed Abstract | CrossRef Full Text | Google Scholar

Zajac, F. (1989). Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control. Crit. Rev. Biomed. Eng. 17, 359–411.

PubMed Abstract | Google Scholar

Zeng, W., Lewicki, K. A., Chen, Z., and Van Citters, D. W. (2021a). The evaluation of reverse shoulder lateralization on deltoid forces and scapular fracture risk: A computational study. Med. Nov. Technol. Devices 11, 100076. doi:10.1016/j.medntd.2021.100076

CrossRef Full Text | Google Scholar

Zeng, W., Mukherjee, S., Caudillo, A., Forman, J., and Panzer, M. B. (2021b). Evaluation and validation of thorax model responses: A hierarchical approach to achieve high biofidelity for thoracic musculoskeletal system. Front. Bioeng. Biotechnol. 9, 712656. doi:10.3389/fbioe.2021.712656

PubMed Abstract | CrossRef Full Text | Google Scholar