Impact Factor 4.677 | CiteScore 5.4
More on impact ›

ORIGINAL RESEARCH article

Front. Neurosci., 17 January 2020 | https://doi.org/10.3389/fnins.2019.01337

Gait Generation and Its Energy Efficiency Based on Rat Neuromusculoskeletal Model

  • 1Department of Life Sciences, Graduate School of Arts and Sciences, The University of Tokyo, Tokyo, Japan
  • 2Department of Aeronautics and Astronautics, Graduate School of Engineering, Kyoto University, Kyoto, Japan
  • 3Department of Physiology and Biological Information, School of Medicine, Dokkyo Medical University, Tochigi, Japan
  • 4Department of Mechanical Engineering and Intelligent Systems, Graduate School of Informatics and Engineering, The University of Electro-Communications, Tokyo, Japan

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.

1. 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.

FIGURE 1
www.frontiersin.org

Figure 1. Oxygen consumption of horses in walking, trotting, and galloping.

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., 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.

2. Method

2.1. 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.

FIGURE 2
www.frontiersin.org

Figure 2. Musculoskeletal model.

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 Fm (m = SSP, SPD, BR, TRIL, BIC, TRI, IP, GM, VL, TA, SO, BF, and GA) through contractile and passive elements, which is given based on Aoi et al. (2013a) by

Fm=Fmmax(amFmlFmv+Fmp)    (1)

where Fmmax is the maximum muscle tension, am is the muscle activation (0 ≤ am ≤ 1), Fml is the force-length relationship, Fmv is the force-velocity relationship, and Fmp is the passive component. The muscle lengths were normalized by lmmax, which was set so that all uniarticular muscles had a length of 85% of lmmax and all biarticular muscles had a length of 75% of lmmax 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.8lmmax.

The muscle activation am is determined through

τacta˙m+{τactτdeact+(1-τactτdeact)um}am=um    (2)

where τact and τdeact are respectively, activation and deactivation time constants (11 and 18 ms, respectively) and um is the motor command determined in the motor control model.

2.2. 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 um is the summation of the two components from the movement generator and the movement regulator, namely umSyn and umReg, respectively.

um=umSyn+umReg    (3)

2.2.1. Movement Generator

The movement generator is based on the hypothetical two-layer 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 ϕij (i = left, right, j = fore, hind) for the oscillator phase (0ϕij<2π), which follows the dynamics given by

ϕ˙leftfore=2πT-K1 sin(ϕleftfore-ϕrightfore-π)-K2 sin(ϕleftfore-ϕlefthind+Δ)ϕ˙rightfore=2πT-K1 sin(ϕrightfore-ϕleftfore-π)-K2 sin(ϕrightfore-ϕrighthind+Δ)ϕ˙lefthind=2πT-K1 sin(ϕlefthind-ϕrighthind-π)-K2 sin(ϕlefthind-ϕleftfore-Δ)ϕ˙righthind=2πT-K1 sin(ϕrighthind-ϕlefthind-π)-K2 sin(ϕrighthind-ϕrightfore-Δ)    (4)

where T is the gait cycle duration and K1 and K2 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 (Ivanenko et al., 2006; Dominici et al., 2011; Markin et al., 2012; Catavitello et al., 2015; Rigosa et al., 2015). Specifically, we used four rectangular pulses pi (i = 1, …, 4) for each limb, which are given by

