Based on the results of previous studies (Ogihara and Yamazaki, 2001; Aoi et al., 2010; Ichimura and Yamazaki, 2019), we constructed a two-dimensional musculoskeletal model including the head, arms, torso (HAT), thighs, shanks, and feet (Figure 1). We determined the segment length and inertia parameters of the model based on past findings (Jo and Massaquoi, 2007; Ichimura and Yamazaki, 2019). Each joint was modeled as a pin joint with a linear viscous element. The viscosity coefficients of the hip, knee, and ankle joints were 1.09, 3.17, and 0.943 Nms/rad, respectively (Aoi et al., 2010; Ichimura and Yamazaki, 2019). The knee and ankle joints were locked to avoid hyperextension and hyperflexion which are unrealistic joint ranges of motion. The heel or toe received ground reaction forces (GRF) when they contacted the ground. The GRF is modeled by a linear spring and damper system which was employed in previous studies (Ogihara and Yamazaki, 2001; Aoi et al., 2010; Ichimura and Yamazaki, 2019) and could mimic the measured GRF. The coefficients of the spring and damper were 5.0 × 10³ N/m and 1.0 × 10³ Ns/m in the horizontal direction and 2.5 × 10⁴ N/m and 5.0 × 10² Ns/m in the vertical direction, respectively (Ichimura and Yamazaki, 2019). Nine primary muscles were used in each leg for the muscle model (Figure 1): gluteus maximus (GM), iliopsoas (IL), biceps femoris long head (BFL), rectus femoris (RF), biceps femoris short head (BFS), vastus (VA), gastrocnemius (GC), soleus (SO), and tibialis anterior (TA). Muscles receive signals from the corresponding α-motoneurons and generate muscle tension through force-length and force-velocity relationships. We used the following mathematical model described by Ogihara and Yamazaki (2001), which included contractile (CE), passive elastic (PE), and passive damping (PD) elements, respectively: F m = F ¯ m CE ⋅ k ξ m ⋅ h η m ⋅ α m + F m PD + F m PE , k ξ m = 0.32 + 0.71 ⁡ exp − 1.112 ξ m − 1.0 sin 3.722 ξ m − 0.656 , h η m = 1 + tanh 3.0 η m , F m PD = c m PD L ˙ m , F m PE = k m PE exp 15 L m − L ¯ m − 1.0 , (1) where F m is the muscle tension generated by the mth muscle, F ¯ m CE is the maximum muscle tension due to the CE, k ξ m is the force-length relationship, h η m is the force-velocity relationship, α m is the stimulus signal from the corresponding α-motoneuron (0 ≤ α m ≤ 1), and F m PD and F m PE are the forces generated by the damping and elastic elements, respectively. ξ m and η m are the normalized muscle length and contraction velocity divided by the muscle optimum length L ¯ m and the muscle maximum contraction velocity L ˙ ¯ m , that is, where ξ m = L m / L ¯ m , η m = L ˙ m / L ˙ ¯ m , L m , and L ˙ ¯ m are the muscle length and contraction velocity, respectively. c m PD is the viscosity coefficient, and k m PE is the coefficient of the elastic element. These parameters were also used by Ogihara and Yamazaki (2001), Aoi et al. (2010), and Ichimura and Yamazaki (2019).

