Implementation and validation of the extended Hill-type muscle model with robust routing capabilities in LS-DYNA for active human body models.

BACKGROUND: In the state of the art finite element AHBMs for car crash analysis in the LS-DYNA software material named *MAT_MUSCLE (*MAT_156) is used for active muscles modeling. It has three elements in parallel configuration, which has several major drawbacks: restraint approximation of the physical reality, complicated parameterization and absence of the integrated activation dynamics. This study presents implementation of the extended four element Hill-type muscle model with serial damping and eccentric force-velocity relation including [Formula: see text] dependent activation dynamics and internal method for physiological muscle routing.
RESULTS: Proposed model was implemented into the general-purpose finite element (FE) simulation software LSDYNA as a user material for truss elements. This material model is verified and validated with three different sets of mammalian experimental data, taken from the literature. It is compared to the *MAT_MUSCLE (*MAT_156) Hill-type muscle model already existing in LS-DYNA, which is currently used in finite element human body models (HBMs). An application example with an arm model extracted from the FE ViVA OpenHBM is given, taking into account physiological muscle paths.
CONCLUSION: The simulation results show better material model accuracy, calculation robustness and improved muscle routing capability compared to *MAT_156. The FORTRAN source code for the user material subroutine dyn21.f and the muscle parameters for all simulations, conducted in the study, are given at https://zenodo.org/record/826209 under an open source license. This enables a quick application of the proposed material model in LS-DYNA, especially in active human body models (AHBMs) for applications in automotive safety.

Background

It is hard to imagine that the development of modern vehicles can be done without wide utilization of computer aided engineering (CAE), including state-of-the-art numerical simulations software. These tools not only offer opportunities to the car manufacturers for saving time and costs during the design phase, but also to model and predict future product lifecycle. One of the most demanding and regulated domains is vehicle safety and therefore crash simulations. For more than a quarter of a century, complete vehicles are modelled virtually as finite element models with all significant details, including material and geometrical properties. The same approach was then applied to the human body, and in the last decade several detailed finite element human body models (HBMs) were presented [1-3]. Joint simulations with a combined application of car and human body models allow for prediction of in-crash behaviour and possible injuries for occupants or pedestrians with a sufficient accuracy throughout all stages in development, replacing expensive crash tests using car prototypes and Anthropomorphic Test Devices (ATDs).

These so-called virtual testing methods will gain importance in the future, as the current trends of active safety systems and autonomous vehicles become available on the market. Active safety systems, in contrast to traditional passive safety systems, react preventively prior to a crash to avoid or mitigate a possible impact. This requires active HBMs (AHBMs) that are able to reproduce human behaviour in normal driving situations, as well as the behaviour in the in-crash phase. The same requirements exist for the second trend of autonomous vehicles, where generic driving or sitting positions no longer exist, but where occupants can move freely. From these requirements, three challenges for AHBM modelling arise. The first being the implementation of active muscles as mathematical models of the muscle-tendon complex (MTC) including the activation dynamics, which will be addressed in the contribution. The second challenge is to model biologically relevant neural controllers to enable accurate forward dynamics (FD) simulations of human reactions and voluntary motion in all kind of traffic scenarios. The third challenge is the choice of parameters for AHBMs, as only a correct representation of both the passive components and active components will result in an accurate representation of a living human. Most current HBMs have a passive stiffness which is too high, see e.g. [4, 5]. On the one hand, this is to compensate for the active components still missing. On the other hand, this is because the source for the parameters are almost exclusively post mortem human subjects, where tissue modulus and no-load strain differ from living tissue depending on the postmortem time and post-mortem rigor [6].

In the state of the art finite element AHBMs [7, 8] for car crash analysis the LS-DYNA material named *MAT_MUSCLE (*MAT_156) is used for modeling active muscles. This material is an advanced version of the previous model *MAT_SPRING_MUSCLE [9] for discrete elements, that is no longer being supported. *MAT_156 represents a Hill-type muscle model which consists of three parallel elements: contractile element (CE), parallel elastic element (PEE) and a parallel damping element (PDE), see also Fig. 1a. The implementation was done by Dr. J. A. Weiss based on prior studies and reviews on different Hill-type model element configurations by [10-12]. The implemented configuration was chosen due to its simplicity, ease of parameters derivation from the experiments and computational efficiency. However, in the publication [12] it was pointed out, that an element configuration with better approximation of the physical reality should be used in simulations if possible. Such an extended Hill-type muscle model should have a clear separation between muscle fibres and tendon structures. For a correct representation of the MTC dynamics an additional internal degree of freedom is required to decouple active muscle fibre and elastic tendon dynamics. Subsequent studies investigating the role of the serial elastic element have shown, that such simplifications and assumptions can lead to instabilities produced by force-velocity or force-length relation formulations [13], incorrect energy storage and release in the interaction with the environment [14, 15], unrealistic high-frequency oscillations [16] and differences in muscle force magnitude [17]. All these effects, mentioned in publications above, directly influence the explicit integration scheme used in LS-DYNA thus impacting speed, accuracy and robustness of simulations with AHBMs.

Fig. 1

Schematic structures of the Hill-type muscle models. With CE contractile element, PEE parallel elastic element, PDE parallel damping element, SEE serial elastic element, SDE serial damping element. a Structure of the LSDYNA muscle model *MAT 156. b Structure of the implemented four element muscle model

