Fiber-based modeling and simulation of skeletal muscles

This paper presents a novel fiber-based muscle model for the forward dynamics of the musculoskeletal system. While bones are represented by rigid bodies, the muscles are taken into account by means of one-dimensional cables that obey the laws of continuum mechanics. In contrast to standard force elements such as the Hill-type muscle model, this approach is close to the real physiology and also avoids the issue of wobbling masses. On the other hand, the computational cost is rather low in comparison with full 3D continuum mechanics simulations. The cable model includes sliding contact between individual fibers as well as between fibers and bones. For the discretization, cubic finite elements are employed in combination with implicit time stepping. Several validation studies and the simulation of a motion scenario for the upper limb demonstrate the potential of the fiber-based approach.


Introduction
The simulation of the musculoskeletal system is of great interest in a variety of fields, e.g., in the study of human movements, the design of effective treatments after injuries and the design of ergonomic workplaces. Consequently, strong efforts are currently being made towards the development of realistic digital human models, and the first corresponding software packages are available now. Here we mention among others the commercial packages AnyBody 1 and Biomechanics of Bodies, 2 and the open-source packages OpenSim 3 and MSMS. 4 In such models, the skeleton is typically represented by a rigid multibody system. This is a widely accepted assumption due to the large differences in stiffness of hard bone tissue and soft muscle tissue. Concerning the modeling of the muscles, the state-of-the-art is the use of muscle models that date back to the early work of Hill [8]. These models are computationally cheap but may lack realism in simulations. In particular, the mass is not accounted for, which led to the introduction of so-called wobbling masses [5]. Detailed three-dimensional muscle models from continuum mechanics, on the other hand, are prohibitively expensive in a multibody context.
For these reasons, we introduce here a new muscle model that features mass, shape and elasticity and which is at the same time computationally attractive. The key idea is based on a one-dimensional cable structure that incorporates large deformation and thickness change. Such a fiber-based model may represent a complete muscle or it can be bundled with other fibers in order to mimic the real composition of muscles. Its spatial discretization requires much fewer degrees of freedom than a three-dimensional approach, and in this way the fiber-based model balances the need for more realism with computational efficiency.
Before further proceeding in this direction, we summarize the state-of-the-art in digital human modeling and related fields. Looking at Hill-type muscles in more detail, we observe that the muscle paths are often assumed to be the straight line between the muscle origin and the insertion point. In order to account for a more realistic geometry, the muscle paths can be enhanced by so-called via points [1,3,10,19]. Another possibility to enhance the fidelity is to consider wrapping surfaces. Then the surfaces of the bones and other structures are additionally modeled and define constraints for the muscle paths. This approach leads to a shortest distance problem between the origin and the insertion point considering the wrapping surfaces as constraints [4,13,22]. Mechanically, Hill-type approaches represent line of action models where the material properties of the muscle-tendon complex are represented by lumped parameter systems [33]. For further studies on Hill-type muscles, we refer to [6,14,15,27,30,31] High fidelity models of muscles are based on nonlinear continuum mechanics. They are able to represent three-dimensional geometry, spatially varying features, the multiscale architecture of muscles, and multiphysics effects [2,18,26]. However, the increased modeling complexity and computational cost pose restrictions on the usability of continuum mechanics models. Therefore, many research works on continuum-mechanical models focus on specific aspects in single muscles, e.g., the mechanical behavior of tendon tissue [12], the effect of micromechanical features on the overall mechanics [24,25], or the extension of continuum-mechanical models to include chemo-electro-physiological aspects [7,17]. The coupling of microvibrations and macroscale deformation, which is an emerging topic in today's orthopedic therapies, is discussed in [29]. While investigating these effects in isolation is important, the ultimate goal is to conduct high fidelity simulations of the musculoskeletal system as a whole.
In the present paper, we model a muscle as a three-dimensional continuum located around a one-dimensional curve in space. Starting from this geometric setting, we derive a onedimensional cable model that incorporates large deformation and thickness change. The total stress tensor is additively decomposed into a passive, an active and a prestress contribution. In order to model the passive response, an incompressible Kelvin-Voigt material law for finite strains is used. The active stress contribution is modeled by the relations given in [18]. Considering the coupling of the skeleton and the muscles at the insertion points, we use classical rigid body coordinates for the bones and absolute displacement coordinates for the cable degrees of freedom. Thus there is no floating reference, as opposed to the typical approaches in flexible multibody dynamics [23,28]. Instead, the large deformation model includes the rigid body modes directly. Furthermore, the contact dynamics on the lateral surface of the muscles and the rigid bodies are taken into account.
The paper is organized as follows. Based on nonlinear 3D continuum mechanics, we derive a cable model incorporating large deformations, incompressible hyper-elastic and viscous material response, prestressing, active stresses and thickness change of the crosssection in Sect. 2. In addition, the contact between two fibers is formulated by means of a penalty approach. In Sect. 3 we introduce a flexible multibody system framework and describe the coupling between the fibers and the rigid bodies. The discretization in space and a subsequent discretization in time of the governing equations are performed in Sect. 4. Furthermore, in this section we discuss the numerical treatment of the contact mechanics and present a safeguard algorithm in order to solve the nonlinear system of equations in an robust way.

