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.

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.

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).

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.