Usually, muscles and tendons wrap around bones or joints in both steady state conditions and while performing movements, consequently a physiological muscle path representation (muscle routing) is essential for FD simulations [18, 19]. Slight changes in the muscle line of action will lead to inaccurate muscle forces and resulting moments due to incorrect lever arms and muscle length. To model physiological muscle paths in finite element HBMs different muscle routing methods can be used. Fixed lever arms or the via-point method [20, 21] are the most simple options and the usage of contact detection [22] would be the most sophisticated method. According to [23] a via-point method should be preferred for *MAT_156 using the *ELEMENT_BEAM_PULLEY keyword in LS-DYNA. However, it is unclear if this method is applicable, as there exist so far no successful implementation of this method in AHBMs to the authors knowledge.

Additional disadvantages result from the way parameters are set for *MAT_156 in LS-DYNA. For a number of parameters predefined curves are required, e.g. muscle activation level vs. time or stress vs. the stretch ratio. These curves need to be defined beforehand or might be calculated during the runtime through the *DEFINE_CURVE_FUNCTION keyword and PIDCTL [24] options. This approach is limited, cumbersome and error prone. Instead, muscle parameters and constants found in anatomico-physiological literature should be used directly. Also, some disadvantages exist for muscle activation dynamics. Predefined muscle activation level vs. time curves cannot represent the activation dynamics correctly and to model the dynamics more precisely the activation level has to depend on the muscle length. The complete activation dynamics can be included efficiently in the material model for the muscle itself.

This study addresses the problems mentioned above and presents an implementation of the extended Hill-type muscle model with serial damping and eccentric force-velocity relation proposed by Haeufle et al. [25] into LS-DYNA code. The muscle model is additionally extended by activation dynamics and a method for physiologial muscle routing. The complete description of the model, its implementation, verification and validation are given in the next sections.

Methods

In this section the complete model, its implementation in LS-DYNA, and the verification and validation set-up are described.

Muscle model

One of today’s most popular and widely used macroscopic muscle model was proposed by Hill in 1938 [26] on the basis of experiments with frog muscles. The most important feature is a direct relation between muscle force and contraction velocity. Furthermore the model is also referred to as a first order muscle model, which means, that the muscle elements have neither mass nor inertia, and only an axial force is applied on the skeletal model between origin and insertion point of the muscle. During the past years some disadvantages of the original Hill model were found and there have been many publications with the aim of further developing and improving this model.

In the publication of Haeufle et al. [25], a modified Hill-type muscle model was proposed with improved serial damping and eccentric force-velocity relation. This model consists of four simple mechanical elements: an active contractile element (CE), which is controlled by the activation level q; parallel (PEE) and serial (SEE) nonlinear spring elements and a serial damping element (SDE). The model’s structure was based on a previous study by Günther et al. [16] which determined, that the model with a force-dependent SDE provides the best results in a comparison with constant parallel, constant serial, and force-dependent parallel damping elements. The structure of the MTC of the Haeufle model is shown in Fig. 1b, with two clearly separated parts modeling the active muscle fibres (CE + PEE) and the passive tendon and aponeurosis structures (SEE + SDE). The main equations of the muscle model are presented in the following. They were taken from [16, 19, 25, 27], where more detailed explanations can be found if needed. Furthermore, a comprehensive study on the influence of individual parts and their model formulation is given in [28].

As shown in Fig. 1b the muscle model features an internal degree of freedom which is described by  . The lengths of the passive elements are equal to Then the total MTC length isThe force equilibrium at point P between the muscle fibre and the tendon part is described in [16] as:

Contractile element (CE)

The contractile element represents the active fibre bundles of the MTC. The force of the contractile element  is therefore dependent on the muscle activity q, the contraction velocity  as well as the length-dependent isometric force . It is expressed by the equationThe factors and are so-called normalized Hill parameters, where is normalized with the maximum isometric force and with the optimal fibre length  [16,  p. 64]. The subscript ‘rel’, for relative, indicates the normalization. The optimal muscle fibre length at which the isometric force reaches the maximum value is  . The isometric force  depends on the length of the contractile element and is calculated as follows:This equation represents the bell-shaped force-length relationship of the CE element. The width of the normalized bell curve  and the exponent  may be chosen differently for the ascending and descending limb of the force-length curve.

When calculating the Hill parameters, it is distinguished between an eccentric (lengthening fibres) and a concentric case (shortening fibres). Please note, that in physiological muscle experiments, where shortening work of muscle fibres is examined, the sign convetion for the contraction velocity is usually the opposite (shortening fibres have a positive velocity) to ensure that the work of shortening muscles is positive. In the concentric case, the Hill parameters are: The auxiliary variables for the calculation of the Hill parameters are divided into length- and activation-dependent components. The length-dependent parameters are defined as:and the activation-dependent as: The equations in the eccentric case can be derived from Eq. (5) as mentioned in [25] and thus they are formulated as:and

Parallel elastic element (PEE)

The PEE represents passive properties of the muscle fiber and the collagenous connective tissue surrounding the muscle belly. As soon as the length of the contractile element exceeds the resting length of the parallel elastic element , it also contributes to the force developed by the MTC. Mathematically this is expressed as:The spring stiffness is influenced by the optimal fibre length, the width of the bell curve and the maximum isometric force. It is calculated by:The resting length is defined as , hence is the resting length normalized by  , is width of on a descending limb.