Fiber-based model of skeletal muscle
In this section, we derive the fiber-based model, starting with the standard nonlinear equations of motion of a visco-elastic solid. To that end, we introduce additional kinematic and kinetic assumptions according to a cable model. In the next step, the weak formulation of the governing equations is stated. Finally, the contact between two fibers is formulated by means of a penalty approach.

Fundamentals of continuum mechanics
We model the homogenized muscle-tendon complex by means of continuum mechanics. Thus, the motion of a muscle is governed by the nonlinear elasticity problem. Let V 0 ⊂ R 3 be the initial configuration occupied by a muscle at the beginning of the simulation. We are interested in the deformation at all times t ∈ [0, T ], where T > 0 is the simulation end time. For a material point X ∈ V 0 , its position x(X, t) in the deformed configuration V t is expressed in terms of a displacement vector field u : With these definitions, the deformation gradient F : where ∇ X is the gradient operator with respect to the initial configuration V 0 and I the second order identity tensor. In the following we also use the right Cauchy-Green deformation The motion of the muscle is governed by the linear balance of momentum, which reads in Lagrangian form Here, S : V 0 × [0, T ] → R 3 ⊗ R 3 is the second Piola-Kirchhoff stress tensor, b 0 : V 0 → R 3 the body force acting on the volume of the muscle, and ρ 0 : V 0 → R the mass density. Furthermore,ü is the second time derivative of the displacement field. We assume that the second Piola-Kirchhoff stress tensor can be additively decomposed into a passive, an active, and a prestress part, In the present paper, the passive contribution is determined by an incompressible viscoelastic Kelvin-Voigt material response [32] between muscle (M) and tendon (T) material. The material parameters μ i ∈ R and η i ∈ R with i = M, T are the shear modulus and viscosity, respectively. Furthermore, p : V 0 → R is the volumetric stress response, which is not determined by the material law, due to the assumption of an incompressible material. For the active part, we follow [18] and use where α : [0, T ] → R is the activation parameter, λ f : V 0 × [0, T ] → R is the fiber stretch, S max ∈ R is the maximal stress which can be produced at the optimal fiber stretch λ opt f ∈ R, M ∈ V 0 → R 3 ⊗ R 3 is defined by the dyadic product of the fiber direction in the initial configuration [9], and the parameters ν ∈ R, W ∈ R influence the shape of the active response function. Finally, the prestress S pre in (4) is assumed to be self-equilibrated, i.e., Div(F · S pre ) = 0 in the initial configuration V 0 . In order to set up an initial-boundary value problem, the field equations above are supplemented by homogeneous initial conditions and boundary conditions on 0 = ∂V 0 . In the present paper we assume that the muscles are at rest at simulation start t = 0, i.e., the initial conditions read The displacement boundary conditions for a single muscle fiber are Here,û : D 0 × [0, T ] → R 3 is a given displacement field on the Dirichlet boundary D 0 . In order to state the Neumann boundary conditions, we introduce the Cauchy traction t : Fiber-based modeling and simulation of skeletal muscles with n 0 : 0 → R 3 being the unit exterior vector normal to the boundary surface 0 . Then, the traction boundary condition states wheret : N 0 × [0, T ] → R 3 is the given traction vector on the Neumann boundary N 0 .

