Gait Generation and Its Energy Efficiency Based on Rat Neuromusculoskeletal Model

Changing gait is crucial for adaptive and smooth animal locomotion. Although it remains unclear what makes animals decide on a specific gait, energy efficiency is an important factor. It has been reported that the relationship of oxygen consumption with speed is U-shaped for each horse gait and that different gaits have different speeds at which oxygen consumption is minimized. This allows the horse to produce energy-efficient locomotion in a wide speed range by changing gait. However, the underlying mechanisms causing oxygen consumption to be U-shaped and the speeds for the minimum consumption to be different between different gaits are unclear. In the present study, we used a neuromusculoskeletal model of the rat to examine the mechanism from a dynamic viewpoint. Specifically, we constructed the musculoskeletal part of the model based on empirical anatomical data on rats and the motor control model based on the physiological concepts of the spinal central pattern generator and muscle synergy. We also incorporated the posture and speed regulation models at the levels of the brainstem and cerebellum. Our model achieved walking through forward dynamic simulation, and the simulated joint kinematics and muscle activities were compared with animal data. Our model also achieved trotting by changing only the phase difference of the muscle-synergy-based motor commands between the forelimb and hindlimb. Furthermore, the speed of each gait varied by changing only the extension phase duration and amplitude of the muscle synergy-based motor commands and the reference values for the regulation models. The relationship between cost of transport (CoT) and speed was U-shaped for both the generated walking and trotting, and the speeds for the minimum CoT were different for the two gaits, as observed in the oxygen consumption of horses. We found that the resonance property and the posture and speed regulations contributed to the CoT shape and difference in speeds for the minimum CoT. We further discussed the energy efficiency of gait based on the simulation results.

Changing gait is crucial for adaptive and smooth animal locomotion. Although it remains unclear what makes animals decide on a specific gait, energy efficiency is an important factor. It has been reported that the relationship of oxygen consumption with speed is U-shaped for each horse gait and that different gaits have different speeds at which oxygen consumption is minimized. This allows the horse to produce energy-efficient locomotion in a wide speed range by changing gait. However, the underlying mechanisms causing oxygen consumption to be U-shaped and the speeds for the minimum consumption to be different between different gaits are unclear. In the present study, we used a neuromusculoskeletal model of the rat to examine the mechanism from a dynamic viewpoint. Specifically, we constructed the musculoskeletal part of the model based on empirical anatomical data on rats and the motor control model based on the physiological concepts of the spinal central pattern generator and muscle synergy. We also incorporated the posture and speed regulation models at the levels of the brainstem and cerebellum. Our model achieved walking through forward dynamic simulation, and the simulated joint kinematics and muscle activities were compared with animal data. Our model also achieved trotting by changing only the phase difference of the muscle-synergy-based motor commands between the forelimb and hindlimb. Furthermore, the speed of each gait varied by changing only the extension phase duration and amplitude of the muscle synergy-based motor commands and the reference values for the regulation models. The relationship between cost of transport (CoT) and speed was U-shaped for both the generated walking and trotting, and the speeds for the minimum CoT were different for the two gaits, as observed in the oxygen consumption of horses. We found that the resonance property and the posture and speed regulations contributed to the CoT shape and difference in speeds for the minimum CoT. We further discussed the energy efficiency of gait based on the simulation results. Keywords: rat, walk, trot, energy efficiency, central pattern generator, muscle synergy, neuromusculoskeletal model

