A Whole-Body Musculoskeletal Model of the Mouse

Neural control of movement cannot be fully understood without careful consideration of interactions between the neural and biomechanical components. Recent advancements in mouse molecular genetics allow for the identification and manipulation of constituent elements underlying the neural control of movement. To complement experimental studies and investigate the mechanisms by which the neural circuitry interacts with the body and the environment, computational studies modeling motor behaviors in mice need to incorporate a model of the mouse musculoskeletal system. Here, we present the first fully articulated musculoskeletal model of the mouse. The mouse skeletal system has been developed from anatomical references and includes the sets of bones in all body compartments, including four limbs, spine, head and tail. Joints between all bones allow for simulation of full 3D mouse kinematics and kinetics. Hindlimb and forelimb musculature has been implemented using Hill-type muscle models. We analyzed the mouse whole-body model and described the moment-arms for different hindlimb and forelimb muscles, the moments applied by these muscles on the joints, and their involvement in limb movements at different limb/body configurations. The model represents a necessary step for the subsequent development of a comprehensive neuro-biomechanical model of freely behaving mice; this will close the loop between the neural control and the physical interactions between the body and the environment.


I. INTRODUCTION
Terrestrial animals exhibit a variety of complex motor behaviors, some rhythmic (e.g., locomotion, grooming) and some non-rhythmic (e.g., reaching and grasping, crouching, posture-control). These motor behaviors result from complex interactions between the neural circuits in the brain and spinal cord, the musculoskeletal system, and the environment (Fig.  1) [1]- [5]. Investigating the underlying mechanisms of these motor behaviors in awake, behaving animals is highly challenging. Computational modeling is a powerful tool to complement experimental studies of neural control. It provides a platform to investigate the underlying neural mechanisms, allowing investigators to reproduce existing experimental observations, propose mechanistic explanations for the observed behaviors, suggest new experiments to test the proposed mechanisms, and propose possible approaches for treatment of different injuries or disorders [6], [7].
Recent advances in mouse molecular genetic approaches have led to substantial progress in the studies of spinal and supra-spinal networks involved in motor control [2], [7]- [15]. Specifically, it is now possible to dissect the neural network by identifying and manipulating specific neural populations and to relate them to specific behavioral outcomes [2]. The ubiquity of molecular genetic tools available for the mouse made it the preferred experimental animal for the study of neural control of movement.
Computational modeling of motor control in freely behaving mice, particularly modeling of locomotion, requires the development of detailed models of both the neural and biomechanical components and, importantly, their interactions during movements [16] (Fig.  1). Such detailed musculoskeletal models have been developed of the human [17]- [21], the guinea fowl [22], the ostrich [23], and the dog [24]. Yet, no detailed whole-body three-dimensional musculoskeletal model of mouse currently exists; and most of the existing neuromechanical models (models incorporating neural controllers with the musculoskeletal system) of mammals have been limited to two-dimensional movements or had significantly simplified musculoskeletal or muscular components of the system [7], [25]- [30].
Here, we present an open-source, configurable whole-body musculoskeletal model of the mouse. We digitized all the bones of the mouse skeletal system and identified the corresponding joints to have a fully articulated rigid body model. Musculature for hindlimbs and forelimbs were modeled as Hill-type muscles. We compared and validated the hindlimb musculature of the model with the previously published single hindlimb model of mouse developed by Charles et al. [31]. The three-dimensional (3D) model allowed us to explore the relationship of muscle moment-arms and moments on the hindlimb and forelimb joints in a comprehensive manner. Our analysis on the operational range of muscle-fibers for limb muscles give insights into the role of active and passive forces under isometric conditions. Finally, we evaluated the sensitivities of parameters to highlight a reduced set of crucial parameters in the complex whole-body mouse model.
In the forelimbs, the shoulder joints were modeled as a spherical joint with three rotational DoFs: retraction-protraction (rotation about the transversal axis), abduction-adduction (rotation about the coronal axis), and external-internal rotation (rotation about the sagittal axis). Elbow joints were modeled with two DoFs: extension-flexion (rotation about the transversal axis) and supination-pronation (rotation about the sagittal axis). Wrist joints were modeled with two DoFs: extension-flexion (rotation about the transversal axis), and abduction-adduction (rotation about the coronal axis). Only a single DoF of flexionextension (rotation about the transversal axis) was modeled between each digit in the forepaws.
In total, the mouse model consists of 225 rotational joints. Table 1 lists the joints in the model.

3) JOINT RANGE-OF-MOTION AND LIMITS-Joint
range-of-motion defines the extent to which a joint can rotate between a minimal and a maximal angle. It is an important attribute for articulated rigid body simulations. For animals, joint ranges are imposed due to a combination of factors such as ligaments, muscles, and elastic forces of tissues. From previous experimental measurements on mice and rats, we have compiled and reported the joint limits for the DoFs of the model in Table 1. In this work, we modeled joint limits to mimic elastic ligaments that engage beyond the specified joint range by applying a resisting torque (modeled as a torsional spring and damper system as described in Opensim v4.2 [34]).

4) ZERO-POSE-Zero
-pose defines the reference coordinate system. In the model, the zero-pose was defined as the pose at which all the joint angles are set to zero with respect to the corresponding coordinate frame. A non-anatomical pose in which some joints go beyond natural joint limits is chosen as the zero-pose to facilitate the coordinate frame transformations. Joint angles are measured with respect to the zero-pose. The model in the zero-pose is shown in Fig. 2A and a possible rest-pose is shown in Fig. 2B.

5) ESTIMATING INERTIAL PROPERTIES-Inertial
parameters are among the most important parameters for an accurate articulated rigid-body physics simulation. These parameters include the mass, the center-of-mass, and the inertia tensor for every bone (and surrounding soft tissues). To compute the inertial parameters, we assumed an uniform density of water (1000 kg/m 3 ) along the body [39]- [43]. The volume around the bones was estimated based on the convex hulls that represent the net volume encapsulated by the skin, soft-tissues, muscles, and the bone [44], [45]. Masses were computed based on the density and volumes of the convex hulls. The center-of-mass and inertia tensors were computed based on the shape of the convex hulls and estimated mass assuming uniform density. The mouse model is 16.90 cm long from the head to the tip of the tail in the zero-pose and weighs 34.32 grams.