Serial elastic element (SEE)

Since structures similar to muscle tissue are also present in the tendon, their elastic properties are similar. The serial elastic element has a nonlinear or linear spring behaviour depending on the deflection . When the tendon is relaxed and does not generate any force. In the range of it has a nonlinear characteristic, and a linear characteristic for :The length  of the SEE at the transition from nonlinear to linear characteristic, the exponent , and the nonlinear and linear stiffness factors  and are defined by the following formulas:The complete description of these independent parameters can be found in [16, Fig. 4, p. 69].

Serial damping element (SDE)

The force-dependent serial damping element reduces unphysiological high-frequency oscillations in the tendon part of the muscle model. As a side effect this also increases numerical efficiency [16]. The force-dependent damping of the material damping element is calculated as:with the maximum absorption value ofusing the dimensionless scaling factor and minimum damping value  .

Contraction dynamics

Inserting Eq. (5) into Eq. (4) for the force equilibrium yields a quadratic equation of the form:This equation must be solved for the contraction velocity at each time step. Subsequent integration gives the solution for the internal muscle model degree of freedom—the length of the contractile element . Since the coefficients and are always less than zero for our configuration, the solution for the contraction dynamics is given as:In this equation, the index e denotes that eccentric Hill parameters must be computed from Eqs. (12, 13). The coefficients , and are determined as follows:with the additional coefficient

Activation dynamics

In the application of muscle models, not only the muscle dynamics itself, but also muscle activation dynamics needs to be considered. Activation dynamics is the link between stimulation input from the nervous system and the activity level of a muscle. For the proposed muscle model, two different muscle activation strategies are implemented: one depending only on the neural activation level (STIM) by Zajac [11] and another, which takes into account length-dependent sensitivity of  level change by Hatze [29]. These two activation dynamics are outlined below.

The first implemented activation dynamics by Zajac [11] was extended by [16] by adding a minimum muscle activity level to represent the fact that in reality a muscle is never physiologically completely inactive (). The differential equation for the activation dynamics therefore is noted as:In this equation STIM is the input. It is the neural stimulation that emanates from the nervous system and varies from 0 to 1. The output is q, the CE element activation level with a possible range of . It represents the concentration of free  ion in the muscle. is the activity time constant and is the ratio between time constant on activation and deactivation. Thus, for the deactivation time constant is less than that of the activation.

The second activation dynamics implemented is a two-step approach introduced by Hatze [29]. In this approach, the activity level q depends on both the length of the contractile element  and the free  ion concentration. The activity level is calculated as follows:The  ion concentration is accounted for in the differential equation as :and the relative CE length is included in :Here m, c and are constants and is the ratio between the contractile element length and the optimal fibre length . Thus the length-dependent  ion sensitivity is taken into account, namely the relation that the longer the contractile element the higher the sensitivity. In other words, stretched muscles produce a larger force at the same stimulation level compared to an already contracted muscle [30]. In addition, the sensitivity contributes to low-frequency stiffness of the muscle, which is defined as the change in the equilibrium muscle force relative to a change in the equilibrium length with constant stimulation [31].

Muscle length offset and muscle routing

To enable physiological muscle path representation for the extended Hill-type muscle model several routing methods could be considered. In advanced modelling frameworks the muscle path is usually redirected either by specific points, so called via-points, or by surfaces of geometrical objects (e.g. in OpenSim [32]). See [33] additionally for an in-depth review and comparison of routing methods in biomechanical models. It was decided to use the via-point approach as described in [20], because it is possible to implement this method with the standard routing elements available for seatbelts in LS-DYNA. This method has proven to be reliable, as it is used in almost every crash simulation involving occupant models. To implement this, it is necessary to divide the MTC into muscle element and seatbelt elements, as only the latter can be routed. Therefore, an offset in length  is introduced, defined as the difference between the actual length of the muscle beam element length  in the model and the length of the entire MTC :If necessary an offset can be added on both ends of the muscle beam element, to allow for two or more via-points, see Fig. 2. Standard seatbelt elements can be attached to the end of the muscle beam element and all the standard routing methods of LS-DYNA e.g. sliprings can be used. The seatbelt elements can move through a slipring node freely, while at the same time, the muscle model internally works with the correct length and dynamics of the entire MTC. To preserve the muscle dynamics it is required, that the stiffness of the seatbelt elements is orders of magnitude higher compared to the stiffness of the muscle elements. For the example in "Application in the ViVA OpenHBM Arm with routing" a stiffness of  N/m was used for the seatbelt elements.

Fig. 2

Comparison of the length relations for the muscle routing approach with the via-point method. a A full beam element with Hill-type muscle material, where the beam element represents the entire MTC. b A shortened beam element with Hill-type muscle material extended by seatbelt elements, and c the via-point routing method with two via-points. For the latter two, the muscle force is still calculated based on the entire MTC length (muscle + seatbelt), however, it is acting only in the beam element. This approach allows to use the slipring routing method of LS-DYNA

LS-DYNA implementation

It is possible to include self-written code in the LS-DYNA FE solver through so-called ‘User Subroutines’. These subroutines have to be written in FORTRAN and can, among other options, be used to define user materials [24]. The muscle model described in "Muscle model" was implemented in LS-DYNA as a user material for truss elements to simulate the active contraction behaviour, as well as the passive spring and damping effects of human muscles. The FORTRAN code is available at https://zenodo.org/record/826209.