INTRODUCTION
Animals can generate adaptive and smooth locomotion in various conditions. One important strategy for such locomotion is the use of different gaits. For example, quadruped animals walk, amble, trot, pace, canter, and gallop. Although gait is the motor outcome of a complicated and redundant musculoskeletal system controlled by the central nervous system, it is largely unclear what makes animals decide on a gait. One important factor for deciding gait is the energy efficiency of locomotion; that is, animals want to minimize the cost of transport (CoT). In particular, it has been reported that the relation between oxygen consumption and speed is U-shaped for each horse gait and that different gaits have different speeds at which oxygen consumption is minimized (Figure 1) (Hoyt and Taylor, 1981). Walking, trotting, and galloping are energy-efficient at low, middle, and high speeds, respectively. Walking and trotting share a common speed range, as do trotting and galloping. Therefore, horses can produce energy-efficient locomotion over a wide speed range by changing their gait. However, the underlying mechanisms making the oxygen consumption U-shaped and the speeds for minimum consumption different between gaits remain unclear.
Locomotion is generated through interactions between the central nervous system, musculoskeletal system, and environment. It is difficult to fully analyze the locomotor mechanism with animal data alone. Recently, modeling studies have attracted attention because physiological findings and hypotheses can be used to develop reasonably realistic motor control models, and biomechanical and anatomical findings can be used to construct detailed musculoskeletal models (Ivashko et al., 2003;Yakovenko et al., 2004;Ekeberg and Pearson, 2005;Nishii, 2006;Aoi et al., 2013a;Fukuoka et al., 2015;Hunt et al., 2015;Aoi and Funato, 2016;Markin et al., 2016;Fujiki et al., 2018). Motor control and musculoskeletal models are integrated to produce locomotion through forward dynamics simulation. This allows the locomotor mechanism to be examined from a dynamic viewpoint.
In this study, we investigated the energy efficiency of gait using a rat neuromusculoskeletal model. Specifically, we constructed a musculoskeletal model composed of the trunk, forelimbs, and hindlimbs based on anatomical data. This model is an improvement on our previous rat hindlimb model (Aoi et al., FIGURE 1 | Oxygen consumption of horses in walking, trotting, and galloping. 2013a). We also improved our previous motor control model to control the rat four-limb model. The motor control model was developed based on the hypothetical two-layer central pattern generator (CPG) model at the spinal cord level (Burke et al., 2001;Rybak et al., 2006) and the muscle synergy hypothesis (Tresch et al., 1999;Todorov and Jordan, 2002;d'Avella et al., 2003;Ting and Macpherson, 2005;Ivanenko et al., 2006;Drew et al., 2008;Takei et al., 2017), which describes a simple control strategy for redundant motor systems. Furthermore, we incorporated movement regulation models at the levels of the brainstem and cerebellum through brainstem descending pathways. We simulated the walking of our model and compared the simulation results with animal data. In addition, we simulated trotting and changed the speed of each gait using simple motor control strategies. We calculated the CoT of walking and trotting for the generated speeds and, in this paper, we discuss the energy efficiency of gait based on the simulation results.

Musculoskeletal Model
We developed a rat musculoskeletal model based on our previous model, which focused on the hindlimbs without incorporating the forelimbs (Aoi et al., 2013a). The skeletal part of the model consists of eleven rigid links representing the trunk (including the head), forelimbs (two links), and hindlimbs (three links), as shown in Figure 2. This model is two-dimensional, and the walking behavior is constrained to the sagittal plane. When the brachium and antebrachium are in a straight line and perpendicular to the trunk, the shoulder angle is 120 • and the elbow angle is 180 • . When the thigh, shank, and foot are in a straight line and perpendicular to the trunk, the hip angle is 120 • and the knee and ankle angles are both 180 • . The joint angles increase as the joints extend. We modeled the contact between the limb tips and the ground using viscoelastic elements. We derived the equations of motion using Lagrangian equations and solved them using the fourth-order Runge-Kutta method with a time step of 0.02 ms.
For the muscle part of the model, we used six principal muscles for each forelimb: four uniarticular, namely shoulder extension (supraspinatus, SSP), shoulder flexion (spinoltoideus, SPD), elbow flexion (brachioradialis, BR), and elbow extension (triceps lateral head, TRIL), and two biarticular, namely shoulder extension and elbow flexion (biceps, BIC) and shoulder flexion and elbow extension (triceps, TRI), as shown in Figure 2. We used seven principal muscles for each hindlimb: five uniarticular, namely hip flexion (iliopsoas, IP), hip extension (gluteus maximus, GM), knee extension (vastus lateralis, VL), ankle flexion (tibialis anterior, TA), and ankle extension (soleus, SO), and two biarticular, namely hip extension and knee flexion (biceps femoris, BF) and knee flexion and ankle extension (gastrocnemius, GA). The moment arms of the muscles around the joints are constant, regardless of joint angles. Each muscle generates muscle tension F m (m = SSP, SPD, BR, TRIL, BIC, TRI, IP, GM, VL, TA, SO, BF, and GA) through contractile and passive where F max m is the maximum muscle tension, a m is the muscle activation (0 ≤ a m ≤ 1), F l m is the force-length relationship, F v m is the force-velocity relationship, and F p m is the passive component. The muscle lengths were normalized by l max m , which was set so that all uniarticular muscles had a length of 85% of l max m and all biarticular muscles had a length of 75% of l max m at a neutral posture with the shoulder at 60 • , the elbow at 85 • , the hip at 70 • , the knee at 90 • , and the ankle at 100 • . In addition, 2 • of joint motion corresponded to 1% of muscle length change, except for BIC and GA (4.5 • at the shoulder for BIC, 1.5 • at the ankle or 4.5 • at the knee for GA). The muscle contractile velocities were normalized by 1.8l max m . The muscle activation a m is determined through where τ act and τ deact are respectively, activation and deactivation time constants (11 and 18 ms, respectively) and u m is the motor command determined in the motor control model.

Motor Control Model
We developed a motor control model based on our previous work (Aoi et al., 2013a). It consists of the following two components: 1. a movement generator, which produces motor commands in a feedforward fashion at the spinal cord level to create periodic limb movements based on the muscle synergy hypothesis and 2. a movement regulator, which creates motor commands to regulate locomotor behavior in a feedback fashion at the brainstem and cerebellar levels based on proprioceptive and somatosensory information. The motor command u m is the summation of the two components from the movement generator and the movement regulator, namely u

Movement Generator
The movement generator is based on the hypothetical twolayer CPG model composed of a rhythm generator (RG) network, which produces rhythm and phase information for motor commands, and a pattern formation (PF) network, which produces spatiotemporal patterns of motor commands (Burke et al., 2001;Rybak et al., 2006). For the RG model, we used four simple phase oscillators, each of which produces a basic rhythm and phase information for the corresponding limb. We used φ j i (i = left, right, j = fore, hind) for the oscillator phase (0 ≤ φ j i < 2π), which follows the dynamics given bẏ where T is the gait cycle duration and K 1 and K 2 are gain parameters. The second term on the right-hand side ensures that the left and right limbs move in antiphase to maintain interlimb coordination. The third term on the right-hand side ensures that the ipsilateral limbs move in relative phase of to maintain interlimb coordination.
For the PF model, we determined the motor commands necessary to produce periodic limb movements in accordance with the corresponding oscillator phase based on the muscle synergy hypothesis, which suggests that the linear combination of only a small number of basic signals produces a large portion of motor commands in animal locomotion Dominici et al., 2011;Markin et al., 2012;Catavitello et al., 2015;Rigosa et al., 2015). Specifically, we used four rectangular pulses p i (i = 1, . . . , 4) for each limb, which are given by where i and i (i = 1, . . . , 4) are the onset and end phases of the pulse, respectively, and we omitted the suffix of φ. p 1 , p 2 , p 3 , and p 4 contribute to early extension, late extension, early flexion, and late flexion, respectively, as shown in Figure 3 [extension and flexion phases start at 1 (= 0 rad) and 3 , respectively]. We used the same values of i and i for the four limbs irrespective of whether they were forelimb or hindlimb. The motor command u Syn m of the movement generator is given by where w m,i (i = 1, . . . , 4) is the weighting coefficient.