C. MUSCLE SYSTEM
Muscles were modeled using the Hill-type formalism [46]. Modeling each muscle required identification of attachment points in the hind and forelimbs and the estimation of their parameters.
1) HILL-TYPE MUSCLE MODEL-Hill-type models make use of passive elements such as springs and dampers and experimental data to represent the active and passive dynamics of a muscle [46]. The contractile element (CE) and the parallel elastic element (PE) represent the muscle fibers; the series elastic element (SE) represents the total series elasticity of the muscle-tendon complex. It is important to keep in mind that due to this definition even muscles with short or no tendon still have non-zero tendon lengths in the parameterization of the Hill-type muscle model. Fig. 3 shows the formulation of the muscle model used in this work. To simplify the dynamics of the muscle model, we assumed rigid tendons [47].
The Hill-type muscle model illustrated in Fig. 3 was characterized by the following parameters: the optimal fiber length (l m o ; the length of muscle-fiber at which the muscle produces maximal active force), the tendon slack length (l t s ; the length of the tendon below which the muscle transfers no force to the attached bones), the maximum muscle-fiber velocity (v m max ), the pennation angle when the fiber length is at its optimal (α o ), and the maximum isometric force (F m 0 ). The formal description of the Hill-type muscle model used [47] is described in Appendix A.
2) MUSCLE ATTACHMENT POINTS-Muscles are attached to the bones via tendons. The attachment at the most proximal bone is called the origin (that bone tends to move less during muscle contraction) and the one at the most distal bone is called the insertion (that bone tends to move most during muscle contraction). The complex muscle paths were approximated as a polyline, a sequence of straight lines starting at the origin and ending with the insertion and are connected by waypoints in between (Fig. 4). The polyline approximation is a common approach to describe the muscle paths [18], [26], [34], [48].
(Refer to Fig. 4 for an example showing the polyline muscle path description around the knee joint). The identification of the muscle attachment points and waypoints for hind and forelimbs used in this study is described below.
a: HINDLIMB ATTACHMENT POINTS: The attachment points for mouse hindlimb muscles have been previously identified by Charles et al. [31], [49]. However, because of the differences in the bone geometry between our model and Charles model [31], it was necessary to transfer the muscle attachment points from Charles model to the bone surfaces of our mouse model. To automate the transfer process and limit the errors, we first performed mesh registration over the two bones from both models (mesh registration involves identifying appropriate landmarks between each bone segment of the hindlimb of our model and theirs). Based on the chosen landmarks, a coordinate transformation matrix was computed to describe an affine transformation from the bone mesh in Charles model to our model. The affine transformation obtained from the mesh registration was then used to transfer the attachment points and waypoints between the models. The abovedescribed process was carried out using the open-source mesh software CloudCompare [50]. Identification of the landmarks on each bone from the two models was done manually. This process was carried out for the pelvis, femur, tibia, and pedal (tarsus, metatarsus, phalange) bones of the mouse hindlimb (Fig. 5A).
b: FORELIMB ATTACHMENT POINTS: Modeling the forelimb muscles was challenging because of the lack of prior studies identifying attachment points of forelimb muscles. Delaurier et al. [51] has previously developed a 3D forelimb atlas of the mouse embryo at embryonic day 14.5 using Optical Projection Tomography and digital segmentation. Although the model was based on data from a mouse embryo, it still provides useful information on muscle attachments to the bones. This data allowed us to identify attachment points and necessary waypoints for 17 forelimb muscles (Fig. 5B). was scaled in proportion to the muscle-tendon lengths (i.e., the ratio of muscle-tendon length measured by Mathewson et al. [54] to the muscle-tendon length in our model for each forelimb muscle was used as the scaling factor for l m o ). The scaling procedure used for hindlimb muscle parameters could not be employed here because there was no previous biomechanical model of the mouse forelimb. F m 0 was computed by multiplying the physiological cross-sectional area (PCSA) of the muscle by isometric stress (σ) taken to be 0.3 N/mm 2 , the same used by Charles et al. [31] for the hindlimb muscles.

3) ESTIMATION OF MUSCLE PROPERTIES-As
The remaining muscle property, l t s , was estimated using an adaptation of the numerical optimization technique formulated by Manal and Buchanan [55]. We removed the restriction that muscle fibers only operate in the ascending region of the force-length curve. We extended the algorithm to consider all the DoFs a muscle operates on while optimizing for l t s . Refer to the Appendix C for further details.

4) COMPUTATION OF MUSCLE-TENDON LENGTH AND MOMENT-ARM-
Muscle-tendon length (l mt ) is the distance from the origin of the muscle on the proximal bone to the insertion point on the distal bone. Based on the polyline approximation for describing the muscle paths (Fig. 4), where P n are the muscle attachment points on the polyline withP 0 being the origin, P N the insertion point and waypoints in between P i : i = 1…(N − 1).
Muscle moment-arms (r) were computed using the perturbation method described in Sherman et al. [56], r = ∂l mt / ∂θ (2) where θ is generalized coordinate representing the DoF of interest across which the muscle moment-arm is to be computed.

D. ANALYSIS OF SENSITIVITY OF JOINT MOMENT-ARMS AND MOMENTS TO CHANGES IN MUSCLE PARAMETERS AND ATTACHMENT POINTS
Like any model, the mouse biomechanical model was developed based on several assumptions and simplifications. To highlight and study the influence of various characteristics of the biomechanical model on the overall motor behavior, we analyzed the sensitivity of joint moment-arms and moments to changes in muscle attachment points and model parameters.
For both analyses, we used the Sobol method [57], a variance-based global sensitivity method. The Sobol method decomposes the proportion of model output variance caused by each individual parameter [58]. This method also allows for the study of the inter-parameter effect on the model's output variance, but we restricted our analysis to first-order indices (main effects). First-order indices "measure the direct contribution of each input factor to the output variance" [58]. A value of 1.0 for the first-order indices indicates that the parameter is solely responsible for all the variance in the models' output, whereas a value of 0.0 represents no influence on the models' output variance. We used the root-meansquare (RMS) value of the joint moment-arm or moment over the DoFs range-of-motion to represent the scalar value necessary to evaluate the sensitivities. The Sobol sensitivity analysis was performed using SALib v1.4 [59], an open source python library for the sensitivity analysis.
The polyline approximation of the muscle paths determines the computation of muscletendon length (l mt ) and thereby the moment-arm (2). For the polyline approximation, change in l mt can be reduced to just two attachment points of the segment that cross the DoF of interest (refer to Fig. 4 showing an example). Therefore, for the analysis of sensitivities to muscle attachments, we referred to the proximal attachment point of the segment as P WO and to the distal one as P WI . The sensitivities to these attachment points were computed by varying their 3D position within a range of ±0.5 mm in x, y, and z directions. With P WO and P WI we had six parameters (x, y, and z coordinates for each point) to study the sensitivity of the joint moment-arm. We limited our analysis to study how the 3D location of the two attachment points influence the joint moment-arms. As the total sum of the first-order indices is equal to 1.0, we reported the sum of the first-order indices along the individual 3D coordinates (x-y-z) for each attachment point. These values represent a direction-independent measure of moment or moment-arm sensitivity to the attachment points. Data on the individual sensitivities are available in the supplementary material.
Muscle force production depends not only on the geometric relationship defined by the muscle-path but also on the muscle dynamics. The Hill muscle dynamics are parameterized by four main parameters (as described in section II-C1).

