How muscle stiffness affects human body model behavior

Background Active human body models (AHBM) consider musculoskeletal movement and joint stiffness via active muscle truss elements in the finite element (FE) codes in dynamic application. In the latest models, such as THUMS™ Version 5, nearly all human muscle groups are modeled in form of one-dimensional truss elements connecting each joint. While a lot of work has been done to improve the active and passive behavior of this 1D muscle system in the past, the volumetric muscle system of THUMS was modeled in a much more simplified way based on Post Mortem Human Subject (PMHS) test data. The stiffness changing effect of isometric contraction was hardly considered for the volumetric muscle system of whole human body models so far. While previous works considered this aspect for single muscles, the effect of a change in stiffness due to isometric contraction of volumetric muscles on the AHBM behavior and computation time is yet unknown. Methods In this study, a simplified frontal impact using the THUMS Version 5 AM50 occupant model was simulated. Key parameters to regulate muscle tissue stiffness of solid elements in THUMS were identified for the material model MAT_SIMPLIFIED_FOAM and different stiffness states were predefined for the buttock and thigh. Results During frontal crash, changes in muscle stiffness had an effect on the overall AHBM behavior including expected injury outcome. Changes in muscle stiffness for the thigh and pelvis, as well as for the entire human body model and for strain-rate-dependent stiffness definitions based on literature data had no significant effect on the computation time. Discussion Kinematics, peak impact force and stiffness changes were in general compliance with the literature data. However, different experimental setups had to be considered for comparison, as this topic has not been fully investigated experimentally in automotive applications in the past. Therefore, this study has limitations regarding validation of the frontal impact results. Conclusion Variations of default THUMS material model parameters allow an efficient change in stiffness of volumetric muscles for whole AHBM applications. The computation time is unaffected by altering muscle stiffness using the method suggested in this work. Due to a lack of validation data, the results of this work can only be validated with certain limitations. In future works, the default material models of THUMS could be replaced with recently published models to achieve a possibly more biofidelic muscle behavior, which would even allow a functional dependency of the 1D and 3D muscle systems. However, the effect on calculation time and model stability of these models is yet unknown and should be considered in future studies for efficient AHBM applications.

1. The effect of isometric contraction of lower extremity muscles on occupant safety was analyzed by defining different stiffness states for the volumetric muscle elements of THUMS Version 5. Altering the muscle stiffness influenced the overall active human body model (AHBM) behavior during frontal crashes, as well as the predicted injury outcome. 2. Altering the muscle material stiffness via the scaling factor of the ordinate value (SFO) of the load curve referenced in all muscle and soft tissues of THUMS in the LS-DYNA material model MAT_SIMPLIFIED_FOAM has no significant impact on the computation time. 3. Stiffness parameters of the human body model THUMS ™ Version 5 are identified, varied and compared to experimental literature data for verification. Results are in general correlation, but verification by precise numerical simulations of the experimental setups should be followed in future studies. 4. Based on strain-dependent injury prediction, the injury risk for muscle and soft tissues during frontal impact is reduced with higher muscle stiffness. Further, based on contact force evaluation, increasing muscle stiffness leads to higher probability of hip fracture or dislocation and to lower probability of knee injuries. Higher risk with increasing muscle stiffness was also found for effective plastic strain evaluation, while first principal strain results rather show an arbitrary influence of muscle stiffness on the injury risk of cortical and spongy bones. Further studies are necessary to investigate the exact influence of muscle stiffness on bone injury risk in detail.

Background
This work shall provide major explications and additional information of a study previously presented by the authors [1].

Development of active human body models
In the past decades, several changes have been implemented to bring finite element models for the mechanical behavior of human bodies closer to biological reality. These human body models (HBM) now contain, e.g., facial bones, internal organs and a detailed spine model. The latest features are contractible, one-dimensional Hill-type muscle elements, allowing the HBM to perform musculoskeletal movements comparable to a living human being. Based on this implementation, the term  20:53 'active human body model' (AHBM) was established. Hill-type muscle elements follow an activation curve, which generates forces at the joints, thus enabling musculoskeletal movement and increasing the joint stiffness during simulation [2,3]. The use of seatbelt elements to represent tendons and different combinations of sliprings allows the definition of complex paths, where each muscle is adjustable separately by individual activation curves [3]. This mechanism was implemented in the latest Total HUman Model for Safety (THUMS ™ ) Version 5 1 [2] and to a certain extent in other HBMs, e.g., for neck muscles in the open source VIVA model [4] or the models of the GHBMC family [5]. THUMS is one of the few models in which all major skeletal muscles (except facial muscles) were incorporated based on the 1D Hill-type muscle system (Fig. 1a) [6].