The explicit integration scheme in LS-DYNA, shown in Fig. 3, is updating the element strain  in each timestep based on the nodal displacement . Material models translate the strain  to stress , which yield nodal forces . These forces result in nodal acceleration , which are integrated to nodal velocity  and displacement  for the next timestep.

Fig. 3

Implementation of the user material subroutine into LS-DYNA workflow

It should be pointed out, that material subroutines require an element stress as a return value. In the concept of the muscle model, only forces are calculated. Since truss elements can only have axial stress, the stress was calculated from the muscle force and the element cross-section area viaIf the material card for user-defined material models is specified in an input deck, LS-DYNA internally calls the routine usrmat, which starts the corresponding element routine, depending on the element type. In the case of beam elements this is urmatb and for truss elements urmatt. Finally, the actual material routine is called, which the user can program himself. It is possible to have up to ten user materials defined in the subroutines umat41 to umat50. The user can implement arbitrary material models in these routines and, among other things, access the material parameters specified in the material card. In addition, the programmer may call further subroutines, which then return, for example, nodal coordinates or various element properties.

Verification and validation set-up

For the Hill-type model parameters identification, a general test procedure requires three experimental set-ups: (a) concentric contraction, (b) isometric contraction and (c) quick release [12]. They are depicted in Fig. 4 and explained in detail below. A complete set of muscle parameters is almost never found in a single source since it is hard to perform all three tests in a row with the same muscle specimen, so a short literature survey is always required.

Fig. 4

Illustration of the a concentric and b isometric contraction and the c quick release experiments

The validation conducted for the muscle model is based on mammalian muscle experiments. As there is no experimental data available for actual human muscle tissue, the model validation is based on piglet [16], cat [34] and rat [35] muscle experiments. The verification is done in comparison with an already existing implementation of the same muscle model in the Matlab based multi-body code Neweul-M2 [36]. Additionally, a comparison with the *MAT_156 muscle material from LS-DYNA is shown for the concentric contraction experiment. Data sets from all three experimental set-ups are available for the piglet muscle. For the other two species, a specific set-up for an isometric contraction experiment is applied, which was shown to be sufficient to determine all necessary Hill-type parameters for simulations [37]. An overview of all experimental set-ups is presented in Table 1. Furthermore all model parameters are given in tabular form for each validation case including references.

Table 1

Overview of all set-ups used for validation

	Piglet	Cat	Rat
Concentric contraction	X
Isometric contraction	X	X	X
Quick release	X

Concentric contraction experiment

In a concentric contraction, the muscle is shortened, which means that the distance between muscle origin and insertion point decreases, e.g. elbow flexion to lift a weight, see Fig. 4a. In the simulation set-up, the muscle element is orientated vertically and the upper node is fixed. Masses between  kg and  kg are attached to the lower node in accordance with the experimental studies. At the beginning of the test, the muscle is relaxed and the mass is resting on a plane. Then the muscle is stimulated () and starts to contract. At first no external motion is recorded as only the internal length is decreasing. Once the mass is lifted and the contraction velocity is recorded and compared to the experimental results from the piglet muscles.

Isometric contraction experiment

In an isometric contraction the muscle force is increased, while the length of the muscle is kept constant, see Fig. 4b. This contraction mode occurs, for example, when attempting to hold a heavy weight. In the simulation set-up both nodes of the muscle element are fixed. At the beginning of the test, there is no stimulation (), thus the muscle experiences only a minimal activity . Starting from 0.1 s, the muscle is stimulated completely () and relaxed again completely after 1.1 s () in the piglet [16] and cat [34] experiments. In the rat experiments the muscle is only activated for a shorter time period of 300 ms [35]. In the piglet muscle experiments the isometric contraction was carried out for different fixed muscle lengths between 5.1 and 6.6 cm around the anatomical resting length of . In the other experiments the isometric contraction was only tested for the anatomical resting length. In the results, the force vs. time curves are compared and also the differences resulting from the two activation dynamics implemented are analyzed.

Quick release experiment

The quick release is a combination of isometric and concentric contraction. In this set-up, the muscle is fixed at both ends at the beginning, it is then stimulated (STIM = 1) and isometric contraction occurs Fig. 4c. After 1 s the lower end of the muscle carrying a mass is released and is pulled up quickly due to the force built up during the isometric contraction. After a total time of 1.5 s the stimulation is switched off again (STIM = 0). As in the concentric case, the influence of the different masses is examined (m = 200, 400, 600, 800, 1000, and 1500 g) for the piglet muscles only [16].

Simulation results

The verification and validation simulation results are shown in the following sections for piglet, rat and cat muscles. Using the piglet data, an additional comparison to the muscle model *MAT_156 already existing in LS-DYNA is shown and the differences resulting from the two distinctive muscle activation dynamics implemented are illustrated. To demonstrate the application of the model in AHBMs an example illustrating the routing capabilities is given using an elbow model extracted from the ViVA OpenHBM [3].

Piglet calf muscle