III. RESULTS
The parameters of all muscles in the model are specified in Table 2. Because of the differences in the available experimental data, muscle parameters for the hindlimb and forelimb muscles were obtained using different approaches. For hindlimb muscles, the maximum isometric force F m 0 and pennation angle (α o ) were directly taken from Charles et al. [31] while tendon slack length (l t s ) and optimal fiber length (l m o ) were scaled using the method described in [53] to account for the differences in model dimensions. For forelimb muscles, F m 0 and l m o were scaled from Mathewson et al. [54] while keeping α o the same. l t s was optimized using a modified version of the algorithm originally described by Manal and Buchanan [55] (see section II-C3.b for details).

A. MOMENT-ARM ANALYSIS FOR THE HINDLIMB MUSCLES
Because of the differences in bone geometries, we transferred the muscle attachment points from the hindlimb model developed by Charles et al. [31], [49] to our model. The mesh registration technique described in section II-C2.a was used to transfer the attachment points. The attachment points of the muscle have a direct influence on the muscle momentarms ( Fig. 4) and consequently on the moments on the joints they influence. To compensate for the lack of wrapping surfaces, additional waypoints were introduced along the muscle path when necessary. To validate the muscle attachment process, moment-arms of two flexor muscles of the hip, knee, and ankle joints from Charles model were compared with the moment-arms in our model (Fig. 6). For this comparison, moment-arms were normalized to their respective thigh (femur) lengths, since the two models are of mice of different age and size. The moment-arms from the two models are in good quantitative and qualitative agreement throughout the range-of-motion for the three hindlimb joints. Fig. 7 presents a comprehensive overview of each muscle's influence over every DoF it spans. We present the functional grouping of the muscles based on the moment-arm and moment. Each DoF was sub-divided into two functions, representing the possible directions in which the muscle can influence the DoF. Thus, for every DoF it spans over, a muscle has the possibility to apply a moment either on one of the DoF functions' or both (zerocrossing). Moment-arms were computed in the default pose shown in Fig. 5; moments were computed assuming maximal muscle activation (a(t) = 1.0). Fig. 7A,B shows the functional grouping of hindlimb muscles based on the moment-arm and moment, respectively. For hip flexion, iliacus, psoas major, psoas minor, and rectus femoris have the maximal moment-arms. However, since rectus femoris has the largest maximum isometric force, it exhibits the highest moment towards hip flexion. For hip extension, adductor longus, adductor brevis, quadratus femoris, biceps femoris posterior, semimembranosus, and semitendinosus have the dominant moment-arms. However, grouping them by moment reveals that quadratus femoris, semimembranosus, and semitendinosus are the dominant muscles. This classification also highlights that all hip flexors except psoas major, psoas minor, iliacus, and rectus femoris have zero-crossings (Fig. 7A). While gluteus maximus muscle has the highest moment-arm to abduct the hip, semimembranosus along with gluteus maximus also strongly contributes to hip abduction moments as their maximum force is larger. For hip adduction, adductor brevis, adductor longus, adductor magnus, gracilis anterior, and gracilis posterior are the dominant adductors when we consider the moment-arms. Considering moments, quadratus femoris and rectus femoris dominate over other muscles because of their large F m 0 . For hip external rotators quadratus femoris has the most significant moment-arm and moment in the group. Also, gluteus maximus has the largest moment and moment-arm among the hip internal rotators.