Movement Regulator
At the levels of the brainstem and cerebellum, locomotor behavior is regulated based on proprioceptive and somatosensory information (Takakusaki, 2017). For the rat, it is crucial to maintain body height and forward speed during locomotion (Figure 4). For simplicity, we focused on these two factors.
For the body height, we used simple feedback control for the standing limb. For the forelimbs, we used the BR and TRIL muscles to maintain the shoulder height. The motor command p height m (m = BR and TRIL) is given by where h Shoulder andḣ Shoulder are the shoulder height and its rate, respectively, h Shoulder 0 is the reference height, and K  (8) where h Hip andḣ Hip are the hip height and its rate, respectively, and h Hip 0 is the reference height. For the forward speed, we used simple feedback control for the standing limb. We used the SSP, SPD, IP, GM, TA, and SO muscles to maintain speed. The motor command p speed m (m = SSP, SPD, IP, GM, TA, and SO) is given by where v is the forward speed, v 0 is its desired value, and K speed m is a gain parameter.
The summation of these elements produces the motor command of the movement regulator. Because regulation is managed at the brainstem and cerebellar levels, the command signals are delayed and the motor command u Reg m of the movement regulator is given by and τ Delay (= 15 ms) is the delay between receiving the transmission of proprioceptive and somatosensory information at the brainstem and cerebellar levels and sending the motor command to the spinal cord level.