For the piglets calf muscles results for concentric and isometric contraction and quick release experiments are available in [16]. As this is the most complete data set, it was also used for verification and a comparison with *MAT_156 and a comparison of the different activation dynamics available in the extended Hill-type muscle model. In Table 2 the parameters used for the piglet simulations are listed. The material card for LS-DYNA is found in Appendix "Material card for piglet simulations".

Table 2

Muscle parameters for the piglet simulations. See [16, Table 2, p. 68]

Activation dynamics(Zajac [11])	q_0[] 1.0e−4	τ_q[s] 0.025	β_q[] 0.5
Activation dynamics(Hatze [29])	q_0[] 5.0e−3	c[ ]1.373e−4	η[] 5.27e−4	k[ ]2.9	m[ ]11.3
Isometric force	F_max[N] 30.0	l_CE,opt[m] 0.015	ΔW_des[] 0.14	ν_CE,des[] 3.0	ΔW_asc[] 0.57	ν_CE,asc[] 4.0
Force-velocity hyperbola	A_rel,0[] 0.1	B_rel,0[1/s] 1.0	S_e[] 2.0	F_e[] 1.8
PEE	L_PEE,0[] 0.9	ν_PEE[] 2.5	F_PEE[] 1.0
SEE	l_SEE,0[m] 0.045	ΔU_SEE,nll[] 0.1825	ΔU_SEE,l[] 0.073	ΔF_SEE,0[N] 60.0
SDE	D_SDE[] 0.3	R_SDE[] 0.01

Concentric contraction

The numerical and experimental results are presented in Fig. 5. All curves are shifted in time so that the mass is pulled up or occurs at . As shown in the figure, the simulation results from LS-DYNA are very consistent with the experiments and the simulations with the muscle model in Neweul-M2. A comparison with the simulation data from [16] would give even better results. The differences between the simulation results from Neweul-M2 and LS-DYNA can presumably be attributed to different computational accuracies and integration methods for the differential equations. These simulations were both run with explicit integrators and a constant time step. Consequently, we can state that we have a correct muscle model implementation for the concentric contraction case.

Fig. 5

Concentric contraction velocity of the MTC of a piglet over time at different muscle loads. Full line LS-DYNA, dots Neweul-M2, dashed line experimental results from [16]

In addition a comparison with the muscle material model *MAT_156, already existing in LS-DYNA, is shown in Fig. 6. The initial contraction velocity provides similar results. Also, the maximum force value for low masses up to about 200 g is well approximated. The *MAT_156 material, however, shows significant weaknesses in speed decay and in the correct representation of the damping properties. At this point, it should be noted that further optimization of the material parameters might achieve better results. For these simulations the *MAT_156 parameters were derived from the previous work of [38]. This comparison shows the potential of the newly implemented muscle model and shows how it can help to deliver more realistic simulation results.

Fig. 6

Comparison of the concentric contraction velocity between the extended Hill-type muscle model and *MAT_156. Full line extended Hill-type muscle model, dashed line *MAT_156. Colours are identical to Fig. 5

Isometric contraction

In Fig. 7 the numerical results for both the LS-DYNA and the Neweul-M2 simulations are depicted together with the experimental results. In the piglet experiments different starting length of the MTC  were tested. In Fig. 7 the tests are differentiated by the stretch ratio of the starting length and the anatomical resting length .

Fig. 7

Force output of the MTC plotted versus time for different fixed stretch ratios h in isometric contraction. Full line LS-DYNA, dots Neweul-M2, dashed line experimental results

The comparison of the LS-DYNA results with the experimental data shows, that the muscle force for the inactive muscle in the time intervals and is underestimated if the muscle is lengthened considerably relative to the anatomical resting length (). Also, a clear deviation in the force increase for the stretch ratios and exists, while the final force at is met. Very similar differences are also present in the simulation results from [16]. According to this source, the proposed muscle model does not represent all internal dependencies of the Hill parameters correctly. Also potential history effects visible in the experimental curves, namely non-steady force plateaus, are made responsible for the differences. Furthermore, possible deficits in the identification of the parameters for the activation dynamics and the rise of play a role. The comparison to Neweul-M2 shows high agreement with only slight deviations in the muscle activation interval for the ratios and 0.91. As it was the case in the concentric contraction, the differences in the results are larger for higher dynamics, which can probably be attributed to the different integration method for solving the differential equations.

The most important point illustrated by the isometric contraction in Fig. 7 is the strong dependence of the maximum isometric force on the muscular length. This finding is decisive for the application in AHBMs, since in this example deviations of approximately 15 mm in muscle length lead to differences in muscle force of more than 30 N.

Comparison of activation dynamics for isometric contraction

In the extended Hill-type muscle model two different approaches to describe the activation dynamics are implemented. In Figs. 8 and 9 these methods by Zajac [11] and Hatze [29] are compared in an isometric contraction.

Fig. 8

Comparison of muscle forces for Zajac and Hatze activation dynamics. Colours identical to Fig. 7

Fig. 9

Comparison of the activation level q for Zajac and Hatze activation dynamics. Colours identical to Fig. 7