1) HINDLIMB MUSCLES-
For knee flexion, semitendinosus, biceps femoris posterior, gracilis anterior, gracilis posterior, and semimembranosus have the strongest moment-arms, but semimembranosus and semitendinosus show the largest moments followed by lateral gastrocnemius. For knee extension, rectus femoris, vastus intermedius, vastus lateralis, and vastus medialis have the largest moment-arms, while rectus femoris exerted the largest moment, followed by vastus lateralis and lateral gastrocnemius.
Several muscles show strong moment-arms for ankle plantarflexion with soleus and plantaris having the largest moment-arms. Considering the moments, lateral gastrocnemius overshadows all the other ankle plantarflexors. For ankle dorsiflexion, tibialis anterior exhibits the strongest moment-arm and moment. Among ankle evertors, peroneus digiti quarti and peroneus longus show the largest moment-arm but flexor digitorum longus and peroneus longus, due to their large maximum isometric force, also significantly contribute to the moments. For ankle inversion, medial gastrocnemius has the highest moment-arm and moment among all the muscles in its group. Peroneus digiti quarti, peroneus longus, and peroneus brevis have large moment-arms for ankle abduction, while lateral gastrocnemius, flexor digitorum longus, and peroneus longus Can exert large moments for ankle abduction. Medial gastrocnemius, plantaris, and soleus have significant moment-arms for ankle adduction group, and lateral gastrocnemius and medial gastrocnemius have the highest moments.
2) FORELIMB MUSCLES-Figs 7C,D show the functional grouping of forelimb muscles based on the moment-arm and moment, respectively. Anatomically, the forelimb is more complex than the hindlimb, the scapula is a floating link with all six DoFs and elbow joints have complex joint rotations. The kinematic complexities are reflected in the identification of muscle functions. Most shoulder joint muscles originate from vertebral processes. Here, we have not included those muscles that originate from the spinal segment due to the lack of available experimental data; the model only includes a limited number of proximal muscles in the forelimb. Shoulder retraction is actuated by the long head of biceps brachii whereas shoulder protraction is actuated by coracobrachialis and the lateral head of triceps brachii. Both latter muscles have similar maximum moment-arms, but the maximal isometric force of the lateral head of triceps brachii is much larger compared to coracobrachialis, resulting in a larger moment. Coracobrachialis, the long head of biceps brachii, and the lateral head of triceps brachii contribute similarly to shoulder abduction both in terms of moment-arm and moments, while the long head of biceps brachii and the lateral head of triceps brachii contribute to shoulder adduction. This shows that shoulder adduction is possible in the model only due to the zero-crossing of the long head of biceps brachii and the lateral head of triceps brachii. For shoulder external rotation, coracobrachialis and the lateral head of triceps brachii act on the joint function with the lateral head of triceps brachii dominating the joint function both in moment-arm and torque magnitudes. The long head of biceps brachii alone acts on the shoulder-internal-rotation joint function, leaving us to draw no further inferences about this joint function.
Among the many forelimbs muscles that span over the elbow joint, the elbow-flexors, the long head of biceps brachii and extensor carpi radialis longus have the largest moment-arms followed by short head of biceps brachii, brachialis, extensor carpir adialis brevis, and flexor carpir adialis. When we consider the moments, long head of biceps brachii has the largest maximal moment followed by short head of biceps brachii. Other muscles that have a large moment-arm are weak and exert smaller moments. For elbow-extension, all the triceps muscles (triceps brachii lateral, triceps brachii medial, and triceps brachii long head) have relatively strong moment-arms, with the long head of triceps brachii exerting the largest maximal moment. In the current model, we did not include any elbow-pronator muscles. For elbow supination, extensor carpi radialis longus has the largest moment-arm, followed by short head of biceps brachii, extensor carpi radialis brevis, and flexor carpi ulnaris. The short head of biceps brachii becomes the most dominant muscle for elbow supination when we consider moments, followed by extensor carpi radialis longus and flexor carpi ulnaris.
All the wrist flexor muscles have similar maximum moment-arms across the joint function, with extensor carpi radialis brevis, extensor carpi radialis longus, and extensor carpi ulnaris having the most dominant moment-arms. For moments, extensor carpi ulnaris has the strongest moment followed by extensor carpi radialis longus and extensor carpi radialis brevis. Three muscles act on wrist-flexion, of which peroneus longus has the highest moment-arm followed by flexor carpi ulnaris; when we consider moments, however, the dominance is flipped. While flexor carpi radialis, flexor carpi ulnaris, and peroneus longus have the strongest moment-arms for wrist abduction, the strongest moments are exerted by the flexor carpis ulnaris, followed by flexor carpi radialis. Among the three muscles that influence wrist adduction, extensor carpi ulnaris has the strongest moment-arm followed by extensor carpi radialis brevis and extensor carpi radialis longus. Extensor carpi radialis longus has a zero-crossing such that it influences both wrist abduction and adduction. For moments, extensor carpi ulnaris has the highest moment followed by extensor carpi radialis brevis.

C. OPERATIONAL RANGE OF MUSCLE-FIBER LENGTH
The muscle-fiber length (l m ) determines the working region of the muscle in the force-length (FL) curve (Fig. 8B). More specifically, the fiber-length determines if the muscle produces force purely based on muscle contractions or a combination of muscle contraction and passive forces. The range of operation in the FL relationship is determined by the muscle attachment points and parameters. In Fig. 8A we show the working ranges of normalized muscle-fiber lengths (l m ) across all the DoFs each muscle spans. Muscles that span over a single joint tend to have shorter range of operation in the FL relationship. Most muscles in the mouse forelimb and hindlimb operate within the active (l m ≤ 1.0) and passive force (l m > 1.0) regions of the FL curve. Gluteus maximus in the hindlimb is the only muscle that operates completely in the active region of the FL curve. Short head of biceps brachii, brachialis, extensor carpi radialis, and extensor indicis proprius muscles in the forelimb exhibit operation only in the active region. These muscles in the hindlimb and forelimb have no passive/elastic force contributions during movement. In the hindlimb, extensor digitorum longus, extensor hallucis longus, peroneus digiti quarti, and peroneus tertius are the muscles that operate only in the passive region of the FL curve. All the mentioned muscles span over the ankle joint. Pronator quadratus is the only muscle in the forelimb that operates only in the passive region of the FL curve. This is because the optimization of the parameters of forelimb muscles was performed under the constraint that every muscle must have some operation range in the active section of the FL curve. A muscle whose operational range is only in the passive region always produces an elastic force on the joints and the force increases exponentially as the muscle is stretched. Passive forces are useful as they consume no energy to produce movement but if the desired movement is against the passive force, then it requires additional energy to overcome it. Hence, there is an interesting optimum that could be reached to minimize energy consumption. All the other muscles in the hindlimb and forelimb operate across the active and passive regions depending on the pose/configuration of the joints.

D. SENSITIVITY ANALYSIS
Muscle attachments and muscle parameters together determine the overall function and moments produced by the muscle. The high number of muscles and the DoFs each muscle influences makes it a challenging task to perform a comprehensive sensitivity analysis. Here, we use the Sobol method, a variance-based global sensitivity method to perform the analysis (see Section II-D).
It is important to note that in Fig. 9, results of two separate sensitivity analyses are displayed, one for the muscle attachments and one for the muscle parameters. Thus, sensitivity indices for the two cases should be interpreted independently.

1) SENSITIVITY OF JOINT MOMENT-ARMS TO CHANGES IN MUSCLE
ATTACHMENT POINTS-The first observation we can make from the sensitivity analysis to variation in attachment points is that the majority of muscle joint moment-arms are highly sensitive to either P WO (attachment to the parent bone of the joint) or P WI (attachment to the child bone of the joint) and only very few to both. Fig. 10(A-D) shows an example of the variation of moment-arms when the 3D positions of P WO and P WI are varied between ±0.5 mm for gracilis anterior muscle across the 4 DoFs it spans over the hip and the knee joints. For hip flexion-extension and abduction-adduction, moment-arms of gracilis anterior are highly sensitive to P WO . For hip internal and external rotation and knee flexion-extension DoF, the moment-arm of gracilis anterior is more sensitive to changes in P WI than P WO . The same observations are reflected in the comprehensive representation shown in Fig. 9.
Referring to Fig. 4, we can interpret the geometric reason for why a muscle-joint pair is sensitive to either P WO , P WI , or sometimes both. When the attachment points are perturbed to perform the sensitivity analysis, points further away from the joint rotation centers will result in larger changes in muscle length within the range-of-motion of the joint. As we have seen from (2), larger change in muscle-tendon length for the same change in DoF motion results in larger moment-arms. Thus, the proximity of the attachment point to the joint's center-of-rotation will influence its sensitivity. From the analysis (Fig. 10), we can identify those attachment points that are most important and interesting to explore for a given muscle-DoF pair.