Cable model for muscle fiber
In this section a cable model for the continuum mechanics setting introduced above is developed. To this end, we introduce the parameter domain P = [0, 1] × [0, 1] × [0, 2π] with the cable radius r and assume that the initial configuration V 0 of one muscle is parametrized by X 3D : P → V 0 with where X(θ 1 ) : [0, 1] → 0 is a parametrization of the initial centerline curve 0 , N : [0, 1] → R 3 and B : [0, 1] → R 3 are the normal vector and the binormal vector to 0 defined by and r 0 : [0, 1] → R is the initial radius of the cross-section. Furthermore, X ,θ 1 and X ,θ 1 θ 1 denote the first and second derivative with respect to θ 1 . The shape of the cross-section is described by φ 1 , φ 2 : [0, 1] × [0, 2π] → R. In the present paper we consider only circular cross-sections by setting The current configuration V t is parametrized by x 3D : P → V t as where x(θ 1 ) : [0, 1] → is a parametrization of the deformed centerline curve , and r : [0, 1] → R is the current radius of the cross-section. The normal vector n(θ 1 ) and the binormal vector b(θ 1 ) are defined analogously to (12). Thus, the base vectors G i = ∂X 3D ∂θ i with i = 1, 2, 3 are given for the initial configuration by with μ 0 = 1 − θ 2 κ 0 sin θ 3 , where κ 0 = A 1,θ 1 · N is the curvature, and ϕ 0 = B ,θ 1 · N the torsion. Furthermore, A 1 = X ,θ 1 is the tangent vector to the centerline curve in the initial configuration. Therefore, the components of the metric tensor Due to (14), the base vectors in the current configuration are with μ = 1 − θ 2 κ sin θ 3 , κ = a 1,θ 1 · n, and ϕ = b ,θ 1 · n. Furthermore, a 1 = x ,θ 1 is the tangent vector to the centerline in the current configuration. Thus, the components of the metric are given by ⎡ ⎢ ⎢ ⎢ ⎣ g 11 g 12 g 13 The volume elements are given by the determinants of the respective metric coefficient matrices, Next, we formulate the incompressibility constraint on the cross-section level. This reads 2π Inserting (18) and (19), we immediately obtain the relation This relation can be interpreted as follows. The expression ||a 1 || ||A 1 || is the ratio of the lengths of the tangent vectors to the centerline curves in the initial and the current configurations and therefore represents the stretch of the centerline curve. In case of an elongation of the centerline curve, the cross-section is supposed to shrink. Therefore, the inverse of the stretch appears in (21). Due to circular cross-sections, the area scales quadratically with the radius. This leads to the square root in (21). Furthermore, we note that, once the tangent vector in the current configuration a 1 is known, the radius of the circular cross-section in the current configuration can be computed.

Passive stress response
In order to proceed with the passive stress response in local coordinates, we introduce some more notation. The dual basis vectors to G i are denoted by G j and fulfill Furthermore, the contravariant coefficients of the metric are G ij = G i · G j . The right Cauchy-Green deformation tensor is given in local coordinates as C = In the present model, we assume that the radii r 0 (θ 1 ) as well as r(θ 1 ) do not vary rapidly along the muscle fiber such that r 0,θ 1 = r ,θ 1 ≈ 0. Furthermore, bending and torsion is neglected. Thus, the right Cauchy-Green deformation tensor is approximated by and the inverse is given by The time derivative of the right Cauchy-Green deformation tensor is given bẏ due to The passive contribution of the second Piola-Kirchhoff stress tensor is given by with Due to the slenderness of a muscle fiber, it is natural to introduce the common assumption in cable models of vanishing stress in radial direction, i.e., S p,22 ≈ 0. The enforcement of this assumption allows us to statically condensate the volumetric response Due to (21), substitution of (29) into (28) yields Active stress We assume that the active stress is generated such that it acts only along the direction of the centerline curve. Therefore, we have with where λ opt ∈ R is the optimal centerline stretch and the stretch of the centerline λ : → R is defined by Prestress For the prestress contribution, we have where σ 0 ∈ R is an input parameter. Thus, in total the second Piola-Kirchhoff stress tensor is given by with