The muscle force results with Hatze and Zajac activation differ mainly during muscle deactivation after , see Fig. 8. The forces of muscles with a high stretch ratio decrease significantly later than the muscles that are shortened. The concentration of free  ions evolves similar to the activity level described by the differential equation of Zajac. increases slightly slower and takes about 0.1 s longer to decay. The larger influence is the dependence of the activation on the CE length through for the Hatze activation dynamics, see Fig. 9. By stretching the muscle, is significantly longer for high h-ratios. As a result, will increase and the activity is rising faster for the stretched muscles. The stretch ratio and thus also affect the maximum activity level reached in the simulations. It can be seen in Fig. 9, that the maximum activity of the muscle with is only about . For Zajac’s activation dynamics only one curve is found in Fig. 9, this is because Zajac’s activation dynamics is length-independent and therefore all activation curves are identical.

As the formulation of activation dynamics by Hatze is the more biofidelic and superior option [37], the comparison for rat and cat experiments will be done only with this dynamics.

Quick release

The quick release experiments are a combination of the two experiments above. Here the force produced by the MTC is analyzed versus the contraction velocity. In Fig. 10 the numerical results from LS-DYNA and Neweul-M2 as well as the experimental data is shown.

Fig. 10

Force output of the MTC plotted versus contraction velocity at different muscle loads in quick release experiments. Full line LS-DYNA, dots Neweul-M2, dashed line experimental results

The isometric muscle force at zero velocity, i.e. before the muscle is released, is about 2 N lower than in the experiments. However, the results clearly approach the respective maximum contraction velocities. In [16] it is stated, that this is due to history effects within the tendon in the experiments that are not represented in the muscle model. The best agreement is achieved for a mass of 1000 g, where according to [16] the history effects were absent. The difference with Neweul-M2 is once again negligibly small for the bigger masses and slightly larger for the high velocities or smaller masses.

Rat gastrocnemius medialis muscle

The experiments were done by Siebert et al [35] on the rat (Rattus norvegicus, Wistar) M. gastrocnemius medialis muscle. The parameters for the Hill-type muscle model are listed in [35, Table 4, p. 222] and the experimental set-up description is provided on page 218 of the same publication. Optimal Hatze activation dynamics parameters are listed in [37, Table 2, column 6, p. 278].

For convenience, they are collected in Table 3 and the LS-DYNA material card is given in Appendix "Material card for rat simulations". As seen in Fig. 11 the simulation results are in good agreement with the experimental results, being a little faster in the muscle deactivation slope.

Table 3

Muscle parameters for the rat simulations

Activation dynamics(Hatze [29])	q_0[] 6.0e−3	c[ ]1.373e−4	η[] 5.27e4	k[ ]2.9	m[ ]22.54
Isometric force	F_max[N] 11.2	l_CE,opt[m] 0.0148	ΔW_des[] 0.35	ν_CE,des[] 1.5	ΔW_asc[] 0.35	ν_CE,asc[] 3.0
Force-velocity hyperbola	A_rel,0[] 0.06	B_rel,0[1/s] 1.42	S_e[] 0.99	F_e[] 1.35
PEE	L_PEE,0[] 3.12	ν_PEE[] 2.5	F_PEE[] 2.0
SEE	l_SEE,0[m] 0.0123	ΔU_SEE,nll[] 0.0425	ΔU_SEE,l[] 0.017	ΔF_SEE,0[N] 4.48
SDE	D_SDE[] 0.3	R_SDE[] 0.01

Fig. 11

Isometric experimental and simulation results for a rat muscle

Cat soleus muscle

Mörl et al. [34] conducted the experiments on the cat soleus muscle. The parameters for the Hill-type muscle model are found in [34, Table 1, p. 5] with optimal Hatze activation dynamics parameters once more taken from [37]. They are also collected in Table 4 and a material card for LS-DYNA is provided in Appendix "Material card for cat simulations". The corresponding simulation results depicted in Fig. 12, are in excellent agreement with the experimental results, this time being a little faster in the muscle activation slope.

Table 4

Muscle parameters for the cat simulations

Activation dynamics(Hatze [29])	q_0[] 1.0e−4	c[ ]1.373e−4	η[] 5.27e4	k[ ]2.9	m[ ]22.54
Isometric force	F_max[N] 10.0	l_CE,opt[m] 0.053	ΔW_des[] 0.35	ν_CE,des[] 1.5	ΔW_asc[] 0.35	ν_CE,asc[] 3.0
Force-velocity hyperbola	A_rel,0[] 0.07	B_rel,0[1/s] 0.2	S_e[] 2.0	F_e[] 1.5
PEE	L_PEE,0[] 0.9	ν_PEE[] 2.5	F_PEE[] 2.0
SEE	l_SEE,0[m] 0.060	ΔU_SEE,nll[] 0.0425	ΔU_SEE,l[] 0.017	ΔF_SEE,0[N] 4.0
SDE	D_SDE[] 0.3	R_SDE[] 0.01

Fig. 12

Isometric experimental and simulation results for a cat muscle

Application in the ViVA OpenHBM Arm with routing

The extended Hill-type muscle model is applied in an arm model extracted from the ViVA OpenHBM [3]. The model includes bones, modelled as rigid bodies, and flexible flesh and skin of the upper extremity. We added the main flexors (biceps long and short heads, brachialis, brachioradialis, pronator teres and extensor carpi radialis) and extensors (triceps long, lateral and medial heads) of the elbow joint, which we idealized as a revolute joint. Here the via-point routing method is compared to simple direct line and lever approaches. Moreover a complete description of the set-up of the elbow model [39] and the choice of parameters for all muscles at the elbow [32] is out of scope for this publication.