2) SENSITIVITY OF JOINT MOMENT TO CHANGES IN MUSCLE
PARAMETERS-In addition to the muscle attachment points, the muscle parameters determine the dynamics of the muscle and the moments it generates on the joints. Estimation of muscle parameters is a challenging process that could potentially lead to modeling errors. This applies to our forelimb model, where the muscle parameters were estimated based on several different data sources (See Section II-C3.b). The sensitivity analysis provides an overview of the parameters that have most influence on the moment generation for each muscle-joint pair. We performed this analysis for four muscle parameters [maximum isometric force F m 0 , optimal fiber length (l m o ), tendon slack length (l t s ), and pennation angle (α o )] within a range of ±10% of their original values. The sensitivities of joint moments to changes in muscle parameters are reported in Fig. 9. We do not show the results for α o , since no joint moment exhibited a significant sensitivity.
From the Hill-type muscle force (see (3) in the Appendix) F m 0 has a linear influence on the muscle force while l m o and l t s have complex, non-linear relationships.
Among the hindlimb muscles that span the hip joint, F m 0 is the parameter influencing the joint moments the most; l m o being the next parameter in a few muscles in this group. No muscle in this group exhibits sensitivity for the choice of l t s parameter. Next, sensitivity of moments caused by muscles that span hip and knee joints are equal due to the changes in

IV. DISCUSSION
We presented a whole-body three dimensional (3D) musculoskeletal model of the mouse with a fully articulated skeletal system actuated by identified musculature for both hindlimbs and forelimbs. Using the model, we performed a systematic and comprehensive analysis of the limb musculature to study their influence on limb joints. We first studied how the muscles influence joint function based on the moment-arm and moments they exert. The analysis gives a comprehensive view to characterize muscle function. Our results reveal that many muscles that span multiple degrees-of-freedom (DoFs) tend to have zero-crossing (i.e., change their function over the DoFs range-of-motion). Examining the muscle-fiber length range showed how the limb muscles distribute their force production in terms of active and passive forces over the joints' complete range-of-motion. We then performed a sensitivity analysis to highlight the crucial parameters in the model and showed how different parameters affect on each muscle-DoF pair in the model. Although the model was based on a number of simplifications and assumptions, it is an important step in the direction of building complex biomechanical, and ultimately neuromechanical models to study motor behaviors and their underlying neuronal control.

A. MUSCLE SYSTEM DEVELOPMENT
Identification and characterization of the muscles operating in complex musculoskeletal systems is a challenging task. It involves a laborious process of extracting individual muscles from the animal, carefully identifying the attachment landmarks such as origin and insertion locations of each muscle and identification of major muscle parameters.
These steps were traditionally performed directly on cadavers [35]; later studies, used microCT scans along with digital segmentation, which lead to more detailed identification about muscle geometry and attachments, especially for deep muscles [49]. In this work, we developed the muscle system of both hindlimb and forelimbs by incorporating data from several studies on the mouse.

1) TRANSFERRING HINDLIMB MUSCLE ATTACHMENT POINTS FROM
CHARLES MODEL-Development of the hindlimb musculature was largely based on the OpenSim single mouse hindlimb model of Charles et al. [31], [49]. Since, this model was only of a single hindlimb, we had to transfer the muscle system to our full mouse model which had a different bone geometry. The first step in this process was to identify the appropriate landmarks of origin and insertion points of each muscle on the new model. To accomplish this task, we setup an automatic process to identify a coordinate transformation between the bones of two models using mesh-registration technique. With this, every attachment point of amuscle defined in a particular coordinate of a bone was transferred to the coordinate frame of the same bone in the current model. The model in Opensim had incorporated muscle wrapping surfaces to better describe the muscle paths. However, in our current framework, muscle paths were approximated as linear polyline paths similar to [18], [26], [48]. To compensate for this approximation, we had to manually introduce additional waypoints to describe the muscle path closer to the original model. With the use of polyline method it was possible to faithfully describe muscle paths. In figure 6, we compared moment-arms of six hindlimb muscles from our current model with the model developed by Charles et al. in OpenSim. Among the six muscles, five of them had used wrapping surfaces in Charles model (except tibialis anterior (TA)). The comparison showed excellent qualitative and quantitative agreement between the moment-arms of the two models, highlighting that the approximation of polyline method captured muscle paths well throughout range-of-motion.  [55]. Thus, the l t s parameter has no direct measurable tendon property of the muscle. This is in accordance with the original formulation of the Hill-type muscle models [46]. Also, Charles et al.

2) SCALING OF HINDLIMB MUSCLE PARAMETERS-After
showed that the measured tendon lengths were either longer or shorter than the estimated l t s values [31]. Because of these observations, we did not impose any constraints to the algorithm to preserve the ratios between l m o and l t s of Charles' model. During future iterations of the model improvement, l t s parameter could be estimated from animal experiments [22].

3) MODELING THE MOUSE FORELIMB IS MORE CHALLENGING THAN THE
HINDLIMB-Developing the forelimb muscle system was more challenging due to the lack of available biomechanical studies of the mouse or even rat forelimbs. To the best of our knowledge, Mathewson et al. [54] was the only published work that measured some of the muscle properties necessary to model Hill-type muscles. But, since the goal of their work was not building a simulation model of the mouse forelimb, information about muscle attachments were not reported.
Delaurier et al. [51] developed a 3D model of mouse embryonic forelimb using Optical Projection Tomography and digital segmentation. We transferred the attachment points for each muscle using their 3D atlas to our model. The choice of muscle was based on muscle data reported in [54]. This limited the forelimb model to mostly distal muscles. Proximal muscles around the shoulder that had origins from the spine were omitted from the model.
Mathewson et al. [54] reported l m o and α 0 for the forelimb muscles. Unlike in the hindlimb case, we could not employ the algorithm to scale length related parameters from Modenese et al. [53] because of a lack of information about muscle lengths at different model poses.
The l m o parameter was thus scaled based on the ratio of average muscle-tendon length for while α 0 was used as reported. The final missing parameter l t s was estimated based on the extended algorithm by Manal and Buchanan [55]. (Refer to Appendix C for details on the changes incorporated to the algorithm.) Determining l t s is one of the biggest bottlenecks in Hill-type muscle parameterization as there is no method to experimentally to estimate it [31]. The closest experimental method to estimate l t s is described by Cox et al. [22]. However, performing such experiments are extremely difficult in the mouse because of their relatively small size. The modifications we proposed to the numerical method by Manal and Buchanan [55] should improve the reliability of these estimates for the mouse as well as for other musculoskeletal model of animals.
To the best of our knowledge there is no published work that characterizes either the muscle properties or muscle attachments for forelimb muscles attaching to the scapula or the spine; neither for mice nor rats. Hence, these muscles were considered beyond the scope of this work.