Differences between muscle models and real muscles
Despite this complex modeling approach, certain features of biological muscles were hardly taken into consideration in the latest AHBM. One reason for this might be that biological muscles are too complex to implement every detail in the corresponding FE model, especially for multiple muscle systems; or that some of these might not be relevant for the analysis of mechanical muscle behavior [7] or AHBM applications. The features include anatomical, physiological, and material properties of muscles. Anatomical differences comprise the number and shape of the volumetric muscles. Approaches that address the issue in the context of HBM to a certain degree are the detailed GHBMC [8] and the Active THUMS [9]. Physiological differences include all biochemical and electrophysiological pathways and molecular interaction that provide the basis for muscle contraction and relaxation. To a certain extent, this issue was addressed in an 'extended Hill-type muscle model' for multiple Hill-type muscle models applicable for AHBM, considering Ca 2+ levels [10] and a material model from the LS-DYNA library (MAT_ANISO-TROPIC_HYPERELASTIC), considering, e.g., calcium concentration for the determination of active stresses of volumetric muscle elements [11]. Other models even considered single actin-myosin interactions on the molecular level [7,12,13], but due to their high complexity, these models are less suitable for multiple muscle [7] or HBM applications. Differences regarding mechanical properties are the main topic of this work. Based on the modeling approaches in AHBM, we need to distinguish between two systems: On the one hand, a 1D system of Hill-type muscle truss elements, where passive and active muscle behavior from volunteer tests were implemented and optimized over the years [2,6,10]. The 1D system (Fig. 1b-red truss elements) can generate forces at the joints using predefined activation curves, thus enabling musculoskeletal movement and adjusting the joint stiffness during simulation. On the other hand, a 3D system of volumetric (solid) muscle elements, where the Page 4 of 39 Trube et al. BioMed Eng OnLine (2021) 20:53 underlying material data were obtained from characterization experiments with Post Mortem Human Subjects (PMHS), as it is the case for the THUMS model [2,14]. The 3D system depicts the volumetric shape of muscles, where multiple tissues were summarized mainly into two simplified parts (muscle and soft tissue) with a layered build-up, rather than representing an anatomically correct depiction of the shape of real individual muscles. For example, in case of the thigh muscles, two volumetric barrel-shaped parts surround the femur, representing muscle and soft tissues, which were modeled using identical material parameters (Fig. 1b, inner and outer layer of transparent skin-colored tissues). This layered build-up was used for nearly all volumetric muscle and soft tissues in THUMS. This work will only focus on the 3D system of volumetric muscle elements. The data implemented in the corresponding LS-DYNA material models (e.g., MAT_ SIMPLIFIED_FOAM) of volumetric muscle elements rather represent a cadaver than a living human, based on the validation data from PMHS tests. If a braced condition of an occupant is modeled, the stiffness of the 3D tissue (MAT_SIMPLIFIED_FOAM) should be higher than the cadaver tissue stiffness due to isometric contraction [15] to represent a living human realistically. The question is simply how much the stiffness should change between an entirely relaxed, a partly tensed and a tetanically contracted muscle. To answer this question, several studies on single muscles were conducted in the past and yet, no complete consensus was found. Different measurement techniques and muscle samples led to a vast amount of varying results, e.g., for the Young's moduli of several orders of magnitude [16], as shown in Table 1 and Fig. 21. A consensus in the field of human muscle material characterization regarding measurement techniques, measurement parameters and muscle samples would be necessary to further exploit the capabilities of HBM. This would allow further development of multiple muscle models in the field of occupant safety by providing in vivo material data. The use of model populations to determine mechanical muscle properties for FE simulation can be a helpful tool in the process [17].