Weak formulation
The weak formulation of the balance of momentum (3) reads as follows: is fulfilled for all appropriate test functions ξ ξ ξ : V 0 → R 3 . Here, V 0 is the cable volume in the undeformed configuration, whereas is the boundary surface in the current configuration and consists of the cross-sections at the cable ends and the lateral surface. For the weak formulation of the cable model, we consider only test functions which are constant over the cross-section. This allows us to preintegrate in (37) over the cable cross-section in the undeformed configuration and obtain where S 11 is given by (36) and now ξ ξ ξ : 0 → R 3 . In the following we will present our treatment of the integral over the boundary , leading to coupling and contact terms. To this end, we consider the decomposition On the free boundary f ree , we have the Neumann boundary conditiont = 0, i.e., we do not have to consider the integration over f ree further. The coupling of the cable end cross-sections c and a rigid body system will be presented in Sect. 3.2. The treatment of the contact surfaces R−F between the lateral cable surfaces with the rigid body system is described in Sect. 3.3, whereas the contact surfaces F −F between different lateral fiber surfaces are treated in Sect. 2.4.

Contact between fibers
In this section we discuss the contact between two muscles with centerlines 1 and 2 . We remark that, for simplicity, only a single muscle was considered in the previous sections. To consider multiple fibers, the weak form of all individual fibers have to be summed up. When the muscles are treated as three-dimensional continua, the two-dimensional fiberfiber contact surface is given by F −F = 1 ∩ 2 . In order to simplify the situation, we use a curve-curve contact that considers the cross-sections. Furthermore, instead of the classical master-slave approach, we use the so-called two-half-pass approach first presented in [20,21]. Therefore, in the weak form related to the fiber with centerline 1 , we use the approximation where the integration over the contact surface F −F 1 is replaced by a part of the centerline F −F 1 ⊂ 1 and the surface traction t is replaced by the line forcet : MM 1 → R 3 . In order to specify F −F 1 andt, we define the dimensionless gap function g : 1 → R, where d : 1 × 2 → R + is the distance of a point x ∈ 1 to the closest point y ∈ 2 on the second fiber involved in the contact, The radii of the muscles in the current configuration at x ∈ 1 and y ∈ 2 are denoted by r x ∈ R and r y ∈ R, respectively. In view of the kinematics in (14), we note that the development of a contact zone between two cables is not possible. Thus, we allow that the muscle fibers penetrate each other by using an approach which can be interpreted as a penalty method. To this end, the contact curve is given by F −F 1 = {x ∈ 1 | g(x) < 0} ⊂ 1 and the contact forcet is determined byt (42) Here, the direction of the contact force is given by whereas the magnitude is given by with the penalty parameter κ ∈ R which has the dimension of force per length, and the dimensionless parameter γ ∈ R determining the maximal possible penetration. In (44) we have chosen a rational functional relation for the contact force in dependence of the gap function. It has the desired property that F → ∞ for 1 + γ g → 0, e.g., the contact force becomes infinitely large if the penetration approaches a certain value. This excludes the nonphysical behavior that two fibers could cross each other, which could occur when using a linear functional relation for the contact traction. Due to the nominator in (44) configurations with too large penetration, e.g., when d < r x + r y 1 − 1 γ , are not admissible, see Fig. 1.

Coupling the cables with a rigid multibody system
In this section we describe our strategies to couple the cable model for muscles with a rigid body model representing bones. We start with the introduction of the rigid body system. Afterwards, the coupling at the muscle insertion points is considered. Finally, the contact along the muscles and the rigid bodies is treated. Though the coupled system represents a flexible multibody system, our model does not use the wide-spread floating reference frame approach. Instead, minimal coordinates for the rigid body motion are combined with the large deformation of the cable model, which includes rotation and translation for the elastic fiber, cf. the absolute nodal coordinates, Shabana [23], and the approach taken by Lang et al. [11] for geometrically exact Cosserat rods.

Rigid body system
Let the state of the rigid-body system be described by the time-dependent configuration vector q(t) ∈ R n using minimal coordinates. The Lagrangian related to the rigid-body system There is no floating reference frame required. The initial configuration (black) of the system is described by (q 1 (0), q 2 (0)) and X(θ 1 ), whereas the current configuration (red) at time t is described by (q 1 (t), q 2 (t)) and is then given by It is composed of the kinetic energy T RB (q) : R n → R and the potential energy V RB (q) : R n → R. We introduce the action Then, the Lagrange-d'Alembert principle states where F is an external force acting on the system. Thus, the forced Euler-Lagrange equations are In the present paper we assume that the geometry of each rigid body is a right circular cylinder.