B. MUSCLE MOMENT-ARMS AND MOMENTS
The model includes 59 distinct forelimb and hindlimb muscles. It is a challenging task to provide a comprehensive analysis of a complex model such as this. Conventionally, the practice is to report and describe the relationship between moment-arm and joint movement for each muscle individually. This means that it is necessary to assume the function of a muscle a priori. Instead of making this assumption, we reported a global view of the possible roles of a muscle based on moment-arm and moment (Fig. 7). With this representation one can quickly identify the function of a muscle and observe the contribution of different muscles to a particular degree of freedom.
Previous musculoskeletal modeling studies have reported the behavior of zero-crossing of muscle moment arms in several animals such as cat [61], mouse [31], rat [35] and ostrich [23]. Young et al. [61] speculated that muscles with a zero-crossing moment-arms could intrinsically stabilize the joints around which they change sign without any need for extra neural commands. For the mouse hindlimb model [31] several muscles were reported to have a zero-crossing. But their analysis was limited to the assumption of functional roles assigned to each muscle a priori. Here, we identified more muscles that have zero-crossings. For example, previously only pectineus muscle was reported to have a zero-crossing for hip flexion-extension. From Fig. 11A,B we can observe that in addition to pectineus we have adductor brevis, gemellus, gluteus maximus (ventral), obturator externus and internus and quadratus femoris muscles have zero-crossing. Similarly for the different joints in forelimb and hindlimb, we can observe many muscles exhibiting the zero-crossing behavior.
A shortcoming of the representation we proposed is that it does not show joints' angle information: it is not possible to see at what joint angle the zero-crossing occurs. This limitation is only in the visual representation in this article, and detailed moment-arms of all the possible muscle and joint combinations are available. We encourage the readers to use the data and plotting tools in the code-repository (see section V) to extract the detailed plots.
Along with moment-arms, we also reported the muscle moments across each DoF (Fig. 7). By presenting both moment-arms and moments together, we can quickly observe how the role and importance of a muscle in actuating a particular DoF can change. Muscles with large maximum isometric force naturally become the dominating muscle for the DoF. While the moment-arm is only dependent on the muscle attachments and joint position, moments also depend on the dynamics of the muscle-most importantly muscle activation. From a neural control point of view this is very interesting: the nervous system can operate with a single strong muscle and/or utilize several different weaker muscles together to produce the same movement. Co-activation of multiple synergistic muscles with respect to a specific function could also result in a stabilization of the joint in the other DoFs.
The current muscle moment-arm and moment results were computed while the joint of interest was rotated within the range of motion and while keeping all other joints in their default position. This has a strong influence on the results. It limits the scope of the analysis to a particular pose of the model. This is often the case in biomechanical model analysis as it becomes very complex to interpret and represent the relationship of different joints with muscle moment and moment-arm. Young and colleagues [26], [61] studied the coupled effect of a bi-articular muscle on joints. Also, we see from Fig. 7 that many muscles operate on more than two degrees of freedom, especially muscles that span joints with multiple degrees of freedom (e.g., hip, ankle, shoulder, or wrist). With the condensed representation ( Fig. 7), it is possible to generate plots to study the relationships at various model poses.

Muscle-fiber length (l m ) is a state variable in describing the muscle dynamics (see Appendix
A). It depends on the muscle-tendon length (l mt ), the muscle contraction dynamics, and the muscle activation. l m is thus a very interesting variable to study. In this paper, we reported the operational range of normalized muscle-fiber length (l m ) for each muscle over all the joints it spans in Fig. 8. (Unlike the moment-arm and moment representation in Fig. 7, Fig.  8 incorporates all the possible joint poses for a single muscle.) The range of l m reflects the choice of the tendon slack length l t s and optimal fiber length (l m o ) parameters. In future, if we can obtain experimental data of fiber-lengths at the different limb postures. This data can be used to validate the model by making sure that the experimental measurements fall within the predicted ranges of l m .
Analyzing the l m is also useful in understanding the behaviors by being able to characterize if a muscle is operating in the ascending, plateau or the descending region in the force-length curve (Fig. 7B). In smaller animals, such as mice, due to their low weight and small inertias, the effect of gravity on the body is negligible [62], [63]. With lower influence of gravity, the forces produced by the passive elements become more significant. This has two important consequences on the neural control circuits. One, it is less essential (compared to larger vertebrates) for the neural system to monitor the direction of gravity on the limbs during movements. Second, the momentum of the limbs is less useful during the swing phase during locomotion. For further discussion on the influence of body size on neural control refer to Hooper [62]. Thus, passive forces in mice play an important role in movement generation and it would be informative to use the muscle-fiber length range plots to explore it.
In [64], the authors reported the l m of human leg muscles during different speeds of walking and running. They observed that the l m has a wide range of operation in the force-length curve for different speeds of walking and running. Fig. 8 for the mouse model also highlights that most muscles with the current muscle parameters operate both in the ascending and descending region of the force-length allowing for the kind of variability observed in humans.

D. SENSITIVITY ANALYSIS
While developing a complex model like the one described here, it is important to identify the critical parameters that influence the overall performance of the model. To this aim, we performed a variance based global sensitivity analysis using the Sobol method to systematically study the influences of muscle attachment points and muscle parameters on the moment-arms and moments.
In our analysis of sensitivity of moment-arms to muscle attachment points, we considered only those attachment points (P WO and P WI ) that influence the muscle-length effect on the joint of interest. In Charles et al. [31], the analysis was performed directly at the origin and insertion points, and the perturbations applied to the attachments were along a particular direction. In contrast, in our work, we applied perturbations to the attachments in all three directions. Similar to the observations made by Charles et al. in [31], we unsurprisingly observed that the attachment points located farther from the joint rotation centers had the greatest effect on the joint moment-arms. With the polyline approximation it therefore becomes very important to identify the intermediate points as accurately as possible to describe the muscle paths.
We further analyzed the sensitivity of joint moments to changes in Unsurprisingly, moments across many muscle-DoF pairs were most sensitive to F m 0 with l t s being the next most significant parameters. Our observations presented in Section III-D2 agree with the sensitivity analysis by Charles et al. [31]. Charles et al. observed that muscles with larger l t s compared to their l m o had their moments more s sensitive to l t s than F m 0 . In Table   2 we reported the ratios of . This in accordance with the previous studies not only in mouse [31] but also for chimpanzees [65] and humans [48], [66].
The same applied to the forelimb muscles as well. Although, in the forelimb we have very few muscles whose ratios of Overall, with our comprehensive analysis we have highlighted the most important parameters in the model. This allows future researchers to identify and work with these important parameters more critically. The identification of the critical parameters also allows for adapting and fine tuning the model in case of numerical optimization for specific behaviors.