The via-point method allows the selection of anatomical origin and insertion nodes for the muscles. As a result the modelled muscle length is almost identical to the anatomical muscle-tendon length. This enables the usage of anatomical data from literature for the muscle parameters. The routing is done in LS-DYNA using sliprings fixed to bones in certain positions in space. The routing parameters, i.e. the offset length of the muscle and the position of the via-point, can be chosen independently of the muscle parameters to match the anatomy. This approach makes it possible to model the muscle-tendon dynamics correctly and at the same time making the lever arm vs. joint angle curve and thus the resulting elbow torque more realistic.

In Fig. 13 different strategies for modeling the triceps are shown. A direct line of action approach, lever arms of 10 and 20 mm, and the via-point routing method are shown. In contrast to the other methods, the application of the via-point method can deliver correct lever arms for the complete range of motion of the elbow and fits the experimental corridor from [40] best, see Fig. 14. Additionally, the proposed via-point routing method improves the numerical stability of the model as it provides correct force application directions. Most importantly, the muscle dynamics is independent from the actual length of the combined elements (muscle + seatbelt, Fig. 2) and their path complexity in FE AHBMs. This is because of the separation of muscle and routing parameters in the FE model.

Fig. 13

Different modeling strategies for the triceps at the elbow

Fig. 14

Lever arms for different modeling of the triceps. Colours are identical to Fig. 13. Corridor taken from [40]

In comparison with the LS-DYNA *MAT_156 muscle model a 10 times speed-up is achieved for the validation simulations with one single muscle elements. As this set-up is clearly not very realistic, the ViVA arm simulations are repeated with *MAT_156 muscles to be comparable to the simulations with the extended Hill-type muscles. Here no speed-up is achieved, since most of the CPU time is used for the processing of the volumetric elements and the time needed to process truss elements is insignificant in comparison.

Conclusion and outlook

The upcoming challenges in the field of automotive safety, namely active safety systems and autonomous driving, will require and benefit greatly from AHBMs. The Hill-type muscle model already existing in LS-DYNA has a limited accuracy because it lacks an internal degree of freedom and in addition is difficult to parameterize. Here, an extended Hill-type muscle model was implemented, verified and validated successfully. The source code, parameters and an example set-up for LS-DYNA are provided at https://zenodo.org/record/826209. The verification and validation was done in comparison with experimental data sets from piglet, cat and rat muscles. The results are in very good agreement with the experiments and the new muscle model improves the accuracy available for AHBMs in LS-DYNA considerably. Moreover, the muscle model incorporates the activation dynamics, essential for correct simulations on small time horizons of dynamic active movements. Additionally, a convenient option for routing the muscle around joints was proposed. By introducing an offset to the length of the muscle element, it is possible to route the muscle using e.g. the via-point method, while at the same time the muscle will display the correct dynamics of the full muscle. This also means, that the parameters for the muscle model can be set independently of the routing.

Although the current model allows to predict the gross dynamic contraction characteristics of biological muscles, it has its limitations. For one, it is a force element predicting a scalar force value which is then applied between origin and insertion, or redirected by via-points. Contact forces and resulting shifts in the force direction or their influence on the active muscle force [35] are neglected. Besides that, several physiological effects of the muscle contraction are currently not considered, starting with muscle-morphology specific parameters such as the pennation angle [41] or the fibre composition [42]. Also on the dynamic level, e.g., modelling the force-velocity relation for the eccentric (lengthening) contractions is difficult, as little data is available. Some of the data suggests more complex relations than modelled here for extensive strains [43], which, however, are not reached in our simulations. Furthermore, the experimentally found history effects causing force enhancement and force depression after stretch and shortening are currently not considered, but may be included in more extended approaches [44, 45]. Finally, the muscles model considers no mass or mass distribution, which, however, plays a role in dynamic contractions [46].

To utilize the full potential of the AHBMs, a control strategy for the activation of the muscles is needed. As a controller realization is not in the scope of this work, the authors recommend the review by [47] as a reliable source of information for muscle activations schemes and strategies in AHBMs. In principle, controllers are required which either maintain a desired position against perturbations or allow for the generation of a desired movement. Such controllers can be implemented in the current framework and may be easily added to the code provided in the Appendix.

With this, we provide a comprehensive and valid approach to implement an extended Hill-type muscle model in LS-DYNA, including muscle-tendon properties, biochemical activation dynamics, and muscle routing. By providing the code and the material cards, we hope that this will allow other researchers to work on more biofidelic AHBMs.

Card 1	1	2	3	4	5	6	7	8
Variable	MID	RO	MT	LMC	NHV	IORTH	IBULK	IG
Default	–	1E-6	41	32	15	0	31	32

Variable	Description
MID	Material identification. A unique number or label not exceeding eight characters must be specified
RO^∗†	Mass density. Not used by the material model
MT	User material type. In this case 41 must be defined
LMC	Length of material constant array. For this material 32 must be set
NHV	Number of history variables to be stored. 15 are required for this material
IORTH	EQ.1: if the material is orthotropic
	EQ.2: if material is used with spot weld thinning
	EQ.3: if material is orthotropic and used with spot weld thinning
IBULK	Adress of bulk modulus in material constants array
IG	Adress of shear modulus in material constants array

Card 2	1	2	3	4	5	6	7	8
Variable	IVECT	IFAIL	ITHERM	IHYPER	IEOS	LMCA
Default	0	1	-	–	–	–	–