Coupling between fibers and rigid bodies
We assume a pointwise coupling of the end cross-sections c at the end points of the centerline with the rigid body model, Fig. 2. Thus, the actual shapes of the end cross-sections are not taken into account. Let X c ∈ c 0 be a coupling point in the initial configuration at t = 0, i.e., a material point of a muscle fiber which is attached to the rigid body system. The coupling conditions at X c areĉ RB (X c , q(t)) = X c + u(X c , t), which have to hold for all t ∈ [0, T ]. Here,ĉ RB : R 3 × R n → R 3 is the nonlinear position function of the rigid body system yielding the current position of X c for the configuration q(t). The right-hand side of (49a) is the current position of X c computed by the displacement u(X c , t) of the muscle fiber. Thus, equation (49a) ensures the compatibility of the deformation. In Sect. 4 we will use c =ĉ RB − X c and reformulate (49a) appropriately as a constraint on the end positions of the muscles. Equation (49b) is the Newton's second law, where F RB and F M are the respective forces at the coupling point. It will be considered by means of Lagrange multipliers in Sect. 4.

Contact between rigid bodies and muscles
Similarly to Sect. 2.4 we consider a cylinder-curve contact, where we take the cross-section of the curve and the cylinder into account. To this end we approximate the integral over the contact surface R−F by The gap function is now defined as where d is the distance of a point x ∈ R−F to the closest point y = y 0 − ((x − y 0 ) · t RB )t RB on the cylinder axis. Here, y 0 ∈ R 3 is a point on the cylinder axis and t RB ∈ R 3 is the unit vector in the direction of the axis. In (51) r x ∈ R denotes the radius of the muscle crosssection at x, whereas r RB ∈ R is the radius of the cylinder involved in the contact. The contact tractiont is computed by means of (42).

Numerical realization
In this section, we outline the numerical discretization of the coupled model with rigid bodies and elastic muscle fibers.

Semidiscretization in space
We employ a semidiscretization in space of the cable model by means of cubic Hermite basis functions. Thus, the cable middle curve and the displacement unknowns are approximated by C 1 -continuous functions. On one element we have the decomposition where X i describe the initial geometry of the element, u i are the local unknowns, and φ i (ξ ) are the cubic Hermite element shape functions The element shape functions are pieced together to basis functions N i , i = 1, . . . , N s , and result in the usual way in a semidiscrete substructure for each fiber.
In the next step, the semidiscretized individual muscle fibers with unknown displacement variables u h (t) need to be coupled with the rigid body dynamics q(t) by means of the constraints (49a) and corresponding Lagrange multipliers, which incorporate (49b). For simplicity, we restrict the discussion of the coupled system to a single fiber, for which the resulting equations read The Euler-Lagrange equations (48) for the rigid body system are reformulated in (54a). In (54b) the entries of mass matrix M M of a single fiber are given by whereas the inner force vector is given by and the contribution due to body forces is The coupling conditions (49a) are rewritten in (54c) with a Boolean matrix B that extracts those boundary nodes from u h that are subject to the corresponding constraints. Overall, the semidiscretized system (54a)-(54c) forms a differential-algebraic equation (DAE) with additional Lagrange multipliers λ λ λ and Jacobian matrix C(q) = ∂c(q)/∂q. We remark that all integrals in Eqs. (55)-(57) are approximated by Gaussian quadrature.
The constraint (54c) allows us to directly eliminate the superfluous boundary nodes in the coupling interface. As a next step, we sketch this procedure that transforms the DAE (54a)-(54c) to a state space model. For this purpose, we partition the elastic displacements as u h = (u I , u D ) with independent (or interior, respectively) variables u I and dependent (or Dirichlet boundary, respectively) variables u D . Thus, B = (0, I u D ) with an identity matrix I u D of the size of the boundary nodes, and the nullspace matrix of the constraint (54c) is then In other words, it holds that (C(q), −B) · N(q) = (C(q), 0, −I u D ) · N(q) ≡ 0 for all states q. From the constraint (54c), we deduce the relations u D = c(q),u D = C(q)q,ü D = C(q)q + C q (q)(q,q).
By inserting these relations into the DAE (54a)-(54c) and premultiplying the dynamic equations by N(q) T , we can eliminate the dependent variables and the Lagrange multipliers. To this end, the constant finite element mass matrix M M is partitioned into according to the structure of u h . The state space form then reads The state space model requires the derivative term C q (q)(q,q) from the coupling condition. In our approach we find it more straightforward to use a formulation where the coupling condition is still explicitly given while the Lagrange multipliers are eliminated. This is possible due to the special partitioning of the displacements into u h = (u I , u D ). We introduce velocity variables p =q and v h =u h and conclude from (54b) and (60) that Then we insert this expression into (54a) and obtain, along with the differentiated coupling constraint for the velocities, the DAE formulation where . Note that this formulation is completely equivalent to the original DAE (54a)-(54c) and the state space form (61a), (61b). It has a "GGL-like" structure, and by one-time differentiation of the constraints (63c) and (63d), an ordinary differential equation in all variables is obtained. The index of (63a)-(63d) is thus equal to 1. For a general constrained mechanical system, such a partitioning of the unknowns and also the special structure of the constraint equations do not hold, and consequently (63a)-(63d) cannot be derived.