E. MODEL LIMITATIONS
In this work, we present the most complete musculoskeletal model of the mouse. Yet, as with any model, the model construction was possible only because of some simplifications and assumptions. The skeletal model of the mice was developed from anatomical references rather than from an actual microCT scan. This allowed us to generate a more generic representation of the mouse skeleton but meant that not all anatomical features on the bones were captured. In future iterations, a skeletal model based on CT scans will further allow for better model validation. Here, we calculated the inertias of the bones using a bounding box method and assumed a uniform density of water along the body. We did not consider air cavities with different density in the lungs and head in our model. This introduces variation of inertial properties of certain bones and the overall center-of-mass of the mouse musculoskeletal model. Joint rotations were identified manually in this model and use joints such as revolute or spherical joints to represent the different DoFs in the model. However, in animals, joints are more complex and often difficult to model. Any studies that require validated joint movements should try to extend the current version of the model by including the necessary complexities in joint modeling for the particular study.
For the hindlimb muscles, we were able to build on, compare and validate our model with the previous mouse [31] and rat [35] hindlimb models. However, similar validation was not possible for the forelimb. Because of the lack of previous forelimb studies, we had to estimate both the attachment points and the muscle parameters (F m 0 , l m o , l t s and α o ). The comprehensive sensitivity analysis presented in this work provides information about the important parameters to be critically explored in the model. Further experimental measurements from the mouse forelimb are necessary to improve the forelimb model.

F. MODEL USE AND FUTURE WORK
A musculoskeletal model complements animal experiments by providing the data (EMG, afferent firings, interactions forces between the body and the environment) that is challenging to collect, and they are essential to setup and study closed-loop neuromechanical simulations (Fig. 1). The whole-body model of the mouse presented in this work has the necessary components to contribute to both.
Locomotion is a result of whole-body movement with the neural circuits integrating feedforward and sensory-driven strategies to generate the necessary muscle activation signals. The limb musculature modeled and analyzed in this work presents an excellent platform to use the model to setup predictive simulations to connect computational models of neural circuits to drive the musculoskeletal model. Alternatively, inverse kinematics and inverse dynamics approaches can be used to estimate kinematics (e.g., joint angles), kinetics (e.g., joint moments) and proprioceptive sensory information (e.g., activities of muscle spindles and Golgi tendon organs) from experimental whole-body trajectories made possible now by markerless pose estimation methods such as DeepLabCut [67].
The majority of the previous locomotion studies in mouse have been limited to straight forward locomotion. Exploring other locomotor regimes such as turning is now experimentally possible [68] and the current 3D model has the necessary DoFs to replicate similar behaviors in simulation. The current model also allows to study motor behaviors that does not involve the whole-body movements like reaching and grasping.
While the possible use cases of our model are plenty, it still is only a preliminary step towards a more robust computational model. As mice are one of the most significant experimental animal models to study normal and pathological motor behaviors, it is extremely important to develop computational models that can complement these studies.
The open-source and modular musculoskeletal model presented here offers an opportunity for a community driven approach that can collectively improve and rigorously validate the different components of the model with experimental data. In future iterations of the model, a systematic identification the full set of forelimb muscles along with the spinal muscles will increase the usability of the model in even more complex scenarios and grow towards a more complete model like the ones available for humans.

V. CODE AVAILABILITY
All the resources and code for the development and analysis of the mouse model can be found at https://gitlab.com/paper_submissions/mouse_biomechanics_paper JONATHAN ARREGUIT received the B.Sc. and M.Sc. degrees in microengineering (orientation in robotics and autonomous systems) and a minor in neuroprosthetics from the Swiss Institute of Technology of Lausanne (EPFL). He is currently pursuing the Ph.D. degree with the Biorobotic Laboratory (BioRob), and is involved in a project funded by the Human Frontier Science Program (HFSP) where he mainly studies the combination of central pattern generators with sensory feedback in neuro-mechanical models for amphibian locomotion using physics simulations with complex environments.
During his master's program, he completed a semester project and a master project at BioRob, related to state estimation and motion planning for humanoid robots. His research interests include bio-inspiration, locomotion, connectionism, and simulation of neuro-mechanical models. DIMITRI RODARIE received the degree in computer sciences engineering from Polytech Lyon, France, and the master's degree in artificial intelligence from Universite Claude Bernard, Lyon, France. He is currently pursuing the Ph.D. degree with the Connectomics Team, Simulation Neuroscience Division, Blue Brain Project, EPFL, Geneva. His research interests include data-driven generation, simulation, and analysis of whole mouse brain models. He is also working on embedding detailed microcircuit models in to simplified neuron networks and combining brain models with virtual musculoskeletal bodies.
interests include computational modeling of neurons, neural circuits and neural systems at different temporal and spatial scales, central pattern generators, neural oscillations, neural control, visual perception, and image recognition. The long-term goal of his studies is to understand the key issue of neural control of movement: how different cellular, network, and systems' neural mechanisms are integrated across multiple levels of organization to produce motor behavior and to adapt this behavior to various external and internal conditions. Together with his colleagues, he has developed a series of most comprehensive and realistic models of neural circuits in the brainstem and spinal cord responsible for control of respiration and locomotion. He has been a Professor with EPFL, since 2002, where he is currently the Head of the Biorobotics Laboratory. His research interests include robotics, computational neuroscience, nonlinear dynamical systems, and applied machine learning. He is interested in using numerical simulations and robots to gain a better understanding of animal locomotion, and in using inspiration from biology to design novel types of robots and controllers.Heisalsoinvestigatinghowtoassistpeoplewithlimitedmobility using exoskeletons and assistive furniture.