Variable	Description
IVECT	Vectorization flag (on = 1). A vectorized user subroutine must be supplied
IFAIL	Failure flag
	EQ.0: No failure
	EQ.1: Allows failure of shell and solid elements
	LT.0: |IFAIL| is the address of NUMINT in the material constants array
ITHERM	Temperature flag (on = 1). Compute element temperature
IHYPER	Deformation gradient flag (on = 1 or −1, or 3). Compute deformation gradient, see LSTC Appendix A: LS-DYNA material cards for validation simulations
IEOS	Equation of state (on = 1)
LMCA	Length of additional material constant array

Card 3	1	2	3	4	5	6	7	8
Variable	Act	STIM	q0	tauq/c	bq/eta	k	m	lOffset
Default	–	–	–	–	–	–	–	0

Variable	Symbol	Description
Act		EQ.0.0: Input of activation values
		EQ.1.0: Calculation of activation with Zajac depending on stimulation
		EQ.2.0: Calculation of activation with Hatze depending on stimulation
STIM	STIM	LT.0.0: Constant stimulation or activation level. Depending on Act
		GT.0.0: LCID specifing the stimulation or activation
q0	q_0	Minimum value of activation q
tauq † / c	τ_q / c	If ACT.EQ.1.0: time constant of rising activation
		If ACT.EQ.2.0: Hatze constant c
bq / eta	β_q/ η	If ACT.EQ.1.0: ratio between τ_q and time constant of falling activation
		If ACT.EQ.2.0: Hatze constant η
k	k	If ACT.EQ.2.0: Hatze constant k
m	m	If ACT.EQ.2.0: Hatze constant m
lOffset †	l_Offset	Muscle length offset added to beam length before calculation of the muscle. l_MTC,i=l_Beam,i+l_Offset

Card 4	1	2	3	4	5	6	7	8
Variable	Fmax	lCEopt	dWdes	nCEd	dWasc	nCEa	Arel0	Brel0
Default	–	–	–	–	–	–	–	–

Variable	Symbol	Description
Fmax †	F_max	Maximum isometric force
lCEopt †	l_CE,opt	Optimal fibre length
dWdes	ΔW_des	Width of F_isom(l_CE) on descending limb
nCEd	ν_CE,des	Exponent of F_isom(l_CE) on descending limb
dWasc	ΔW_asc	Width of F_isom(l_CE) on ascending limb
nCEa	ν_CE,asc	Exponent of F_isom(l_CE) on ascending limb
Arel0	A_rel,0	Maximum value of A_rel
Brel0 †	B_rel,0	Maximum value of B_rel

Card 5	1	2	3	4	5	6	7	8
Variable	Secc	Fecc	LPEE0	nPEE	FPEE	lSEE0	dUSnll	dUSl
Default	–	–	–	–	–	–	–	–

Variable	Symbol	Description
Secc	S_e	Step in inclination of F_CE(l˙_CE=0) between eccentric and concentric force-velocity relations
Fecc	F_e	Coordinate of pole in l_CE(F_CE) normalised to F_maxqF_isom(l_CE) for l_CE>0
LPEE0	L_PEE,0	Rest length of PEE normalised to l_CE,opt
nPEE	ν_PEE	Exponent of F_PEE(l_CE)
FPEE	F_PEE	Force of PEE if l_CE is stretched to ΔW_des
lSEE0 †	l_SEE,0	Rest length of SEE
USnll	ΔU_SEE,nll	Relative stretch at non-linear-linear transition in F_SEE(l_SEE)
duSl	ΔU_SEE,l	Relative stretch in linear part for force increase ΔF_SEE,0

Card 6	1	2	3	4	5	6	7	8
Variable	dFSEE0	Damp	Damp1	Damp2	Output	dtOut	IBULK^∗	IG^∗
Default	-	3	-	-	0	0	-	-

Variable	Symbol	Description
dFSEE0 †	ΔF_SEE,0	Force at non-linear-linear transition in F_SEE(l_SEE)
Damp		EQ.1.0: parallel damping
		EQ.2.0: Serial damping
		EQ.3.0: Serial force dependent damping
		Else: No damping
Damp1	d_PE†	If Damp.EQ.1.0: damping coefficient of PE
	d_SE†	If Damp.EQ.2.0: damping coefficient of SE
	D_SDE	If Damp.EQ.3.0: dimensionless factor to scale d_SE,max
Damp2	R_SDE	If Damp.EQ.3.0: minimum value of d_SE normalised to d_SE,max
Output		Definition of desired output content of outputfile fort.(idele) for each muscle element:
		EQ.0. no outputfile
		EQ.1. basic output (idele, tt, ncycle, q)
		EQ.2. basic output and force data
		(idele...,F_MTC,F_SEE,F_SDE,F_CE,F_PEE,F_isom)
		EQ.-1. basic output and length data
		(idele...,l_MTC,l_CE,l˙_MTC,l˙_CE)
		EQ.-2. basic output and force and length data
		(idele...,F_MTC...,l_MTC...)
dtout †		Timestep of outputile fort.(idele)
IBULK^∗		Bulk modulus. Needed by LS-DYNA to calculate time-step automatically
IG^∗		Shear modulus. Needed by LS-DYNA to calculate time-step automatically