Generally, locomotion is considered to be generated and induced by a rhythmic neural network in the spinal cord, called CPG (Grillner and Zangger, 1975). In the present study, we used a mathematical model of the CPG (Matsuoka, 1985; Taga et al., 1991) as follows: τ i u i ˙ = − u i + ∑ j = 1 12 w i   j CPG y j − β v i + u 0 + Feed i θ l s e g } l , { G R F s } s   |   { w k F e e d } k , τ i ′ v ˙ i = − v i + y i , y i = max 0 , u i , (2) where u i is the internal state of the ith neuron, and v i is a variable representing the self-inhibitory effect of the ith neuron. τ i and τ i ′ are time constants, β is a coefficient, and w i   j CPG is a connecting weight from the jth neuron to the ith neuron. u 0 is an external input, and Feed i is the feedback signal from the musculoskeletal system. θ l s e g is the segment angle (l ∈ HAT , thigh , shank , foot for each leg), G R F s is the vertical GRF ( s ∈ left   limb , right   limb ), and w k F e e d is the weight coefficient ( k = 1 , ⋯ , 32 ). The parameter values are listed in the Supplementary Material S1. y i also gives motor commands to α-motoneurons, which activate the muscles. The α-motoneurons also receive feedback signals from various reflexes, such as postural control. The α-motoneuron output α m is given as follows: α m = 2.0 1.0 + exp 0.25 ∑ i = 1 18 w m   i α y i + P m θ j } j , { θ l seg } l , { GRF s } s   |   { w o POS } o − 1.0 , (3) where w m   i α and w o POS are the weight coefficient ( o = 1 , ⋯ , 23 ), P m is a posture control, and θ j is the joint angle (j ∈ { hip , knee , ankle ). These parameter values are also listed in the Supplementary Material S1.

Our model has 51 free parameters ( u 0 ; w k F e e d ; w i   j CPG ; w o POS ) required to be determined to achieve stable locomotion. To search for these parameters, we employed standard genetic algorithms (GAs) (Ogihara and Yamazaki, 2001; Ichimura and Yamazaki, 2019; Oku et al., 2021). We used the evaluation function J to maximize, which is given by the following equation. J = D + P D < 3 m , D + P + 60 C D ≥ 3 m , (4) C = 1.0 TMV ∫ t = 0 T ∑ m = 0 18 E ˙ m d t , (5) where D is the distance until the model falls, P is the penalty given when the model falls, which is set to −3.0, and C is the gross metabolic cost of transport (Koelewijn et al., 2019). T , M, and V represent the locomotion duration, body mass, and walking speed, respectively. E ˙ m is the metabolic energy consumption by all muscles (Minetti and Alexander, 1997).

Initially, 51 free parameters were optimized using GAs to simulate normal locomotion that assumed a kinematically symmetric gait with a symmetric body (Ichimura and Yamazaki, 2019). After the model acquired normal locomotion, we simulated the pathological UTTA conditions. Specifically, the muscles of the unilateral lower leg (TA, SO, and GC) were removed. Subsequently, the weight of the lower leg, moment of inertia, and passive moment of the ankle joint were changed to 65%, 40%, and 400 Nm/rad, respectively, to mimic the lower limb prosthesis based on the findings of a previous study (Russell Esposito and Miller, 2018). The model failed to walk even for one step immediately after this manipulation. Then, we investigated two adaptation rule scenarios for this model. These scenarios optimized the parameters of the amplitude in the CPG signals and the feedback from the musculoskeletal model to the CPG model. The first scenario involved re-searching u 0 of the CPG model and the 16 symmetric feedback parameters w k F e e d for both lower limbs, which we called the ‘symmetric control model.’ This control strategy of locomotion was the same as that of the normal model. The other scenario involved re-searching u 0 of the CPG model and the 32 asymmetric feedback parameters w k F e e d for both lower limbs, which we called the ‘asymmetric control model.’ Such a control strategy for locomotion was different from that of the normal model. These scenarios assumed minimally adaptive locomotion based on the finding that the spinal cord network adapts dynamically (Rossignol et al., 2006). Finally, we analyzed the metabolic cost of transport (Eq. 5), as well as the absolute symmetry index (ASI), an indicator of the asymmetry of the spatiotemporal gait factor (Nolan et al., 2003; Bosch and Rosenbaum, 2010), as a qualitative assessment of locomotion using five steps (from the 3rd to the 8th step). The ASI was calculated by the following equation: ASI = 2.0 R − L R + L × 100 , (6) where R and L are stance time or step length of the right and left limb, respectively.

We implemented GAs utilizing Message Passing Interface (MPI), which is a library for parallel computing. We wrote all programs in the C language and used the fourth-order Runge-Kutta method for the numerical solution of the differential equations. The time step size was set to 0.1 ms.

We performed five simulations with five different random number generator seeds to verify that each simulation’s results were unique. Subsequently, we observed stable bipedal locomotion (walking for continuous 10-s periods), showing no qualitative differences in locomotion patterns due to differences in the seeds.

The model acquired a stable locomotion after 2000 generations of GAs (Figure 2). In Figure 2A, the locomotion pattern of the model qualitatively resembles that of human bipedal locomotion. The locomotion speed of the model was 0.88 m/s. Figures 2B–D show GRFs, joint angles, and muscle activities, respectively. To evaluate the validity of the current results, the cosine similarity (S) and correlation coefficient (R) between simulation data and measured data were calculated (Bovi et al., 2011). Notably, the muscle activations of the IL and BFS were not compared with simulation data due to the lack of measurement data (Figure 2D). However, other studies demonstrated that IL activity was mainly detected in the middle of the gait cycle, whereas BFS activity mainly occurred at the end of the gait cycle, which was qualitatively consistent with the simulation results (Perry and Burnfield, 2010). Thus, normal gait in this model was validated by results that were qualitatively and quantitatively comparable with those of previous studies (Aoi et al., 2010; Song and Geyer, 2015; Aoi et al., 2019). However, the simplified model caused some differences from the human gait pattern. For example, the lack of toe joints led to smaller peaks in the latter phase of the horizontal and vertical ground forces (Figure 2B) and less ankle dorsiflexion (Figure 2C) compared to the measured data. These trends were also observed in Song and Geyer (2015).

To set the UTTA condition, we attached a lower limb prosthesis to the right leg of the normal model. This manipulation caused the model to immediately fall. After 2000 generations of GAs, the symmetric control model and asymmetric control model reacquired stable locomotion, respectively (Figures 3A, 4A).

In the symmetric control model, smaller steps were observed compared to those of the other models (Figure 3A). The locomotion speed of the model was 0.32 m/s. The waveforms of GRFs, joint angles, and muscle activities were similar on the intact limb and prosthetic limb (Figures 3B–D). In addition, the peak activities of hip muscles (IL, GM, and RF) were larger than those of the other muscles, which was consistent with previous findings that hip muscle strength is related to walking ability in lower limb amputees (Nolan, 2012; Crozara et al., 2019).

In the asymmetric control model, the waveforms of GRFs, joint angles, and muscle activities differed between the intact limb and the prosthetic limb (Figures 4B–D, respectively). The locomotion speed of the model was 0.54 m/s. In addition, we observed that the peak activities of the hip muscles (IL, GM, and RF) were larger than those of the other muscles, similar to the symmetric control model. The phase space plots of these models were shown in Supplementary Material S3.

Figure 5 shows the metabolic cost of transport for the qualitative evaluation of locomotion. The normal model had the lowest value (3.23 ± 0.47 J/kg/m) compared to the other models, which was approximately consistent with the findings of previous studies (Russell Esposito and Miller, 2018; Das Gupta et al., 2019; Koelewijn et al., 2019). In the UTTA condition, the symmetric control model had a value of 10.10 ± 0.42 J/kg/m, meaning approximately two times the value in the asymmetric control model (4.97 ± 0.51 J/kg/m). Previous computer simulation studies of locomotion under UTTA conditions demonstrated that the metabolic costs of transport were comparable to or less than those measured in non-amputees (Handford and Srinivasan, 2016; Russell Esposito and Miller, 2018; Miller and Russell Esposito, 2021). These studies assumed that individuals with lower limb amputation can ideally acquire kinematic gait patterns similar to those of non-amputees due to the lack of an embedded neural model. Such findings diverged from clinical practice, in which various pathological gait patterns existed (Ichimura et al., 2022). In the current study, we included a neural model in the musculoskeletal model to generate various gait patterns. Thus, the different values of metabolic costs of transport in the present study may correspond to distinct gait patterns in individuals with lower limb amputation.

Figure 6 illustrates the spatiotemporal gait patterns. In the normal model, both stance time and step length were similar in the left and right limbs (Figure 6A). Under UTTA conditions, stance time was shorter in the prosthetic limb, and step length was shorter in the intact limb (Figures 6B, C). The step length in the intact limb indicated that the intact limb was the leading limb, requiring single support in the prosthetic limb. Thus, our results suggested that prosthetic limbs in the model had weak support, which was generally consistent with measured data (Sanderson and Martin, 1997; Nolan et al., 2003). In the ASI, the UTTA condition displayed an asymmetric gait pattern, showing remarkable differences compared to the normal model (Figure 6D). However, we observed little differences in ASI between the symmetric model and the asymmetric model (stance time ASI [%]: normal model, 0.03 ± 0.76; symmetric control model, 8.14 ± 0.64; asymmetric control model, 9.83 ± 0.76; step length ASI [%]: normal model, 1.53 ± 3.60; symmetric control model, 31.67 ± 4.87; asymmetric control model, 38.37 ± 6.22). The lack of substantial differences suggested that the symmetric and asymmetric control models were difficult to distinguish based on spatiotemporal gait factors. This was consistent with previous studies indicating that functional assessment of gait ability is difficult to achieve using gait asymmetry (Hof et al., 2007; Roerdink et al., 2012).

The simulation of gait reacquisition for individuals with UTTA in this study was based on physiological findings, which has the validity of replicating real-world situations. Darter et al. (2017) demonstrated that individuals with UTTA are capable of adaptations to locomotion comparable to non-amputees, and Rossignol et al. (2006) showed that the spinal cord network could adapt to sensory feedback signals. In addition, motor reorganization occurs in neural circuits after lower limb amputation (Chen et al., 1998), and the neural activity of specially trained athletes, such as Paralympians, exhibited a reorganization that differs from that of able-bodied individuals (Nakazawa, 2022). Such physical changes and training could modify the control strategy of the lower limb. Based on these findings, the present study re-searched only the sensory feedback parameters to the spinal cord model (i.e., the CPG) resulting in the reacquisition of a stable gait in the UTTA condition. The results of this study suggested that the essential factors for reacquiring locomotion under pathological conditions could be found using the computer simulation.

For activities of daily living, individuals with UTTA may require a new gait control strategy with low metabolic costs of transport, such as the asymmetric control model. How can this new behavioral control strategy be acquired? Such acquisitions of control strategies have been reported in hand rehabilitation after brain injury (Murata et al., 2008). In monkeys with brain injury, the pattern of grasping behavior changes with sufficient training, causing temporarily low grasping success rates, then markedly higher rates. In contrast, untrained monkeys could grasp, but their grasping success rates were low, and their grasping behaviors remained poor. These findings suggest that a new behavior control strategy could be acquired through sufficient training (Biernaskie et al., 2004; Sugiyama et al., 2013; Riazati et al., 2022). Similarly, in our study, the acquisition of the new gait control strategy improved the quality of locomotion in terms of low metabolic costs of transport. In other words, the symmetric control model, which is the same gait control strategy as that before the amputation, may be the initial state of training and may change to the asymmetric control model after training. Therefore, these two different gait control strategies may represent separate periods of the training process.

Previous studies (Nolan, 2012; Miller et al., 2017) reported that individuals with UTTA improved their balance and walking ability after training. Most studies, however, have focused on active individuals with UTTA, which means that the factors related to walking acquisition in the early stages of training remain unclear. In addition, fewer than 20% of lower-limb amputees can walk independently (Kamrad et al., 2020), requiring investigations of this factor for locomotion acquisition. To the best of our knowledge, this is the first study to investigate such factors for locomotion acquisition in individuals with UTTA. The results of the present study suggested that differences in gait control strategies modulate locomotion qualitatively, which were difficult to identify with easily observable spatiotemporal gait factors. As reported in previous studies (Hof et al., 2007; Roerdink et al., 2012), assessments for individuals with UTTA need to be tailored towards the initial, transitional, and other periods of rehabilitation, rather than single assessments of gait asymmetry. That is, even if gait asymmetry exists during the locomotor acquisition process, such gait may represent a necessary phase and may not require forceful correction. In contrast, gait asymmetry may increase the risk of osteoarthritis (Norvell et al., 2005; Amma et al., 2021) and may require an adjustment in activity level.

The musculoskeletal model constructed in this study is limited to two dimensions, and the CPG model is mathematically abstracted. For example, the model would not be appropriate for detailed motion analysis in 3D space as reported by Bruel et al. (2022) or for studies of neural activity during locomotion in the cerebrum as reported by Ausborn et al. (2019). Rather, the model represents a minimal closed-loop system of human locomotion, allowing for the generation of essential behavior resulting from some operation, such as a pathological situation. In addition, we employed the same evaluation function to search the parameters in the normal and pathological condition models. For example, landing pain during walking would play a more critical role in locomotion than metabolic costs. Therefore, optimization methods for severe pathological simulation may need to be considered based on the characteristics of the disorder.