pi(ϕ)={1Φiϕ<Ψi0otherwise  i=1,,4    (5)

where Φi and Ψi (i = 1, …, 4) are the onset and end phases of the pulse, respectively, and we omitted the suffix of ϕ. p1, p2, p3, and p4 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 umSyn of the movement generator is given by

umSyn=i=14wm,ipi(ϕ)    (6)

where wm, i (i = 1, …, 4) is the weighting coefficient.

FIGURE 3
www.frontiersin.org

Figure 3. Movement generator. (A) Motor commands as linear combinations of four rectangular pulses in each forelimb and hindlimb. Gray regions indicate extension phases, and others indicate flexion phases. (B) Activated muscles at each pulse.

2.2.2. 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.

FIGURE 4
www.frontiersin.org

Figure 4. Movement regulator based on shoulder height, hip height, and forward speed. BR and TRIL muscles are used for shoulder height, VL, TA, and SO muscles are used for hip height, and SSP, SPD, IP, GM, TA, and SO muscles are used for speed.

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 pmheight (m = BR and TRIL) is given by

pmheight={-Kmheight(hShoulder-h0Shoulder)-Dmheighth˙Shoulderin stance phase0otherwise    (7)

where hShoulder and h˙Shoulder are the shoulder height and its rate, respectively, h0Shoulder is the reference height, and Kmheight and Dmheight are gain parameters. For the hindlimbs, we used the VL, TA, and SO muscles to maintain the hip height. The motor command pmheight (m = VL, TA, and SO) is given by

pmheight={-Kmheight(hHip-h0Hip)-Dmheighth˙Hipin stance phase0otherwise    (8)

where hHip and h˙Hip are the hip height and its rate, respectively, and h0Hip 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 pmspeed (m = SSP, SPD, IP, GM, TA, and SO) is given by

pmspeed={-Kmspeed(v-v0)in stance phase0otherwise    (9)

where v is the forward speed, v0 is its desired value, and Kmspeed 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 umReg of the movement regulator is given by

umReg(t)=umheight(t)+umspeed(t)    (10)

where

umheight(t)=pmheight(t-τDelay)umspeed(t)=pmspeed(t-τDelay)    (11)

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.

2.3. 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 muscle-synergy-based motor command umSyn between the ipsilateral limbs by changing Δ in (4). In particular, we used Δ = π/2 for walking and Δ = π for trotting.

FIGURE 5
www.frontiersin.org

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).

Animals change the gait cycle duration to vary speed, where the duration of the flexion phase for swinging the limb 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 Tex (= βT) using β while keeping the duration of the flexion phase Tfl unchanged (T = Tex + Tfl = Tfl/(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, …, 4), h^0Shoulder, h^0Hip, and v^0 for motor control parameters β, Tex, Φi, Ψi, wm, i (i = 1, …, 4), h0Shoulder, h0Hip, and v0, respectively. The onset phase Φi and end phase Ψi (i = 1, …, 4) of each pulse are given by

Φi={ββ^Φ^ii=1,21-β1-β^Φ^i+2π(β-β^)1-β^i=3,4Ψi={ββ^Ψ^ii=1,21-β1-β^Ψ^i+2π(β-β^)1-β^i=3,4    (12)

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 wm, i (i = 1, 2) of the muscle-synergy-based rectangular pulses for the extension phase (Figure 6B) as

wm,i=1-β1-β^w^m,ii=1,2    (13)

As we changed the locomotion speed, we also changed the reference height (shoulder, hip) and speed for the movement regulator (Figures 6C–E) as

h0Shoulder=h^0Shoulder+αShoulder(β-β^)h0Hip=h^0Hip+αHip(β-β^)v0=v^0+αSpeed(β-β^)    (14)

where αShoulder, αHip, and αSpeed are coefficients.

FIGURE 6
www.frontiersin.org

Figure 6. Regulation of muscle-synergy-based motor command (A) to change speed by changing extension phase duration Tex (= βT) while keeping flexion phase duration Tfl unchanged, where weighting coefficient wm, i/ŵm, i (B), reference shoulder height h0Shoulder (C), reference hip height h0Hip (D), and reference speed v0 (E) depend on β.

2.4. Model Parameters

2.4.1. 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 Fmmax was determined based on the measured PCSA and the moment arms were determined from the center of the range of joint movement during locomotion.

TABLE 1
www.frontiersin.org

Table 1. Physical parameters of skeletal model.

TABLE 2
www.frontiersin.org

Table 2. Physical parameters of muscle model.

2.4.2. Parameters for the Motor Control Model

Based on measured data for rats walking on a treadmill at a speed of 0.4 m/s (Aoi et al., 2013a), we set the durations of the flexion and extension phases for the nominal speed as Tfl = 0.10 s and T^ex=0.16 s, respectively (β^=0.62). We determined the motor control parameters for the nominal speed as follows so that our model achieved steady walking based on our previous results of the hindlimb model (Aoi et al., 2013a): Φ^1=0, Φ^2=0.40π, Φ^3=1.24π, Φ^4=1.42π, Ψ^1=0.33π, Ψ^2=0.89π, Ψ^i=1.42π, Ψ^i=1.71π, ŵSSP, 1 = 0.24, ŵSSP, 4 = 0.20, ŵSPD, 2 = 0.27, ŵSPD, 3 = 0.08, ŵBR, 3 = 0.09, ŵTRIL, 1 = 0.47, ŵTRIL, 2 = 0.57, ŵBIC, 3 = 0.17, ŵBIC, 4 = 0.08, ŵTRI, 1 = 0.27, ŵTRI, 2 = 0.56, ŵIP, 3 = 0.32, ŵIP, 4 = 0.32, ŵGM, 1 = 0.61, ŵGM, 2 = 0.25, ŵVL, 1 = 0.19, ŵVL, 2 = 0.22, ŵTA, 3 = 0.45, ŵTA, 4 = 0.06, ŵSO, 1 = 0.58, ŵSO, 2 = 0.14, ŵBF, 1 = 0.22, ŵBF, 2 = 0.12, ŵBF, 3 = 0.09, ŵGA, 1 = 0.47, ŵGA, 2 = 0.10, ŵm, i = 0 for the other values of m and i, h^0Shoulder=0.033 m, h^0Hip=0.054 m, v^0=0.4 m/s, KBRheight=-2.07, KTRILheight=2.07, KVLheight=12.4, KTAheight=-12.4, KSOheight=12.4, DBRheight=-0.001, DTRILheight=0.001, DVLheight=0.006, DTAheight=-0.006, DSOheight=0.006, KSSPspeed=-0.007, KSPDspeed=0.007, KIPspeed=-0.052, KGMspeed=0.052, KTAspeed=-0.026, KSOspeed=0.026, K1 = 20, and K2 = 10. In addition, we set the coefficients for the regulation of the references in the movement regulator to change the speed as αShoulder = 0.01 m, αHip = 0.01 m, and αSpeed = −4.7 m/s.

2.5. 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.

2.6. 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

ε=WD    (15)

where

W=η++η-η+=TmFm[vm]+dtη-=14TmFm[-vm]+dt

vm 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 um is generated by three elements: rectangular pulses umSyn in the movement generator and motor commands umheight and umspeed to regulate the posture and speed, respectively, in the movement regulator (um=umSyn+umheight+umspeed). Because they determined the muscle activation am in (2), we calculated amSyn, amheight, and amspeed from umSyn, umheight, and umspeed, respectively. Using these values, we calculated the CoTs εSyn, εheight, and εspeed from the three elements to investigate their contributions.

3. Results

3.1. 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.

FIGURE 7
www.frontiersin.org

Figure 7. Simulated joint kinematics (A) and muscle activations (B) in walking compared with measured animal data. Gray regions indicate stance phases.

3.2. 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).

FIGURE 8
www.frontiersin.org

Figure 8. Footfall pattern and stick diagram of simulated walking (A) and trotting (B). Stick diagram shows simulated locomotor behavior between two successive foot contacts of right hindlimb, where bold lines indicate right limbs. See Supplementary Movies S1, S2 for simulated walking and trotting, respectively.

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–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 9
www.frontiersin.org

Figure 9. Generated speed during walking and trotting.

3.3. Cost of Transport

Figure 10A shows the CoT ε of walking and trotting for the generated speeds in the simulation. Both CoT curves are U-shaped. 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.

FIGURE 10
www.frontiersin.org

Figure 10. Cost of transport for speed of simulated walking and trotting. (A) Total CoT, (B) mechanical work, (C) stride length, (D) CoT contribution of muscle synergy-based pulses, and (E) CoT contributions of posture and speed regulators.

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.

4. 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 muscle-synergy-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).

4.1. 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 U-shaped 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.

4.2. 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.

4.3. 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.

Author Contributions

SA developed the study design. MT performed the computer simulations and analyzed the data in consultation with SA, SF, TF, KT, and DY. MT and SA wrote the manuscript, and all the authors reviewed and approved it.

Funding

This study was supported in part by JSPS KAKENHI Grant-in-Aid for Scientific Research (B) JP15KT0015 and Grant-in-Aid for Scientific Research on Innovative Areas JP26120006.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2019.01337/full#supplementary-material

Supplementary Movie S1. Simulated walking.

Supplementary Movie S2. Simulated trotting.

References

Andersson, L. S., Larhammar, M., Memic, F., Wootz, H., Schwochow, D., Rubin, C. J., et al. (2012). Mutations in DMRT3 affect locomotion in horses and spinal circuit function in mice. Nature 488, 642–646. doi: 10.1038/nature11399

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoi, S., and Funato, T. (2016). Neuromusculoskeletal models based on the muscle synergy hypothesis for the investigation of adaptive motor control in locomotion via sensory-motor coordination. Neurosci. Res. 104, 88–95. doi: 10.1016/j.neures.2015.11.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoi, S., Katayama, D., Fujiki, S., Tomita, N., Funato, T., Yamashita, T., et al. (2013b). A stability-based mechanism for hysteresis in the walk-trot transition in quadruped locomotion. J. R. Soc. Interface 10:20120908. doi: 10.1098/rsif.2012.0908

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoi, S., Kondo, T., Hayashi, N., Yanagihara, D., Aoki, S., Yamaura, H., et al. (2013a). Contributions of phase resetting and interlimb coordination to the adaptive control of hindlimb obstacle avoidance during locomotion in rats: a simulation study. Biol. Cybern. 107, 201–216. doi: 10.1007/s00422-013-0546-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoi, S., Ohashi, T., Bamba, R., Fujiki, S., Tamura, D., Funato, T., et al. (2019). Neuromusculoskeletal model that walks and runs across a speed range with a few motor control parameter changes based on the muscle synergy hypothesis. Sci. Rep. 9:369. doi: 10.1038/s41598-018-37460-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Aoki, S., Sato, Y., and Yanagihara, D. (2013). Lesion in the lateral cerebellum specifically produces overshooting of the toe trajectory in leading forelimb during obstacle avoidance in the rat. J. Neurophysiol. 110, 1511–1524. doi: 10.1152/jn.01048.2012

PubMed Abstract | CrossRef Full Text | Google Scholar

Burke, R. E., Degtyarenko, A. M., and Simon, E. S. (2001). Patterns of locomotor drive to motoneurons and last-order interneurons: clues to the structure of the CPG. J. Neurophysiol. 86, 447–462. doi: 10.1152/jn.2001.86.1.447

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappellini, G., Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2006). Motor patterns in human walking and running. J. Neurophysiol. 95, 3426–3437. doi: 10.1152/jn.00081.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Catavitello, G., Ivanenko, Y. P., and Lacquaniti, F. (2015). Planar covariation of hindlimb and forelimb elevation angles during terrestrial and aquatic locomotion of dogs. PLoS ONE 10:e0133936. doi: 10.1371/journal.pone.0133936

PubMed Abstract | CrossRef Full Text | Google Scholar

Chateau, H., Camus, M., Holden-Douilly, L., Falala, S., Ravary, B., Vergari, C., et al. (2013). Kinetics of the forelimb in horses circling on different ground surfaces at the trot. Vet. J. 198, e20–e26. doi: 10.1016/j.tvjl.2013.09.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Clarke, K., and Still, J. (1999). Gait analysis in the mouse. Physiol. Behav. 66, 723–729. doi: 10.1016/S0031-9384(98)00343-6

PubMed Abstract | CrossRef Full Text | Google Scholar

d'Avella, A., Saltiel, P., and Bizzi, E. (2003). Combinations of muscle synergies in the construction of a natural motor behavior. Nat. Neurosci. 6, 300–308. doi: 10.1038/nn1010

PubMed Abstract | CrossRef Full Text | Google Scholar

Deban, S. M., Schilling, N., and Carrier, D. R. (2012). Activity of extrinsic limb muscles in dogs at walk, trot and gallop. J. Exp. Biol. 215, 287–300. doi: 10.1242/jeb.063230

PubMed Abstract | CrossRef Full Text | Google Scholar

Diedrich, F. J., and Warren, W. H. Jr. (1995). Why change gaits? Dynamics of the walk-run transition. J. Exp. Psychol. Hum. Percept. Perform. 21, 183–202. doi: 10.1037/0096-1523.21.1.183

PubMed Abstract | CrossRef Full Text | Google Scholar

Dominici, N., Ivanenko, Y. P., Cappellini, G., d'Avella, A., Mondì, V., Cicchese, M., et al. (2011). Locomotor primitives in newborn babies and their development. Science 334, 997–999. doi: 10.1126/science.1210617

PubMed Abstract | CrossRef Full Text | Google Scholar

Drew, T., Kalaska, J., and Krouchev, N. (2008). Muscle synergies during locomotion in the cat: a model for motor cortex control. J. Physiol. 586, 1239–1245. doi: 10.1113/jphysiol.2007.146605

PubMed Abstract | CrossRef Full Text | Google Scholar

Ekeberg, Ö, and Pearson, K. (2005). Computer simulation stepping in the hind legs of the cat: An examination of mechanisms regulating the stance-to-swing transition. J. Neurophysiol. 94, 4256–4268. doi: 10.1152/jn.00065.2005

PubMed Abstract | CrossRef Full Text | Google Scholar

Elmer, S. J., and LaStayo, P. C. (2014). Revisiting the positive aspects of negative work. J. Exp. Biol. 217, 2434–2436. doi: 10.1242/jeb.092247

PubMed Abstract | CrossRef Full Text | Google Scholar

Farley, C. T., and Taylor, C. R. (1991). A mechanical trigger for the trot-gallop transition in horses. Science 253, 306–308. doi: 10.1126/science.1857965

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujiki, S., Aoi, S., Funato, T., Sato, Y., Tsuchiya, K., and Yanagihara, D. (2018). Adaptive hindlimb split-belt treadmill walking in rats by controlling basic muscle activation patterns via phase resetting. Sci. Rep. 8:17341. doi: 10.1038/s41598-018-35714-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Fukuoka, Y., Habu, Y., and Fukui, T. (2015). A simple rule for quadrupedal gait generation determined by leg loading feedback: a modeling study. Sci. Rep. 5:8169. doi: 10.1038/srep08169

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldenberg, F., Glanzl, M., Henschel, J. R., Funk, S. M., and Millesi, E. (2008). Gait choice in desert-living black-backed jackals. J. Zool. 275, 124–129. doi: 10.1111/j.1469-7998.2008.00417.x

CrossRef Full Text | Google Scholar

Górska, T., Zmysłowski, W., and Majczyński, H. (1999). Overground locomotion in intact rats: Interlimb coordination, support patterns and support phases duration. Acta Neurobiol. Exp. 59, 131–144.

PubMed Abstract | Google Scholar

Goslow, G. E. Jr., Reinking, R. M., and Stuart, D. G. (1973). The cat step cycle: hind limb joint angles and muscle lengths during unrestrained locomotion. J. Morph. 141, 1–41. doi: 10.1002/jmor.1051410102

PubMed Abstract | CrossRef Full Text | Google Scholar

Gustås, P., Johnston, C., and Drevemo, S. (2006). Ground reaction force and hoof deceleration patterns on two different surfaces at the trot. Equine Comp. Exerc. Physiol. 3, 209–216. doi: 10.1017/S147806150667607X

CrossRef Full Text | Google Scholar

Heglund, N. C., and Taylor, C. R. (1998). Speed, stride frequency and energy cost per stride: how do they change with body size and gait? J. Exp. Biol. 138, 301–318.

PubMed Abstract | Google Scholar

Hoyt, D. F., and Taylor, C. R. (1981). Gait and the energetics of locomotion in horses. Nature 292, 239-240. doi: 10.1038/292239a0

CrossRef Full Text | Google Scholar

Hunt, A., Schmidt, M., Fischer, M., and Quinn, R. (2015). A biologically based neural system coordinates the joints and legs of a tetrapod. Bioinspir. Biomim. 10:055004. doi: 10.1088/1748-3190/10/5/055004

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivanenko, Y. P., Poppele, R. E., and Lacquaniti, F. (2006). Motor control programs and walking. Neuroscientist 12, 339–348. doi: 10.1177/1073858406287987

PubMed Abstract | CrossRef Full Text | Google Scholar

Ivashko, D. G., Prilutsky, B. I., Markin, S. N., Chapin, J. K., and Rybak, I. A. (2003). Modeling the spinal cord neural circuitry controlling cat hindlimb movement during locomotion. Neurocomputing 52–54, 621–629. doi: 10.1016/S0925-2312(02)00832-9

CrossRef Full Text | Google Scholar

Margaria, R., Cerretelli, P., Aghemo, P., and Sassi, G. (1963). Energy cost of running. J. Appl. Physiol. 18, 367–370. doi: 10.1152/jappl.1963.18.2.367

PubMed Abstract | CrossRef Full Text | Google Scholar

Markin, S. N., Klishko, A. N., Shevtsova, N. A., Lemay, M. A., Prilutsky, B. I., and Rybak, I. A. (2016). “A neuromechanical model of spinal control of locomotion,” in Neuromechanical Modeling of Posture and Locomotion, eds B. I. Prilutsky and D. H. Edwards (New York, NY: Springer Science), 21–65.

Google Scholar

Markin, S. N., Lemay, M. A., Prilutsky, B. I., and Rybak, I. A. (2012). Motoneuronal and muscle synergies involved in cat hindlimb control during fictive and real locomotion: a comparison study. J. Neurophysiol. 107, 2057–2071. doi: 10.1152/jn.00865.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayer, W. P., Murray, A. J., Brenner-Morton, S., Jessell, T. M., Tourtellotte, W. G., and Akay, T. (2018). Role of muscle spindle feedback in regulating muscle activity strength during walking at different speed in mice. J. Neurophysiol. 120, 2484–2497. doi: 10.1152/jn.00250.2018

PubMed Abstract | CrossRef Full Text | Google Scholar

Nishii, J. (2006). An analytical estimation of the energy cost for legged locomotion. J. Theor. Biol. 238, 636–645. doi: 10.1016/j.jtbi.2005.06.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Ogihara, N., Aoi, S., Sugimoto, Y., Tsuchiya, K., and Nakatsukasa, M. (2011). Forward dynamic simulation of bipedal walking in the Japanese Macaque: investigation of causal relationships among limb kinematics, speed, and energetics of bipedal locomotion in a nonhuman primate. Am. J. Phys. Anthropol. 145, 568–580. doi: 10.1002/ajpa.21537

PubMed Abstract | CrossRef Full Text | Google Scholar

Prost, J. H., and Sussman, R. W. (1969). Monkey locomotion on inclined surfaces. Am. J. Phys. Anthropol. 31, 53–58. doi: 10.1002/ajpa.1330310107

CrossRef Full Text | Google Scholar

Rigosa, J., Panarese, A., Dominici, N., Friedli, L., van den Brand, R., Carpaneto, J., et al. (2015). Decoding bipedal locomotion from the rat sensorimotor cortex. J. Neural Eng. 12:056014. doi: 10.1088/1741-2560/12/5/056014

PubMed Abstract | CrossRef Full Text | Google Scholar

Rybak, I. A., Shevtsova, N. A., Lafreniere-Roula, M., and McCrea, D. A. (2006). Modelling spinal circuitry involved in locomotor pattern generation: insights from deletions during fictive locomotion. J. Physiol. 577, 617–639. doi: 10.1113/jphysiol.2006.118703

PubMed Abstract | CrossRef Full Text | Google Scholar

Santuz, A., Akay, T., Mayer, W. P., Wells, T. L., Schroll, A., and Arampatzis, A. (2019). Modular organization of murine locomotor pattern in the presence and absence of sensory feedback from muscle spindles. J. Physiol. 597, 3147–3165. doi: 10.1113/JP277515

PubMed Abstract | CrossRef Full Text | Google Scholar

Schilling, N., and Hackert, R. (2006). Sagittal spine movements of small therian mammals during asymmetrical gaits. J. Exp. Biol. 209, 3925–3939. doi: 10.1242/jeb.02400

PubMed Abstract | CrossRef Full Text | Google Scholar

Schöner, G., Jiang, W. Y., and Kelso, J. A. S. (1990). A synergetic theory of quadrupedal gaits and gait transitions. J. Theor. Biol. 142, 359–391. doi: 10.1016/S0022-5193(05)80558-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Takakusaki, K. (2017). Functional neuroanatomy for posture and gait control. J. Mov. Disord. 10, 1–17. doi: 10.14802/jmd.16062

PubMed Abstract | CrossRef Full Text | Google Scholar

Takei, T., Confaisa, J., Tomatsua, S., Oya, T., and Seki, K. (2017). Neural basis for hand muscle synergies in the primate spinal cord. Proc. Natl. Acad. Sci. U.S.A. 114, 8643–8648. doi: 10.1073/pnas.1704328114

PubMed Abstract | CrossRef Full Text | Google Scholar

Ting, L. H., and Macpherson, J. M. (2005). A limited set of muscle synergies for force control during a postural task. J. Neurophysiol. 93, 609–613. doi: 10.1152/jn.00681.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Todorov, E., and Jordan, M. I. (2002). Optimal feedback control as a theory of motor coordination. Nat. Neurosci. 5, 1226–1235. doi: 10.1038/nn963

PubMed Abstract | CrossRef Full Text | Google Scholar

Tresch, M. C., Saltiel, P., and Bizzi, E. (1999). The construction of movement by the spinal cord. Nat. Neurosci. 2, 162–167. doi: 10.1038/5721

PubMed Abstract | CrossRef Full Text | Google Scholar

Yakovenko, S., Gritsenko, V., and Prochazka, A. (2004). Contribution of stretch reflexes to locomotor control: a modeling study. Biol. Cybern. 90, 146–155. doi: 10.1007/s00422-003-0449-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yakovenko, S., McCrea, D. A., Stecina, K., and Prochazka, A. (2005). Control of locomotor cycle durations. J. Neurophysiol. 94, 1057–1065. doi: 10.1152/jn.00991.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Zsoldos, R. R., Kotschwar, A. B., Kotschwar, A., Groesel, M., Licka, T., and Peham, C. (2010). Electromyography activity of the equine splenius muscle and neck kinematics during walk and trot on the treadmill. Equine Vet. J. Suppl. 38, 455–461. doi: 10.1111/j.2042-3306.2010.00263.x

CrossRef Full Text | Google Scholar

Keywords: rat, walk, trot, energy efficiency, central pattern generator, muscle synergy, neuromusculoskeletal model

Citation: Toeda M, Aoi S, Fujiki S, Funato T, Tsuchiya K and Yanagihara D (2020) Gait Generation and Its Energy Efficiency Based on Rat Neuromusculoskeletal Model. Front. Neurosci. 13:1337. doi: 10.3389/fnins.2019.01337

Received: 29 June 2019; Accepted: 27 November 2019;
Published: 17 January 2020.

Edited by:

Jun Izawa, University of Tsukuba, Japan

Reviewed by:

Srinivasa Chakravarthy, Indian Institute of Technology Madras, India
Gabriella Lindgren, Swedish University of Agricultural Sciences, Sweden

Copyright © 2020 Toeda, Aoi, Fujiki, Funato, Tsuchiya and Yanagihara. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Shinya Aoi, shinya_aoi@kuaero.kyoto-u.ac.jp; Dai Yanagihara, dai-y@idaten.c.u-tokyo.ac.jp