APPENDIX A: HILL-TYPE MUSCLE MODEL
The force produced by the contractile and parallel elements was modeled as, F m 0 is the maximum isometric force the muscle is capable of producing. a(t) is the muscle activation. f l l m defines the force-length relationship curve during isometric contractions with l m being the normalized muscle-fiber length. f v v m defines the force-velocity relationship curve during isotonic contractions with v m being the normalized muscle-fiber length. Passive forces are generated in the muscle when stretched beyond a threshold length, this is represented by the passive-force-length given by f P E l m . βv m represents the additional damping in the muscle with β being the damping coefficient, set to 0.1 by default.
The activation dynamics (a(t)) and passive-force-length (f PE ) are described based on the implementation from Opensim [69].
where u(t) is the neural-excitation input signal to the muscles bounded within a range of 0 and 1, and τ act is the timeconstant (Note that we used a constant value for τ act instead of varying it as a function of u(t) and a(t) as described by Thelen [69].) f P E l m = exp k P E l m − 1 /ϵ 0 m − 1 exp k P E − 1 (5) The active-force-length (f l ) and force-velocity (f v ) relationships described below are based on the description by De Groote et al. in [70] f l l m = ∑ i = 1 Since, the tendon is assumed to be rigid and inextensible, the force along the tendon (F t ) equates to the force produced by the contractile element described in (3) as,

APPENDIX B: HINDLIMB MUSCLE PARAMETER SCALING
Optimal muscle-fiber length (l m o )and muscle-tendon slack length (l t s ) were scaled from Charles model [31] (referred to as reference model) to our model (referred to as target model). To do the scaling, we used the algorithm proposed by Modenese et al. in [53], steps for which are described below, We solved (9) using the least squares method implemented in Scipy v1.70.

APPENDIX C: MUSCLE-TENDON SLACK LENGTH ESTIMATION
Muscle-tendon slack length for the forelimb muscles were estimated based on the modified algorithm proposed by Manal and Buchanan in [55].
For the purpose of estimating the muscle-tendon slack length, we consider the tendons to be elastic and non-rigid whose force-tendon length f t l t (with l t = l t l t s being the normalized muscle-tendon length) characteristics are defined by the formulation by De Groote et al. in [70] as, f t l t = c 1 exp k T l t − c 2 − c 3 Hill-type muscle model parameters. Substituting (10) in (8), The total length of the muscle-tendon unit is, l mt = l m cos(α) + l t Considering a maximal activation (a(t) = 1.0) under isometric conditions (v m = 0.0), from (10), (11) and (14), l t s can be written as, c 1 exp k T l t − c 2 − c 3 = f l l m + f P E l m cos(α) (15) l t = c 2 + log f l l m + f P E l m cos(α) + c 3 /c 1 k T (16) l t s = l mt − l m l m o cos(α) c 2 + log f l l m + f P E l m cos(α) + c 3 /c 1 /k T Here we uses the tendon force-length relationship formulated by De Groote et al. in [70] instead of the one described Manal and Buchanan in [55].
Given l m at l mt and knowing l m o , we can compute tendon slack length l t s from (17). But, we do not know l m in practice. Manal and Buchanan proposed to use numerical methods to estimate l m at short, long and mid range of l mt such that l t s is to be a constant, which was formulated as a minimization problem.
minimize l t s1 − l t s2 2 + l t s1 − l t s3 2 + l t s2 − l t s3 2 subject to, l m ≤ 1.0 We extended this formulation by Manal et al. for l t s in [55] by including the passive forces (f P E l m ) in (17) so that we include the complete force-length curve for estimation. The formalization described above is limited to the muscles' relationship with only a single degree-of-freedom. We extend this to all the degrees-of-freedom (ndof) where the muscle has a significant moment-arm. This is done by estimating l m at short, long and mid range of l mt for each degree-of-freedom the muscle influences. The new minimization problem is defined as, minimize ∑ n = 1 ndof l t s1 − l t s2 2 + l t s1 − l t s3 n2 + l t s2 − l t s3 2 n (18) subject to, 0.6 ≤ l m ≤ 1.4 (19) We used differential evolution, a stochastic population based optimization technique implemented in Scipy v1.70 to minimize the objective in (18) for each muscle in the forelimb. A population size of 20, maximum iterations of 1000 and relative and absolute tolerances were set to 1e −10 . Components of a closed-loop neuromechanical simulation. Movements in animals arise due to complex interactions between the nervous system, the musculoskeletal system and the environment. A neuromechanical model includes the neural and biomechanical components along with their interactions. Neural models (spinal and brain circuits) produce the necessary instruction signals (motoneuron activations) for a specific movement. Muscle models respond to the neural signals by producing forces that act on the skeletal model and cause movements. The skeletal model and the body interact with the environment to produce reaction forces. The sensory organ models encode the state of the movement (somatosensory afferent feedback signals) and transmit this information back to the nervous system which then adapts the instruction signals to sensed perturbations or external forces.  Hill-type muscle model describing the force generation by the muscles. The contractile element (CE) or the active element produces active contraction forces. The parallel element (PE) prevents the muscle from over stretching the muscle-tendon unit during normal operation. The series element (SE) represents the series elasticity of the muscle, including the muscle-tendon. Contractile element length or fiber-length (l m ) is the length of muscles fibers. Pennate muscles are defined by the pennation angle α o . Series tendon-length (l t ) is the length of series element. The total length of the muscle (l mt ) is defined as l mt = l m cos(α) + l t .   Maximum moment arms (A, C) and moments (B, D) for each muscle and joint function. Grey boxes indicate joint functions for joints the muscle spans but has zero influence over. For example, a pure hip flexor muscle is considered to be spanning over both hip flexion and extension joints, but the corresponding hip extension will be in grey to indicate that the muscle has no influence on hip extension. The moment-arms and moments are normalized by the muscle which has the maximum influence in the group. ϵ indicates a very low non-zero positive value.   Mouse skeletal segments grouping, joint names and joint ranges.  Hindlimb and forelimb muscle parameters.  Maximum isometric force of the muscle. Hindlimb muscle forces are the same as reported in [31]. Forelimb muscle forces are computed by multiplying the physiological cross-sectional area (PCSA)

Muscle-tendon
reported in [54] by isometric stress (σ = 0.3 N/mm 2 ) [31], [60]. (b) Optimal fiber length for hindlimb scaled to the current model from [31] using the scaling based on [53] and forelimb scaled to the current model [54]. (c) Tendon slack length is obtained in the same way as described for the optimal fiber length. (d) Pennation angle for hindlimb is obtained from [31] and for forelimb from [54]. (e) Ratio of optimal fiber length and tendon slack length.
(f) Sorting of muscles into groups based on the joints they span. For example, Biceps femoris posterior (BFP) belongs to Hip-Knee group as it spans over one or more hip and knee joints. The same order is used in the rest of this work.