Semidiscretization in time
We consider a uniform discretization of the time interval [0, T ] : (t 0 , . . . , t Nt ), with t i = i t, i = 0, . . . , N t , where t > 0 is the time-step size and N t the number of time-steps. We denote by u i h the discrete displacement field at time-step t i . This notation applies to all discretized quantities. We use the backward Euler method for time discretization. Then, the fully discrete schema reads: with initial conditions q 0 = q 0 , The nonlinear problem (64) is solved by Newton's method in combination with a safeguard algorithm described in Sect. 4.4. The expressions for the tangent matrix are summarized in the Appendix.

Contact numerics
Here, we describe the numerical treatment of the contact between muscles. As all integrals, the last two integrals in (56) related to contact are approximated by quadrature. Therefore, we have where x k ∈ R 3 and ω k ∈ R are the n Q quadrature nodes and weights. In order to evaluate the contact tractiont between two fibers, we have to find the closest point y(ξ ) on each possible contact partner corresponding to the quadrature point x k . The necessary condition for y(ξ ) to be a minimizer reads In order to obtain y(ξ ) for a given point x k , we remark that the centerline is discretized by Hermite shape functions (53), i.e., by piecewise cubic polynomials. Thus, we have to solve for each possible contact element a quintic equation in ξ of the form with v 1 = y 1 , and the local coefficients y i = Y i + u y i describing the current position of the respective finite element. In order to have a robust method we would like to find all roots of (67). To this end, we use the Matlab function roots. After having found all roots we evaluate also the distance at the element borders and take the point with minimal distance as y(ξ ).

Safeguard algorithm
In each time-step the resulting nonlinear system (64) is solved using Newton's method. However, in order to have a robust method we have to ensure that the contact penetrations do not become too large and not admissible. To this end, we employ a safeguard algorithm and scale the step-size of the Newton method by a factor δ ≤ 1 if necessary. Let the Newton iteration be where u i and u i+1 is the current and next iterate and u the update. In view of (44), the contact traction becomes infinite whenever at any spatial location x. Thus, we have to ensure that (δ) > 0 for the contact between two fibers. However, the evaluation of the function (δ) requires the full evaluation of the contact algorithm. In order to avoid this costly evaluation, let e (δ) be an estimate of (δ) such that (δ) ≥ e (δ). Then, we want to determine δ ∈ [0, 1] such that with chosen parameter 0 <δ < 1. In our implementation we compute the estimate elementwise. Thus, for an element E we have with Here, min E and max E is the respective minimum and maximum on the considered element E, whereas min \ E and max \ E is the respective minimum and maximum on all other fibers except the fiber containing element E.

Numerical results
In this section we present the results of numerical computations.