Changing Gait and Speed
In this study, we focused on two gaits, namely walking and trotting. They are mainly classified by the footfall sequence. Specifically, four limbs move out of phase in walking, and diagonal limbs are paired in trotting (Figure 5). Right and left limbs move in antiphase in both walking and trotting. The major difference between the gaits is the relative phase between the ipsilateral limbs. To change the relative phase of the limb movements, we changed the relative phase of the musclesynergy-based motor command u Syn m between the ipsilateral limbs by changing in (4). In particular, we used = π/2 for walking and = π for trotting. Animals change the gait cycle duration to vary speed, where the duration of the flexion phase for swinging the limb FIGURE 5 | Schematic diagram of footfall sequence for walking and trotting. Color bars indicate stance phase. Right and left limbs move in antiphase. Relative phase between forelimb and hindlimb is π/2 for walking and π for trotting (LF, left fore; LH, left hind; RF, right fore; RH, right hind). remains almost unchanged and the duration of the extension phase for supporting the body and producing the propulsive forces is changed substantially (Goslow et al., 1973;Heglund and Taylor, 1998;Clarke and Still, 1999;Górska et al., 1999;Yakovenko et al., 2005). In this study, we changed the speed by changing the duration of the extension phase T ex (= βT) using β while keeping the duration of the flexion phase T fl unchanged (T = T ex + T fl = T fl /(1 − β)), as shown in Figure 6A. For the nominal speed, which we determined from animal data as explained below, we usedβ,T ex ,ˆ i ,ˆ i ,ŵ m,i (i = 1, . . . We decreased (increased) the extension phase duration to increase (decrease) the speed, which decreased (increased) the duration of pulses of the extension phase. To prevent the model from decreasing (increasing) the speed, we increased (decreased) the weighting coefficients w m,i (i = 1, 2) of the muscle-synergy-based rectangular pulses for the extension phase ( Figure 6B) as As we changed the locomotion speed, we also changed the reference height (shoulder, hip) and speed for the movement regulator (Figures 6C-E) as where α Shoulder , α Hip , and α Speed are coefficients.

Parameters for the Musculoskeletal Model
To determine the physical parameters of the musculoskeletal model, we used seven adult male Wistar rats (body weight: 125 ± 10 g). The rats were deeply anesthetized, and their musculoskeletal features were measured. The experiments were approved by the Ethical Committee for Animal Experiments at the University of Tokyo and carried out in accordance with the Guidelines for Research with Experimental Animals of the University of Tokyo and the Guide for the Care and Use of Laboratory Animals (NIH Guide). For the skeletal model, we measured several physical parameters of the rats, such as masses, joint positions, and distances between joints, and determined the model parameters from these measurements, as shown in Table 1. For the muscle model, we first electrically stimulated individual muscles and determined which joint movements were needed to verify our musculoskeletal model. We measured the attachment, direction, and physiological cross-sectional area (PCSA) for each muscle and determined the model parameters from these measurements, as shown in Table 2, where the maximum muscle tension F max m was determined based on the measured PCSA and the moment arms were determined from the center of the range of joint movement during locomotion.

Comparison With Animal Data
To evaluate our neuromusculoskeletal model, we compared the simulation results for walking with animal data. We used the joint angles of the hindlimbs measured in Aoi et al. (2013a), where rats walked on a treadmill at a speed of 0.4 m/s, and the joint angles of the forelimbs measured in Aoki et al. (2013), where intact rats walked at the average speed of 0.36 m/s in a custom-made runway box (length: 140 cm; width: 14 cm).
We used the electromyographic (EMG) data measured from the muscles of the hindlimbs in Aoi et al. (2013a) and the EMG data measured from two muscles (BIC and TRI) of the forelimbs in Aoki et al. (2013). Because we could not find EMG data for four muscles (SSP, SPD, BR, and TRIL) of the forelimbs of rats, we used EMG data for these muscles in cats, whose gait and joint movements are similar to those of rats, given in Drew et al. (2008), where cats walked on a treadmill at a speed of 0.35-0.45 m/s. In the comparison with the simulation results, we showed the EMG data so that their magnitudes are similar to those of simulated muscle activities.

Evaluation of Cost of Transport
The energetic cost of locomotion for our simulation results for walking and trotting was estimated based on the mechanical energy exerted by muscles. Based on previous work (Ogihara et al., 2011), we calculated the CoT ε as where v m is the contracting velocity of the muscle (positive for contraction), and [x] + is x if x ≥ 0 and 0 if x < 0. η + and η − are the positive and negative mechanical work done by muscles, respectively, for one gait cycle duration. The negative mechanical work was divided by four based on Margaria et al. (1963), Elmer and LaStayo (2014). D is the moving distance of the model for one gait cycle duration, which corresponds to the stride length. In this study, the motor command u m is generated by three elements: rectangular pulses u Using these values, we calculated the CoTs ε Syn , ε height , and ε speed from the three elements to investigate their contributions.

Simulation of Walking
First, we conducted a computer simulation of our neuromusculoskeletal model for the nominal speed of walking using β =β and = π/2 (see Supplementary Movie S1). The generated average speed was 0.2 m/s. Figures 7A,B show the joint angle and muscle activity, respectively, from the simulation compared with animal data. The simulation results show activity patterns similar to those of animals in terms of kinematics and muscle activity levels. However, our model was limited in its ability to accurately reproduce the locomotor behavior observed in animals. In particular, the elbow and knee joints were more extended than those of animals, which resulted in a shorter stride and slower speed than desired. The more extended elbow joint partly occurred because we did not incorporate the hand and wrist in the forelimb, and the forward speed was reduced by large ground reaction forces at the tips of the forelimbs. Similarly, the more extended knee joint partly occurred because we did not incorporate the phalangeal part in the hindlimb. The absence of flexibility of the spine in the trunk is another factor causing the extended posture. In addition, the activity of the SSP muscle appeared in a phase different from that of measured data. The SSP muscle in animals was activated in the same phase as that of the antagonistic SPD muscle so that the shoulder joint stiffness increased. In contrast, the SSP muscle in our model was activated in the same phase as that of the ipsilateral BR and BIC muscles.

Changing Gait and Speed
By changing from π/2 to π, our model achieved steady trotting (see Supplementary Movie S2). Although this gait had activity patterns almost identical to those of walking in terms of joint kinematics and muscle activations (Figure 7), the footfall pattern was different, as shown in Figure 8. The difference of the footfall pattern caused the difference in the trunk movement. In particular, while walking has a slight pitching movement of the trunk, trotting has almost no pitching movement, as shown in Figure 8 (see Supplementary Movies S1, S2).
To change the speed of each gait, we slowly increased or decreased β fromβ while changing the duration of the extension phase, the amplitude of the muscle-synergy-based motor commands, and the reference values for the movement regulator based on β, as in (12)(13)(14). Figure 9 shows the speed of the simulated walking and trotting. Our model achieved speeds of 0.15-0.2 m/s for walking and 0.18-0.22 m/s for trotting. Trotting was faster than walking in each β. Figure 10A shows the CoT ε of walking and trotting for the generated speeds in the simulation. Both CoT curves are Ushaped. The speeds for the minimum CoT for walking and trotting are very different. Walking had lower (higher) CoTs than trotting at slow (fast) speed. The CoT was obtained by dividing the mechanical work W for one gait cycle duration by the stride length D, as in (15). Figures 10B,C show the mechanical work and stride length, respectively, with speed. The mechanical work slightly but monotonically increased in walking and decreased in trotting. In contrast, the stride length shows a single-peaked shape for speed in both walking and trotting. The speeds for the minimum CoT and maximum stride length were almost identical in both walking and trotting. Figures 10D,E, respectively, show the contributions of the muscle synergy-based pulses and posture and speed regulators to the CoT (ε Syn , and ε height and ε speed , respectively). The CoT contribution of the muscle synergy-based pulses was U-shaped in both walking and trotting and was the largest among the three elements. The CoT contributions of the posture and speed regulators were small. While the contribution of the posture regulator remained constant with speed in walking, it decreased in trotting. The contribution of the speed regulator increased in both walking and trotting.

DISCUSSION
In this study, we improved our previous musculoskeletal model of rat hindlimbs (Aoi et al., 2013a) to construct a whole-body rat musculoskeletal model, which consists of the trunk, forelimbs, and hindlimbs. We also improved our motor control model (Aoi et al., 2013a) based on the muscle synergy hypothesis to control the whole-body rat model. Although the motor control model had a large number of motor control parameters, the rat model could be made to walk or trot by changing only the phase difference of the muscle-synergy-based motor commands between the forelimb and hindlimb (Figures 7, 8). Furthermore, the speed of each gait could be varied by changing only the duration of the extension phase, the amplitude of the musclesynergy-based motor commands, and the reference values for the movement regulator (Figure 9). The relation between speed and CoT was U-shaped for both the walking and trotting generated, and the speeds for the minimum CoT were different for the two gaits, as observed in the oxygen consumption of animals (Figure 10).

Characteristics of Cost of Transport
For our simulation, the CoT vs. speed curves were U-shaped for both walking and trotting ( Figure 10A). Walking had lower (higher) CoTs than trotting at slow (fast) speed. The CoTs were the same at a certain middle speed. These results indicate that  walking and trotting are energy-efficient at slow and fast speeds, respectively. These trends are similar to those observed for animals (Figure 1).
The CoT was calculated by dividing the mechanical work of one gait cycle duration by the stride length, as shown in (15). The stride length showed a single-peaked shape against speed ( Figure 10C). The speeds for the minimum CoT and maximum stride length were almost identical. We decreased the extension phase duration to increase the speed, which decreased the gait cycle duration. Because we increased the muscle-synergy-based motor commands during the extension phase as in (13), the stride length increased. However, this increase of the stride length was limited due to the increase of the gait frequency (decrease of the gait cycle duration). The stride length decreased over a critical frequency, which suggests a resonance property of the musculoskeletal dynamics and motor control input. Although these trends were similar between walking and trotting, the maximum stride length differed. These characteristics contributed majorly to the different energy efficiencies of gait. It has been reported that when the locomotor frequency increased, the stretch receptor of the hip prevented the hindlimbs from extending further (Mayer et al., 2018;Santuz et al., 2019). In the future, we would like to incorporate this sensory regulation model to control the stride length to investigate the mechanism of energy-efficient locomotion further.
In our model, the CoT had contributions from three elements, namely muscle synergy-based pulses and posture and speed regulators (ε ≃ ε Syn + ε height + ε speed ). The muscle synergy-based pulses had the largest contribution and determined the basic Ushaped characteristics (Figure 10D). Although the posture and speed regulators had small contributions (Figure 10E), they had specific characteristics. In particular, while the posture regulator for walking remains almost constant with speed, that for trotting decreased, which moved the speed for the minimum CoT to the right (with respect to that for the muscle synergy-based pulses) and increased the difference in speed for the minimum CoT between walking and trotting. This allowed the model to achieve energy-efficient locomotion in a wider speed range. In contrast, the speed regulator increased with speed in both walking and trotting and had a similar shape against speed for walking and trotting, which had a small contribution to the difference in speed for the minimum CoT.

Gait Generation Based on Muscle Synergy
A large portion of motor commands in our model was generated by a linear combination of four rectangular pulses for each limb, where we used the same onset and end phases for the pulses between the four limbs. We changed the relative phase of the pulses between the forelimb and hindlimb to make the gait generation simple. However, a muscle synergy analysis of dogs showed that although a large portion of the muscle activity can be reproduced by a linear combination of four basic patterns for both forelimbs and hindlimbs in walking and trotting, as done for our model, the basic patterns had some differences, especially in the activation timings between forelimbs and hindlimbs and between walking and trotting (Deban et al., 2012;Catavitello et al., 2015). In particular, the basic patterns for the late extension and early and late flexion of the hindlimbs were earlier than those of the forelimbs in walking. The basic pattern for the early flexion of the hindlimbs was earlier than that of the forelimbs in trotting. The control of the activation timings of the muscle synergy patterns could contribute to the gait change (Cappellini et al., 2006;Aoi et al., 2019). In future studies, we would like to measure the trotting of rats and incorporate motor control differences between forelimbs and hindlimbs and between walking and trotting to clarify the gait-generation mechanism further.

Limitations of Our Model and Future Work
Although our simulation results showed features similar to those of animals (Figure 7), our model has limitations that prevent it from accurately reproducing animal locomotion. In particular, we did not incorporate the hand, phalangeal part of the hindlimbs, or flexibility of the spine in the trunk (Schilling and Hackert, 2006). These elements might improve the gait speed and energy efficiency. Furthermore, we confined our musculoskeletal model to two dimensions, which neglected instability in the lateral direction. More contribution of the posture regulator would be required for a three-dimensional model to maintain a stable posture during locomotion. In addition, although head movements are important and specific for gait (Zsoldos et al., 2010), we did not incorporate the neck. We would like to incorporate these features to clarify adaptive motor control mechanisms in animal locomotion further.
Not only the metabolic cost of locomotion but also other factors, such as musculoskeletal forces (Farley and Taylor, 1991), gait stability (Schöner et al., 1990;Diedrich and Warren, 1995;Aoi et al., 2013b), terrain and ground surface conditions (Prost and Sussman, 1969;Gustås et al., 2006;Goldenberg et al., 2008;Chateau et al., 2013), and genetic mutation (Andersson et al., 2012), influence the gait decision of animals. Furthermore, although animals change their gait smoothly when triggered by these factors, the gait transition mechanism also remains unclear.
Our neuromusculoskeletal model will be useful for investigating these mechanisms in the future.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

ETHICS STATEMENT
The animal study was reviewed and approved by the Ethical Committee for Animal Experiments at the University of Tokyo.