Modeling approaches addressing volumetric muscle stiffness changes
Progress regarding contractible 3D muscles was recently made with the first forward dynamic contractible muscle system [18] of a biceps-triceps system. Other works rely on a coupled system, where contractible 1D truss elements are interwoven in 3D solid elements, resulting in a shape and stiffness change of the solid elements due to 1D contraction [9,[19][20][21]. This includes the development of the Active THUMS [9,20]. Another study focused on the deformation processes under varying contraction states of a biceps [22]. These models take muscle stiffness changes into account based on different approaches. In the forward dynamic approach, a shape change resulting from the 3D contraction simulation leads to an increase in stiffness. The Active THUMS contains multiple 3D skeletal muscle reconstructed from MRI data that can change their shape due to contractible, interwoven 1D Hill-type muscle elements that share nodes with the 3D solid muscle elements. In consequence, this also leads to a change in muscle stiffness. This approach was applied to braced steering motion [20] and an entire HBM in pedestrian and occupant load cases [9,21]. Another study connected a 3D muscle model of the lower extremities with the upper body of a Hybrid III dummy model for investigation of possible effects on occupant kinematics and acting forces [23]. However, these approaches do either not consider muscle stiffness on the material card level or might currently not be applicable to entire human body models in an efficient way with reasonable computation time.
The computation time is a crucial aspect for any finite element model and simulation, because it directly reflects the efficiency of the modeling approach. For any biomechanical model, a trade-off needs to be found between biological accuracy on the one hand and calculation speed and data evaluability on the other [7]. Some of these are assumed to be neglectable for mechanical problems in the field of impact biomechanics, such as physiological differences that were described earlier. The THUMS V5 was modeled in an efficient way to allow research and industrial applications. By doing so, certain simplifications were accepted by the developers [2], such as barrel-shaped volumetric muscle and soft tissues (Figs. 1b,19) instead of anatomical muscle paths, muscle-tendon units and individually distributed fat tissues. Further, a hyperelastic material model (MAT_SIMPLIFIED_FOAM), described in detail in "Material model", was used to model muscle and soft tissues of THUMS in a simplified way in LS-DYNA.
Other modeling approaches that showed certain advantages over the default THUMS material model were previously presented. On the one hand, the material model MAT_TIS-SUE_DISPERSED was used [24], which allows the simultaneous input of mechanical passive and active muscle tissue behavior, including the option to use activation curves as input to scale 3D muscle stiffness over time. With this approach, the authors could adequately model the muscle behavior obtained from experiment with a living rabbit [25]. On the other hand, MAT_ANISOTROPIC_HYPERELASTIC is one of the latest additions to the default LS-DYNA library, designed to represent biological soft tissues, specifically muscle tissues. The complexity of the material model can be easily altered according to the desired application and tissue type. Further, additional features, such as activation curve dependent stiffness scaling can be easily implemented [11]. The effect on the computation time of these modeling approaches, when implemented in human body models, is yet unknown.

Aim of this work
The focus of this work addresses the question whether isometric contraction and the resulting change in muscle stiffness can be considered in a biofidelic way using the default material models of THUMS V5 (MAT_SIMPLIFIED_FOAM) without increasing the calculation time in view of efficient research and industrial applications of the altered THUMS model. As one of the most likely application fields of this AHBM adaption is the automotive sector, a simplified frontal impact simulation was chosen to analyze possible effects on occupant safety caused by muscle stiffness changes. Based on a previously presented idea by the authors [1], the material stiffness of the 3D muscle tissue will be predefined. A material stiffness parameter is identified and altered for the thigh and the pelvis to define four different stiffness states. Further, to analyze the effect of a 'worst case scenario' on the computation time, the stiffness of all muscle and soft tissue of THUMS was modified for an exemplary simulation, not only for the pelvis and thigh. Additionally, a fifth state was defined, considering the highly strain-rate dependent passive properties of muscles ( Table 1, Fig. 21), by incorporating data from Myers et al. [25] in all MAT_SIMPLI-FIED_FOAM parts of THUMS by defining a Table ID to compare computation time and material properties.
The validatability of the simulation is limited by two factors. On the one hand, ethical reasons, as an experimental representation of the frontal impact simulations of this work including volunteers is not feasible and was not done to this extent in the past. On the other hand, literature data of PMHS tests cannot be used for validation due to physiological reasons, as isometric muscle contraction and the resulting, voluntary muscle stiffness change is no longer possible. Therefore, related literature [16,[26][27][28][29], where the kinematic and deformation behavior in loading experiments with volunteers was analyzed, was used for validation of the results of the frontal impact simulation. The effects of these stiffness changes on the tissue motion, peak impact force, computation time and the expected injury risk for the buttock, thigh and hip region will be analyzed in detail. As both, muscle and soft tissues, are modeled in a simplified way in THUMS (Fig. 19), both can be assumed to represent muscle tissue. Therefore, the stiffness of muscle and soft tissue of THUMS were scaled in the exact same manner.
The buttock, thigh and hip region (Fig. 19) were chosen for detailed investigation in this work, as they are relevant for lower extremity injuries occurring during frontal impacts due to the direct contact to the driver's seat [30][31][32][33]. Most of the lower extremity injuries that result from automotive crashes occur in the knee, thigh and hip region [33][34][35][36][37][38][39]. According to data from 1993 to 1997 of the National Automotive Sampling System (NASS), 21.5% of all occupant injuries were accounted to lower extremity injuries [30]. Accordingly, they are the second most frequent injuries considered in the NASS database. Although lower extremity injuries are in most cases not life threating, they might result in long term hospitalization and physical disabilities [36][37][38][39]. Besides physical limitations, lower extremity injuries can have psychosocial long term effects, such as depression [38]. Both can cause high associated costs for the social sector. Further, this study was motivated by the hypothesis that one of the main function of muscles is the dissipation of mechanical energy after external impacts [40]. Based on these circumstances, further investigation in the field of lower extremity injuries and force effects on this body region during an automotive impact are necessary to allow the development of safer technologies and adaption of existing occupant safety regulations.
This work provides further insights into the possible effects of isometric contraction and the resulting muscle stiffness changes on injury outcome during frontal crash in the knee, thigh and hip region. Figure 2 shows an overview of the single chapters and key findings of this study.

Employed simulation framework
The finite element method (FEM) is chosen for the numerical investigations of this work. The general method can be applied to find an approximate solution for an otherwise unsolvable structural mechanical problem, such as the frontal crash simulation of this work. In the following, a very brief introduction to the basic principles of the FEM is given and a material model, which is central to the described study, is presented. For more details on the subjects, the reader is referred to standard books in literature [41][42][43][44][45].

Conservation equations and finite element method
For the following numerical investigations, the conservation of linear momentum is the essential differential equation that needs to be solved by numerical means. In its local form it is defined as: with div = divergence operator, σ = Cauchy stress tensor, k = spatial force vector, ρ = density, and a = acceleration vector. This differential equation in its local form cannot be solved for every individual material point of a body and is therefore transferred into a weak formulation. Hence, Eq. (1) will only to be solved for the integral mean of a given domain, using the principle of virtual work (2). It states that the forces of inertia, resulting from the bodies' motion and the internal forces must be in equilibrium with the external forces, with W = internal/external/kinetic energy, grad = gradient operator, δu = virtual displacements, t * = traction (stress acting on the surface of the body), dV = differential of the volume, da = derivative of the acceleration,A t = part of the surface area, where forces apply (Neumann boundary) and ' : ' as the double dot product. Using the finite element method, the domain of interest is now discretized into a number of elements with a finite element size. Each element is defined by a set of nodes with individual coordinates. From the displacements u of these nodes, the deformation of the parental element and in total, the deformation of the whole domain can be derived. Equation (2) can be transferred into a matrices equation, which can be formulated as a system of equations. By the application of shape functions N required for the interpolation of the nodal coordinates and any given point within an element, Eq. (2) can be re-written to become numerically solvable: with B ⊤ = transposed strain-displacement-matrix, ü = second derivative of the nodal displacements (nodal acceleration vector) M = mass matrix,f int = internal forces and f ext = external forces. (3) can be written in its simplified form: (4) is solved for each element in discrete time steps using the explicit time integration method. Hereby, data for the future time step (t + �t) are calculated from the data available at the current time step (t) . Since no convergence criterion has to be fulfilled for explicit time integration, the only limit is being set by the maximum critical time step assuring the conditional stability of the scheme, by using the Courant-Friedrichs-Lewy (CFL) criterion where ω max = critical eigenfrequency, c = speed of sound, L = element's length and Δt krit = maximum value of the time step, respectively [41,43,46].

Material model
Material properties of the soft and muscle tissue that are modeled using the material model MAT_SIMPLIFIED_RUBBER/FOAM were employed and strongly varied in this study. The material model can behave rubber-like or foam-like, depending on the value of the Poisson ratio v [47]. If a Poisson's ratio v between 0 and 0.5 is defined, which is the case in THUMS, then the model acts as a hyperelastic, isotropic, compressible foam-like material (MAT_SIMPLIFIED_FOAM). The Hill functional Hill is used and the Poisson's ratio ( v ) determines the penalty term of pressure in the functional. Further, the bulk modulus K is only used to determine the time step, for contact definitions and for hourglass stiffness.
The Hill functional-not to be confused with the Hill-type 1D muscle model addressed in "Background"-is defined in dependency of the stretch tensor C or the principle stretches 1 , 2 and 3 : with C j , b j and n as material constants and the Jacobian determinant J = 1 2 3 according to [47,48] describing the volumetric change. The parameter n is directly dependent on the Poisson's ratio v and can be calculated by the relation In the functional, the term 1 n J −nb j − 1 defines the penalty term of pressure. From Eq. (6), the principal Cauchy stress components σ a can be derived: The eigenvectors n a result from the spectral decomposition of the stretch tensor C . By calculation of the outer product, the second order stress tensor σ can be calculated from the eigenvectors n a and from the principle Cauchy stresses σ a : Force-displacement data obtained from experiments can be used as input for the material model. If the gauge length, width and thickness are calibrated (equal to 1), engineering stress and strain curves are used as input or tables, which reference different curves for e.g., different strain rates. This material model was used and altered for variations in muscular stiffness in the buttock, thigh and hip joint muscle and soft tissue elements of THUMS [41,[49][50][51].

Results
A general overview of the frontal impact simulation is shown in Fig. 3. The events occurring over time are described in detail in Appendix.
In the following, the effect of different volumetric muscle tissue stiffnesses is shown in comparison with literature regarding resulting stress-strain curves during impact, tissue motion, computation time and injury prediction. To achieve different muscle tissue stiffnesses, the ordinate value (stress value) of the engineering stress-strain curve was scaled, which is defined within the hyperelastic material model MAT_SIMPLIFIED_RUBBER/ σ a n a ⊗ n a = J −1 a ∂� ∂ a n a ⊗ n a . FOAM used to model muscle and soft tissues in THUMS. Different values of this scaling factor of the ordinate value (SFO) were selected to achieve different stiffness states, as further described in "Methods" of this work. Each of the following subsections will be discussed in the respective subsection in "Discussion".

Volumetric muscle element stiffness
As shown in Fig. 4, the slope of the effective stress-strain curve for elements of the buttock and thigh region increases with increasing SFO values. The element stiffness is directly defined by the slope of this curve. The slope of the linear fits, calculated for the timespan of highest loading during impact, approximates the Young's modulus for each curve. The slope increases with the SFO value, but not proportionally. Arrows indicate points of maximum effective stress and strain. Literature data [25] incorporated in the model showed the highest slope in comparison [25].
In Fig. 5, Young's moduli of different human and animal muscles from literature data are compared to those calculated from numerical simulation (Fig. 4) for different material stiffness parameters. The data are subdivided into three regions of entirely independent data. From left to right, data on 'relaxed' , 'partly contracted' and 'contracted' muscle samples are shown. 'Relaxed' includes experimental data from ex vivo experiments from, e.g., isolated bovine and porcine muscle, as well as in vivo human muscle measurements without voluntary contraction. 'Partly contracted' includes muscle stiffness states at various levels (e.g., 20%, 30% voluntary contraction, lifting of 7.5 kg weight, etc.) from human volunteers. Data classified as 'contracted' includes experimental data from human volunteers, during which the muscle was contracted to the maximum voluntary level or where heavy weights were lifted (15 kg). Further details can be found in Table 1 of Appendix. Experiments during which tetanic excitation was reached in rabbit tibialis anterior muscle by nerve excitation was also assigned to the 'contracted' samples [25]. An averaged Young's modulus of different strain rates in the 'contracted' state is shown, as well as separate values of Young's moduli for different strain rates at the 'relaxed' , passive muscle state. Results from SFO0.5 and SFO1 were both classified as 'Relaxed' . SFO2 as 'Partly contracted' and SFO10 as 'Contracted' . As the Young's modulus of SFO10 was much higher than in the literature data, respective values are listed separately. The same accounts for literature values of [25], who performed their experiments at higher strain rates than the other sources, resulting in much higher Young's moduli. Data to the presented injury risk curve [25] were incorporated in MAT_SIMPLIFIED_FOAM via a Table ID and are shown as the stress-strain curve 'Myers 1998' , cf. Fig. 4. Linear fits, approximating the Young's modulus, are in range of the literature data (Fig. 5) obtained during tensile tests [25].

Tissue motion
If the material stiffness parameter (SFO value) is increased, the tissue motion caused by the frontal impact is reduced. This can be observed in Fig. 6 in the distal region of the thigh, close to the knee joint. Only optical evaluation on tissue motion based on Figs. 6 and 7 was performed.  Table 1. Additional literature data on Shear moduli are shown in Fig. 21. The majority of the literature data shown were summarized prior to this study [40]. Apart from the subdivision in different contractile states, there is no correlation on the abscissa of the data points. Distances of data points within one group are arbitrary

Computation time and model stability
Varying the stress-strain load curve (LC) from the default value (SFO = 1) hardly had an impact on the computation time (Figs. 8,9). The largest difference with approximately 8% to the default SFO1 case was found for the SFO0.5 case with a predefined time-step size of 0.001 ms achieved by selective mass scaling. Generally, the current server utilization by other users and processes as well as hyperthreading seemed to have a much higher impact on the computation time than the change in material model parameters. This was found for the four different SFO values and 'Myers' , which were calculated multiple times at various states of server utilization. Depending on the state of server utilization, computation times varied by several hours. Differences in computation time were smallest for the isolated CPU. Some of the simulations with undefined time-step size were error-terminated after a minimum computation time of 60 ms (SFO0.5, Isolated CPU). All computations passed cycle 100,000 and 400,000, which were therefore used for computation time comparison.
Results from the simplified cuboid model showed that any changes in SFO value or stress-strain curve including different strain rates (Myers)   Regarding model stability, all simulations reached the final calculation cycle of explicit time integration (160 ms simulation time) for the mass-scaled solution (dt2ms), except the SFO0.5 simulation, which was error-terminated after 125 ms due to negative volume errors in the knee region. Therefore, reducing the muscle stiffness resulted in model instabilities.

Injury prediction
Results regarding injury prediction will be presented in two subsections: injuries of bone tissues and injuries of muscle and soft tissues. Resultant forces (contact forces), first principal strain and effective plastic strains were used as injury criteria for injury prediction. Fig. 9 Computation time in percent for all stiffness states compared to the default stiffness state (SFO1). The same comparison regarding keyword and hardware settings and stiffness cases as in Fig. 8 is shown. Error terminations occurred for all simulations with undefined time step before the final termination of 160 ms. Therefore, the final cycles could not be compared regarding computation time, but only the cycle 100,000 (~ 20 ms) and cycle 400,000 (~ 55 ms). Error termination also occurred for the SFO0.5 case (dt2ms − 0.001) after 125 ms due to negative volume error in the transition area between knee and tibia. a Simulation was run on a server with undefined time step in the CONTROL_TIMESTEP keyword. b Simulation was run on an isolated CPU with undefined time step. c Simulation was run on a server with a predefined time-step size of 1 ms achieved by mass scaling

Bone tissues
Force-dependent injury prediction The seat bottom was exclusively in contact with the THUMS pelvis and proximal region of the thigh. The knee bolster was exclusively in contact with the knee. Resultant forces were obtained from contact forces as described in subsection "Injury prediction". Contact forces between the THUMS skin and the vehicle parts and the resulting injury probability are shown in Fig. 10. For the buttock, the peak resultant force increases with increasing SFO value, but not by the same factor. The highest peak force between buttock and seat was found after 75 ms for all stiffness values except the Myers stiffness value. For the latter, the peak resultant force is reached after 65 ms. Increasing muscle stiffness accordingly leads to an increasing probability of hip fracture or dislocation according to the presented injury risk curve [52]. The muscle and soft tissue stiffness changes in the hip and thigh region do not have an impact on the probability of knee injuries (Fig. 10c, d, all except purple). For SFO0.5, SFO1, SFO2 and SFO10 the stiffness of knee surrounding muscle and soft tissues were unchanged from the default settings. For the Myers stiffness case, higher tissue stiffness was defined for all volumetric muscle and soft tissue elements modeled with MAT_SIM-PLIFIED_FOAM, including the tissues surrounding the knee joint. Increasing the stiffness of knee surrounding tissues reduced the probability of knee injuries (AIS2+) from 9-11 to 2%.

Strain-dependent injury prediction
The threshold of 3% ultimate strain was hardly exceeded by any cortical or spongy bone elements when analyzed regarding effective plastic strains and first principal strains. Therefore, a lower limit of 1.5% was selected to allow a better comparison of muscle stiffness effects on the predicted injury risk of the bones. The volume of shell and solid elements that failed the threshold of 1.5% effective plastic strain or first principal strain was added up, resulting in a 'failed element volume' . Depending on the strain type and bone tissue, different impact of muscle stiffness on injury risk prediction can be found. Data on effective plastic strains of cortical bones (Fig. 11) show that the predicted injury risk based on failed element volume is increasing with increasing muscle stiffness, as well as the number of locations, where the threshold of 1.5% strain is exceeded. In contrast, the Myers stiffness case which is stiffer than SFO10 (Fig. 4) has a slightly lower injury risk based on failed element volume (Fig. 11). Body regions that exceeded the threshold were the right hipbone 1.5 mm thickness (all simulations), 3rd lumbar vertebrae posterior (all except SFO0.5), left hipbone 1.5 mm thickness (SFO2, SFO10) and the 1st lumbar vertebrae posterior (Myers).
Data on first principal strains of cortical bones (Fig. 12) show an increase in predicted injury risk based on failed element volume in the following order: SFO2, Myers, SFO0.5, SFO1 and SFO10. Affected body regions of THUMS were the left femur upper null shell (all except Myers) and the right femur upper null shell (SFO1, Myers).
Data on first principal strains of spongy bones (Fig. 13) show an increase in predicted injury risk based on failed element volume in the following order: SFO1, SFO2, SFO0.5, Myers and SFO10. SFO10 had by far the highest failed element volume. Affected body regions were the right hipbone (all), left hipbone (all except Myers), left patella (all except Myers).
Therefore, no clear tendency between muscle stiffness changes and predicted injury risk of bones of the lower part of the body based on first principal strain evaluation (Figs. 12, 13) was found in this study.

Muscle and soft tissues
For muscle and soft tissues, effective strains were used for injury prediction, as well as the Cumulative Strain Damage Measure (CSDM) for different body regions and tissues based on effective strains ("Injury prediction").
Effective strain color plots in the dorsal view are shown in Figs. 14 and 15. They were obtained after 75 ms, where for most regions, maximum loading was found. Highest strains in muscle and soft tissue were found in the periphery of the anus and in the dorsoproximal region of the thigh. In the buttock, higher strains were found for the soft tissue (outer layer) compared to the muscle tissue (inner layer). In contrast, high strains were distributed over a larger area for the thigh muscle tissue than for the soft tissue. Peak effective strain values were continuously decreasing with increasing muscle stiffness for all analyzed tissue types and body regions, which is described in detail below. The CSDM value (Figs. 16, 17) also decreases with increasing muscle stiffness for all analyzed tissue types and body regions. Therefore, it becomes evident that increasing muscle stiffness reduces the predicted muscle and soft tissue injury risk in frontal impact scenarios based on strain-dependent data.

Volumetric and muscle element stiffness
With increasing SFO values, the slope of the effective stress-strain curve is increasing for solid elements of the muscle tissue of the buttock (Fig. 4). This confirms the expected direct influence of the SFO scaling factor on the tissue stiffness, as the slope of the effective stress-strain curve is directly correlated with the element stiffness. As shown in Fig. 5, literature data on Young's moduli of muscles increase with increasing voluntary muscle tension within the same studies (indicators), e.g., LEVINSON et al. [27] or SHI-NOHARA et al. [29]. Regarding this fact, data shown in Fig. 4 are in general agreement with literature data. For the highest material stiffness parameter (SFO10), the stiffness value obtained from simulation (Fig. 5, SFO10-green lines) was much higher than the corresponding data of contracted muscle from literature (Fig. 5, 'contracted') obtained by shearwave or sonoelastography measurement (Table 1), where usually only little strain of 0.1 to 2% is required for stiffness measurement [53]. However, other literature sources [25] performed experiments at much higher strain rates (1/s to 25/s) than the other sources of Fig. 5 and Table 1 (e.g., 0.4 s −1 [26]). They show a reverse trend by having much higher Young's moduli for the relaxed (passive) than for the tensed (active) state. Muscle material properties are significantly strain-rate dependent in the passive muscle state (Fig. 5, [25,54]). However, for the active muscle state [25], differences in material properties can only be found between quasi-static and dynamic strain rates (Fig. 5), while there is no significant difference within dynamic strain rates [25]. The data of the SFO10 case are in range of experimental data at dynamic strain rates [25] for the active state, while all SFO cases are much lower compared to the passive moduli (Fig. 5).
If stress-strain data for different strain rates from experiment [25] are incorporated in the MAT_SIMPLIFIED_FOAM material model via a TABLE ID, results are in range of experimental data of tension tests regarding the Young's moduli (Fig. 5, Myers-purple lines). Therefore, strain-rate dependency of muscle tissue should be considered when muscle behavior is analyzed and modeled in future approaches in a biofidelic way. The large variety of experimental setups, muscle samples, measurement techniques and consequently Young's moduli in literature currently makes it difficult to determine specific SFO values for each muscle in THUMS. Measurements at different strain rates for different human muscles, comparable to [25], would be necessary to further develop human body models, but are of course not possible due to ethical reasons. The motivation of this work was to cover a large area of possible material stiffnesses that resembles the stiffness of different contractile states of the literature data. For future optimization, this approach would benefit from consistent experimental data on all skeletal muscles relevant for occupant safety at different strain rates in the passive and active state.

Tissue motion
As observed by Pain and Challis (2002) for a 27 year old volunteer [28], the intrasegmental tissue motion is 50% lower in a tensed arm than in a relaxed arm. Figure 10a shows a decrease in tissue motion with increasing muscle and soft tissue stiffness (SFO value). Therefore, the simulation data are generally consistent with literature data [28].

Computation time and model stability
Based on the results, no significant effect of different material parameter values on the computation time was found. The reason for the small differences in computation time is expected to be based on server utilization and IO processes rather than changes in the material model parameters. Higher computation times were found when calculating the simulation on the server, while lower computation times were found for the isolated CPU in comparison to the default stiffness case (SFO1). For both, no significant change in computation time was found. Also, for the simplified cuboid model, no significant difference in computation time was found, supporting the argument above. The model stability was analyzed based on whether or not the simulations reached the final termination cycle after 160 ms simulation time. The model stability was unaffected by the muscle stiffness changes, except for the SFO0.5 simulation, which ended with an error termination after 125 ms due to negative volume errors in the knee region. Therefore, it can be claimed that decreasing the muscle stiffness from the default value results in model instabilities. However, as the underlying material data used for SFO1 (default) were obtained from the PMHS test, it can be assumed that this should be the lowest material stiffness definition for human muscle tissues and that a lower stiffness would not be feasible.
Other material models from the default LS-DYNA library [11,24] offer a variety of options to model the mechanical behavior of muscle in more detail than MAT_SIMPLI-FIED_FOAM and consider, e.g., passive and active muscle behavior or the anisotropic properties of muscle. However, when used in automotive HBM applications, their influence on the computation time and model stability compared to default HBM material models is still unknown and might very well change, depending on the level of detail at which the microscopic and macroscopic muscle properties are incorporated. As already stated in former studies [7,55], a trade-off needs to be found between biological accuracy on the one hand and calculation speed and data evaluability on the other. It is yet unknown, if a change in muscle material models could improve the THUMS behavior in terms of biofidelity, while maintaining reasonable computation time. This aspect should be considered in future studies in the field of occupant safety and muscle modeling.

Injury prediction
Results from "Injury prediction" are discussed.

Bone tissues
According to literature, hip injuries are more frequent than knee or thigh injuries [34]. This can be confirmed based on the force-and strain-dependent injury prediction of this study (Figs. 10, 11, 12, 13).
The most frequently observed injuries due to frontal crashes are acetabular fractures [34], i.e., fractures in the hip joint. This is in compliance with the results regarding first principal strain of the cortical bones (Fig. 12). However, this cannot be confirmed using effective plastic strain-based injury prediction (Fig. 11). Although high loadings of the acetabulum were observed, the threshold of 1.5% effective plastic strain was not exceeded. Peaks were rather found for the right hipbone and the lumbar vertebrae. The second most frequent fractures are found in the pubic ramus, sacrum and the femoral head and neck region. Respective data are listed in the Crash Injury Research and Engineering Network (CIREN) [34]. Although higher loading of some of these body regions was found, they were not regions with highest peak strain values (comp. Figs. 11,12,13). The reason for this might be that the high degree of simplification of the vehicle has a significant impact on the location of predicted injuries. Results regarding injury risk might therefore not be comparable to accident statistics, as the real vehicle interior and the simplified vehicle of this study differ too much. Further, other boundary conditions, such as vehicle impact speed and deceleration might differ in the accident statistic from the setup of this study, resulting in differences in predicted injury outcome and real injury outcome.
Based on the results of this study, muscle stiffness can affect the number and location of bone fractures, but not their general occurrence during frontal impact. A detailed numerical replication of real-world accidents regarding all boundary conditions using THUMS with different muscle stiffness states would be necessary, to determine which muscle stiffness states leads to a more realistic injury risk prediction.

Force-dependent injury prediction
Based on the force-dependent injury prediction, increasing muscle stiffness leads to a higher probability of hip fracture or dislocation (Fig. 10). In contrast, increasing the muscle stiffness in muscle and soft tissues surrounding the knee joint decreases the probability of knee injuries (Fig. 10, purple line).
Although the injury risk curve was developed for frontal impacts, such as the one simulated in this study, the setup shown in this study is not exactly matching the PMHS test setup from the experimental study, where the respective knee and hip injury risk curves were determined [52]. Further, the effect of the hands partly slipping off the steering wheel in some simulations cannot be exactly quantified. This should be addressed in future studies.
In an experiment where a medicine ball was dropped on the thigh of a volunteer [56], an increase of 11% in peak impact force from the relaxed to the voluntarily contracted muscle state was observed. Likewise, an increase in peak force with increasing muscle stiffness parameter value was found between the THUMS skin and the seat [56]. An increase by 11% of the peak force was determined, but only from the SFO1 to the SFO2 stiffness case. Therefore, the SFO2 stiffness value might be considered as the correct value to model partial voluntary muscle contraction of the human upper thigh. To verify this assumption regarding tissue motion [28] and peak impact force [56], the respective impact case should be replicated numerically.

Strain-dependent injury prediction
In general, hardly any elements of the THUMS cortical and spongy bones did exceed the respective threshold of 3% [57]. Higher impact velocities might be necessary to exceed the threshold. Therefore, a lower threshold of 1.5% was selected to allow a comparison of strain-dependent bone injury risk assessment.
Depending on the selected strain type, different impacts of the muscle stiffness changes on injury risk can be obtained for cortical bone. Effective plastic strain-based injury prediction (Fig. 11) shows a rather clear trend of increasing peak loading on the cortical bone with increasing muscle stiffness. First principal strainbased injury prediction shows an arbitrary correlation of muscle stiffness and peak loading for cortical bone (Fig. 12) and spongy bone (Fig. 13).
To confirm the results from one of the different strain values over the other, additional studies are necessary, where the injury risk could be determined for isolated body parts. A comparison with experimental studies would be essential for verification, although the effect of muscle stiffness on human bone injury risk can only be considered to a certain limit due to ethical reasons. Literature on experimental animal testing might rather be considered for validation of simplified simulation setups, such as Myers et al. [25].
Comparison: force-and strain-dependent injury prediction Results from force-and strain-dependent bone injury prediction correlate to a certain extent. Similarities regarding the impact of muscle stiffness on bone injury can be found for the forcedependent injury probability of the hip (Fig. 10b) and effective plastic strain of the cortical bone (Fig. 11), as both suggest an increase in bone injury risk with increasing muscle stiffness. Results regarding first principal strains are not correlated with the force-based injury prediction.

Muscle and soft tissues
Increasing the muscle and soft tissue stiffness had a clear tendency of decreasing the injury risk in the simulation (Figs. 14,15,16,17). In Figs. 16 and 17, effective strain was chosen as injury criterion with a value of 59%, representing the lower limit of the failure threshold of 95 ± 36% [58]. Therefore, the CSDM calculation might overestimate the muscle and soft tissue injuries. The effective strains and CSDM curves always decrease with increasing muscle stiffness. Only the SFO10 and Myers stiffness reduced the effective strain in the muscle and soft tissue to a value below the threshold of 59%.

Outlook and limitations
Certain experimental studies that were compared to the frontal impact simulation results showed major differences regarding the test setup [28,56,58]. As the topic of muscle stiffness effects on occupant safety has not been investigated experimentally to the extent necessary for comparison with this study, the literature data presented are the closest approximation currently available in this field to the best knowledge of the authors. However, in future studies, the test setup mentioned above should be represented in numerical simulations using AHBM to further investigate muscle material properties and the relevance of muscle stiffness changes in a variety of loading scenarios.
The approach of this work, where muscle stiffness changes were predefined, was tested for THUMS V5 in this study. The same approach would also be applicable to other human body models, if MAT_SIMPLIFIED_FOAM is used for muscle and soft tissue modeling. A comparison between different material models regarding computation time and biofidelity would further be of interest for different muscle stiffness states in AHBM applications. Some of these models allow to define an activation curve dependent stiffness scaling. This would enable the user to define a functional dependency of the 1D and 3D muscle systems on the material card level, which would bring the AHBM muscle model closer to the biological role model. This work showed that isometric contraction and the resulting change in muscle stiffness has an influence on the AHBM behavior and affects the predicted injury outcome. Therefore, muscle stiffness changes should be considered in different fields of biomechanics, such as the automotive sector, powered-two-wheeler safety, medical engineering, ergonomics or seat comfort analysis in the future.

Conclusion
In this study, we presented an approach to consider isometric contraction of muscles and the resulting change in stiffness of muscle tissues for the human body model THUMS. Different stiffness states were predefined. As expected, the scaling factor of the engineering stress-strain curve (SFO value) can be used to scale the stiffness of volumetric muscle elements in THUMS. It was shown that stiffness changes in the buttock and pelvis region have an influence on the occupant kinematics and peak impact forces between THUMS and the simplified seat as well as on the loading and injury risk of bone, muscle and soft tissues of this body region. Although limitations regarding differences of the experimental setup from literature and the numerical setup exist, the results were in good agreement. Changes in muscle stiffness had no significant effect on the computation time and model stability. Server utilization and IO processes seemed to have a greater impact on the computation time than changes of the material model parameters. This was shown for different time-step sizes and hardware setups. Applicability of the approach of predefined muscle stiffness is given for research and industrial applications in terms of computational costs. In the future, comparative studies on the effect of different volumetric muscle material models in AHBM applications regarding biofidelity and computational cost would be of interest.

Simulation setup
The human body model THUMS, representing a 50 percentile average sized American adult male as an occupant, was used in this study. The braced muscle contraction state was predefined, during which the THUMS is pushing itself off the steering wheel and is stepping on the break. The most relevant vehicle parts (steering wheel, seatbelt, inflatable airbag, seat, footrest and kneebolster) were modeled as simplified and rigid, except for the airbag and seatbelt, see Fig. 20 in Appendix. The HBM was restrained by a seatbelt using both 1D seatbelt elements to include sliprings, pretensioner and retractor nodes efficiently and 2D shell elements for the contact to THUMS during the crash pulse. By accelerating the vehicle backwards according to a predefined acceleration curve (Fig. 18), a simplified frontal impact scenario is simulated. Occupant kinematics during impact are shown in Appendix (Fig. 20).

Muscle stiffness parameters
The material model is, among others, defined by the linear bulk modulus (KM), the shear modulus (G), an engineering stress-strain load curve (LC), a damping coefficient (MU) and a limit stress value (SIGF). The LC is only defined as such if the specimen dimensions are calibrated (equal to 1) under uniaxial loading. The variable which showed the most promising results in terms of stiffness alteration was the scaling factor of the ordinate value (SFO) that scales the engineering stress value of the load curve (LC). The default tissue stiffness and, accordingly, the default LC shape of MAT_SIM-PLIFIED_FOAM in THUMS V5 was defined to represent overall PHMS behavior. The variety of scaling factors (SFO values) was selected in a way that a wide range of different muscle stiffness states is covered. The four SFO values were chosen as: (1) SFO = 0.5 for reduced stiffness, (2) SFO = 1.0 for default stiffness (THUMS V5, PMHS data [2,14]), (3) SFO = 2.0 for slightly increased stiffness (possibly representing partial, voluntary contraction) and (4) SFO = 10.0 for highly increased stiffness (possibly representing maximum, tetanic contraction). The muscle stiffness was altered for muscle and soft tissue elements of the buttock and thigh that came into contact with the seat bottom (green and yellow parts in Fig. 19). This body region was analyzed and compared with literature data, mainly Pain and Challis and Tsui and Pain [28,56]. To determine the effect of the SFO value on the material stiffness and for comparison with literature data, the effective stresses and effective strains were analyzed. Additionally, a fifth stiffness case was defined, where strain-rate-dependent data were incorporated in the model (Myers) via a Table ID using literature data obtained from experiments with rabbit muscles [25]. This Myers stiffness case was incorporated in all muscle and soft tissue elements of THUMS, not only for the buttock and thigh. The effective stress, also known as von Mises stress, is defined as: where σ 1 , σ 2 and σ 3 are the principal stresses. The effective strain, expressed in tensorial notation, is defined as:

Simulation environment and computation time
The equation of conservation of momentum (subsection "Conservation equations and finite element method", (1)) has to be solved numerically with constitutive models using the explicit time integration method for the loading scenario described above. By application of the explicit FE code MPP LS-DYNA v971 revision 7.1.2_95028 single precision (LSTC), the numerical solution was calculated using 16 cores on a Supermicro Computer (3.20 Ghz Intel ® Xeon ® CPU E5-1680 v3 processors with 128 GB RAM using Melanox Infiniband) running on CentOS Linux release 7. To achieve reasonable computation times, the time-step size for mass-scaled solutions in the LS-DYNA CON-TROL_TIMESTEP [59] keyword was set to 0.001 ms by selective mass scaling 2 for all simulations.
To analyze the impact of muscle stiffness on the computation time, different keyword and hardware setups were compared. On the one hand, the dt2ms value of the CON-TROL_TIMESTEP was set to -0.001, thus fulfilling the CFL criterion by mass scaling. For all dt2ms = − 0.001 simulations, the mass of 85.29 kg of the entire simulation setup increased by 2.67% (2.28 kg). This was the standard setup for the simulation of this work, used for analysis of the injury outcome, tissue motion and peak impact force and the material behavior. On the other hand, dt2ms was left blank, thus achieving the CFL criterion by altering the minimum time-step size in dependence of the smallest element size. Regarding hardware setups, simulations were calculated on the server, where other users could perform simulations at the same time, using the setup mentioned above with 16 cores. For comparison, simulations were also calculated on an isolated CPU with 8 physical cores without hyperthreading one after another. Besides the four SFO stiffness simulations, the SFO10 value and strain-rate-dependent data (Myers) were applied to all volumetric MAT_SIMPLIFIED_FOAM muscle and soft tissue parts to generate 'worst case scenarios' regarding computation time.
The influence on computation time was further analyzed in simplified simulations, where a cube model consisting of eight cuboid solid elements (edge length: 1 mm) modeled with the MAT_SIMPLIFIED_FOAM model from THUMS buttock and thigh muscle and soft tissues. No loading was applied. The simulation ran on the isolated CPU for 50 ms with no defined time-step size.
The influence of muscle stiffness on the THUMS model stability was analyzed based on whether or not the simulation reached the final calculation cycle of explicit time integration for mass-scaled solutions (dt2ms). The termination time, meaning the predefined simulation time after which the simulation automatically stops, was set to 160 ms. Trube et al. BioMed Eng OnLine (2021) 20:53 Literature data of muscle tissue stiffness measured using elasticity imaging Figure 21 contains supplementary data on Shear Moduli G besides the data on Young's moduli E from Fig. 5. Detailed information, including the type of samples and measurement technology are listed in Table 1. Many of the literature sources presented were previously summarized by Sarvazyan et al. [40].  Table 1. Apart from the subdivision in different contractile states, there is no correlation on the abscissa of the data points. Distances of data points within one group are arbitrary   A variety of measurement technologies that can all be assigned to the field of 'Elasticity Imaging' or 'Elastography' [16] were used to obtain this data. Data were collected and assembled from publications listed in the column 'Source' by the author of this work