Verification of the cable model
In this first example we verify our implementation of the cable model. To this end, we consider a single fiber with properties given in Table 1. The muscle ratio γ M follows the function where θ ∈ [0, 1] is the dimensionless position in the fiber. Following the idea of a manufactured solution, we prescribe the motion of the fiber centerline as For this chosen motion we compute the necessary body force by For the actual computation of (76), we used the symbolic computation capabilities of Matlab. The resulting function is prescribed in our numerical calculation. We have conducted a convergence study of the error ||u h 2 (0.5, t) − u 2 (0.5, t)|| L 2 ([0,T ]) using computations with different spacial and temporal discretizations. Here, denotes the second component of the exact solution (75) in the middle of the cable, whereas u h 2 (0.5, t) is the corresponding numerical solution. The results of the convergence study are illustrated in Fig. 3. We observe the optimal linear convergence rate of the temporal discretization provided the spatial discretization is fine enough.

Verification of the contact integration
In this example we investigate the influence of the contact integration described in Sect. 4.3. To this end we consider the contact between a rigid cylinder and a homogeneous fiber. The model parameters are given in Table 2. We load the fiber by a distributed line load of b z = −1000 N/m. The initial and the deformed configurations are visualized in Figs. 4a and 4b, respectively. We study the effect of the integration and the number of elements on the displacement solution. To this end we observed the minimal vertical displacement (see Fig. 5). As a reference value we have taken the displacement obtained with the most quadrature points and 128 elements. We observe convergence with growing number of quadrature points for a fixed number of elements used. The error due to integration is less than 1% when 4 or more quadrature points are used. Furthermore, also the convergence with respect to the number of elements can be observed. When 32 elements or more are used, the error is below 1%. Therefore, we will use 5 quadrature points per element for the contact integration in the following examples.

Spatial convergence of the contact between fibers
In this verification example we study the spatial convergence of the contact algorithm for the fiber-fiber contact. To this end, we consider the contact between two homogeneous fibers with material parameters given in Table 3. The initial configuration, where no contact is present, is depicted in Fig. 6a. The upper fiber is statically loaded by a uniform line load b z = −5000 N/m in vertical direction such that it deflects and comes in contact with the lower fiber. The deformed configuration after loading is depicted in Fig. 6b. We have solved the problem with different space discretizations and evaluated the displacement in the z-direction u h z (0.5) at the middle of the upper fiber. The results of the convergence study are illustrated in Fig. 7. As a reference value u ref z (0.5), we have taken the result obtained for the finest discretization. We observe the convergence of the error with growing number of elements used.

Upper limb model
We consider a strongly simplified model of the upper limb, consisting of three rigid bodies and five fibers modeling the biceps (two fibers) and triceps muscles (three fibers), see Fig. 8. The rigid body representing the upper arm is assumed to be fixed in space. The second rigid body representing the forearm is assumed to have only one rotational degree of freedom around the y-axis. The third rigid body is an additional wrapping surface in order to force the triceps fibers to bend around the elbow. The biceps and the triceps are modeled by two and three fibers lying (almost) in parallel, respectively. The radii of all muscles and the where θ ∈ [0, 1] is the dimensionless position along the fiber. The remaining model parameters are given in Table 4. We remark that the model parameters have been adapted from [18] and [13]. However, due to the lack of data, the parameters η and σ 0 , as well as the functions in (77), are an ad hoc choice of the authors. The prestress σ 0 has been adapted such that forearm is nearly facing downwards in the initial configuration. Therefore, the results have to be understood in a qualitative way, showing the potentials of the proposed approach. We perform a dynamic forward simulation with the prescribed activation curves α bi (t) = 0.0704 tanh(πt) + 0.0608 sin(πt), for the biceps and triceps fibers, respectively, see Fig. 9. The simulation end time is chosen to be T = 8 s. We have solved the problem with different discretizations in space and time. The discretization parameters and the elapsed computation time are given in Table 5. Here, N e refers to the number of elements per fiber and N t is the number of time-steps. Thus, we have used 40, 80, 160, 320 finite elements resulting in 542, 1022, 1982, 3902 degrees of freedom in the semidiscretized systems. All computations are performed on a personal computer with an AMD Ryzen 7 3700X 8-Core processor, which uses an Ubuntu 18.04 operating system.
The angle of rotation of the forearm is plotted over time in Fig. 10. Since the resulting curves virtually agree we conclude that the discretization error is low. We observe some oscillations in the beginning of the simulation, which are related to the initial conditions.  Due to physical and numerical damping, these oscillations disappear after the first cycle. Figure 11 illustrates the forces in one biceps and one triceps fiber at the coupling point to the upper arm are plotted over time. For the same fibers, the internal forces are displayed over the fiber length for two times in Fig. 12. In Fig. 13 some snapshots in time of the upper limb model are plotted.

Conclusion
A new fiber-based simulation framework for the forward-dynamics of the musculoskeletal system has been presented. All components including the muscles are represented by three-dimensional bodies. Due to their stiffness, bones are modeled as rigid bodies, whereas muscles are described by one-dimensional cables that are derived from continuum mechanics. The advantage of this approach lies in the relatively low computational cost compared to models accounting for full three-dimensional kinematics, without introducing too much assumptions like in lumped parameter models.
In the present approach the rigid bodies are restricted to be cylinders. In future work we plan to incorporate triangulated bone surfaces into our model. Furthermore, we would like to couple the present model with an electrochemical model on the microscopic level that captures the relevant effects of the actin-myosin binding process during muscle contraction [16]. So far, our focus has been on forward-dynamics, like for the simulation of the upper limb model. For realistic applications, however, the control of the biomechanical system  becomes an important challenge. Currently, we are investigating the usage of both classical optimal control and reinforcement learning for this purpose. Finally, so far no attempt of validation against experimental data of the presented approach has been made. This important step will also be part of future work. support by the German Federal Ministry of Education and Research of Germany in the framework of DY-MARA (project number 05M16UKD).
Funding Note Open access funding provided by Graz University of Technology.

Conflict of interest
The authors declare that they have no conflict of interest.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A.3 Linearization of the contact between cables
The contribution to the residuum reads The tangent matrix has the two main parts D ux R f c = θ D uxt · ξ ξ ξ g 1 +t · ξ ξ ξ D ux g 1 dθ, D uy R f c = θ D uyt · ξ ξ ξ g 1 dθ, where u x refers to the displacement of the cable which is integrated and u y refers to the displacement of a second cable involved in the contact. The derivative of the Jacobian is given by The derivatives of the contact traction reads D uxt = F D ux g n + F D ux n, D uyt = F D uy g n + F D uy n, where F is given by (44) and F (g) = 0, g≥ 0, − αβ (β+γ g) 2 , g < 0. (91) The derivative of the gap function g with respect to the displacement u x reads D ux g = D ux d r x + r y − d(D ux r x + D ux r y ) (r x + r y ) 2 , whereas the derivative with respect to the displacement u y is D uy g = D uy d r x + r y − d D uy r y (r x + r y ) 2 .
Now, the derivatives of d can be expressed as Let the change in fiber radius be such that r x (θ ) = (θ )r 0 (θ ) (see (21)). The derivative of the current radius with respect to the displacement u x thus reads D ux r x = − r 0 2 The derivative of the radius r y (ξ ) = (ξ )r 0 (ξ ) with respect to the displacement u x reads D ux r y = ( r 0 + r 0 )D ux ξ, where = − 1 2 The derivative of the radius r y with respect to the displacement u y reads D uy r y = ( r 0 + r 0 )D uy ξ − r 0 2 ||A 1 || ||a 1 || 5 (a 1 · v y ).
The derivatives of the normal vector n are given by Altogether we have D uxt = F −n · v x r x + r y − d(D ux r x + D ux r y ) (r x + r y ) 2 n + F y D ux ξ − v x + n(n · v x ) d , D uyt = F n · v y r x + r y − d D uy r y (r x + r y ) 2 n + F y D ux ξ + v y − n(n · v y ) d . (103)

A.4 Linearization of the contact between rigid bodies and cables
The contribution to the residuum reads R bc (u) = C t(u) · ξ ξ ξ ds x = θt (u) · ξ ξ ξ a 1 dθ.
Thus, the tangent matrix has the form D u R f c = θ D ut · ξ ξ ξ g 1 +t · ξ ξ ξ D u g 1 dθ.
The derivative of the Jacobian is already given in (89). The derivative of the contact force is given by D ut = F D u g n + F D u n.
We define the normal projector P of a cylinder with axis t RB as Then, the normal vector n and the gap function g can be explicitly written as n = P(y 0 − x(θ )) P(y 0 − x(θ )) , where r = r RB + r M 0 and y 0 is a point on the rigid cylinder axis. The derivatives are given by D u n = (I − n ⊗ n) Pv P(y 0 − x) , where the derivative of the radius is analogous to (98).