Investigating the roles of reflexes and central pattern generators in the control and modulation of human locomotion using a physiologically plausible neuromechanical model

Objective. Studying the neural components regulating movement in human locomotion is obstructed by the inability to perform invasive experimental recording in the human neural circuits. Neuromechanical simulations can provide insights by modeling the locomotor circuits. Past neuromechanical models proposed control of locomotion either driven by central pattern generators (CPGs) with simple sensory commands or by a purely reflex-based network regulated by state-machine mechanisms, which activate and deactivate reflexes depending on the detected gait cycle phases. However, the physiological interpretation of these state machines remains unclear. Here, we present a physiologically plausible model to investigate spinal control and modulation of human locomotion. Approach. We propose a bio-inspired controller composed of two coupled CPGs that produce the rhythm and pattern, and a reflex-based network simulating low-level reflex pathways and Renshaw cells. This reflex network is based on leaky-integration neurons, and the whole system does not rely on changing reflex gains according to the gait cycle state. The musculoskeletal model is composed of a skeletal structure and nine muscles per leg generating movement in sagittal plane. Main results. Optimizing the open parameters for effort minimization and stability, human kinematics and muscle activation naturally emerged. Furthermore, when CPGs were not activated, periodic motion could not be achieved through optimization, suggesting the necessity of this component to generate rhythmic behavior without a state machine mechanism regulating reflex activation. The controller could reproduce ranges of speeds from 0.3 to 1.9 m s−1. The results showed that the net influence of feedback on motoneurons (MNs) during perturbed locomotion is predominantly inhibitory and that the CPGs provide the timing of MNs’ activation by exciting or inhibiting muscles in specific gait phases. Significance. The proposed bio-inspired controller could contribute to our understanding of locomotor circuits of the intact spinal cord and could be used to study neuromotor disorders.


Introduction
Limbs' movements result from the complex interaction between brain centers, the spinal cord, and the musculoskeletal system [45].The spinal network is essential in the control, coordination, and modulation of locomotion [32].While there is direct evidence of a central pattern generator (CPG) in mammals and other vertebrates [1,32], the lack of direct experimental access in humans means that there is only indirect evidence [39].Furthermore, sensory feedback pathways may play a major role compared to other mammals, and lower vertebrates in humans [24-26, 28, 34].Different studies suggested that the muscle activity observed during human locomotion may be controlled by five locomotor primitives that could be generated by rhythmic neural circuits [14,30].Investigating these questions about the roles of different spinal components in controlling locomotion is challenging since we possess only partial knowledge of the interactions between the different subsystems involved in this process.On top of that, limited experimental access complicates the observation of the sub-components functions leading to difficulties in model validation.Computer simulations are necessary and have been proven useful in the past to evaluate the contribution of each control component by evaluating different models, and parameters [3,18,19,22,29,41,48].
Various neuromechanical models have been proposed in the past to address these questions.In 1995, Taga proposed a musculoskeletal system controlled by a neural rhythm generator (RG) composed of seven pairs of neural oscillators and simple sensorymotor signals [49].Successively, in 2001, Ogihara and Yamazaki developed a neural controller composed of motoneurons (MNs) receiving inputs from a common CPG and reflexes from stretch and force receptors, where the spindle reflexes had inhibitory inputs to antagonist's muscles [40].In the context of locomotion controlled by CPG mechanisms, Aoi et al (2010) constructed a CPG model based on a two-layered hierarchical network composed of a RG and a pattern formation (PF) layer.The RG model produced rhythmic information using phase oscillators and was regulated by phase resetting based on foot-contact gait events, whereas the PF model generated feedforward commands composed of five motor primitives based on the muscle synergies analysis performed by Ivanenko et al [3,30].On the other hand, Geyer and Herr demonstrated that the kinematics and muscle activation observed in human locomotion could be reproduced without CPG commands by relying purely on sensory feedback activated at specific gait cycle phases [22].The activation of sensory feedback was regulated by a statemachine mechanism (i.e. a set of if-then-else rules) that enabled and disabled specific reflexes depending on the detection of stance or swing phases.A similar controller with partial modifications has then been proposed by Ong et al [41].In these studies, the activation of sensory responses in the gait cycle is regulated by a state-machine mechanism activating and deactivating reflexes in five gait cycle subphases.The necessity of including such state-machine mechanisms in reflex-based controllers hints at the need for a more sophisticated circuit that controls the underlying reflexes.Other studies have integrated CPG commands on top of these purely sensory-based controllers, showing the benefits of rhythmic circuits in gait modulation [19,53].
In this study, we propose a novel bio-inspired controller composed of a feedforward network inspired by Aoi et al consisting of two CPGs that produce the locomotor rhythm and patterns and a new physiologically plausible implementation of spinal reflexes based on neurophysiological studies in locomotion [54,56] without relying on any state machine mechanism.This network controls nine muscle actuators generating torques in a previously assessed musculoskeletal model [17].
The performance of this controller in replicating the behavior of human locomotion and its modulation are investigated and compared with previous experimental and neuromechanical studies.In addition, we investigate the performance of the sensory feedback controller alone to verify whether it is possible to generate human walking behavior with a purely reflex-based controller without relying on state machine mechanisms and to verify the benefit of CPG mechanisms.Finally, we examine the contribution given by pattern generation and reflex circuits to the MNs at slow, intermediate, and fast speeds performing a correlation analysis to identify possible parameters responsible for speed modulation.With these experiments, we aim to address the following questions: • What is the role of CPGs and reflex circuits in the generation of muscle activation in human locomotion?• Can low-level feedback circuits produce stable locomotion without a CPG or state-machine?• Is the contribution of these two neural components changing with increasing of gait speed?
Our results show that the reflex rules implemented in previous models [22,41] could be reproduced into less abstract and more realistic models of neural circuits.The insights given by the proposed controller suggest that spinal reflexes alone could not reproduce rhythmic locomotion without a state machine mechanism regulating the activation of reflexes in specific phases of the gait cycle.CPG networks appear to play the role of state machines in previous models and to be necessary to promote muscle activation in specific gait cycle phases.In addition, the modulation of CPG frequency seems necessary to modulate step duration.The modulation of either reflexes, CPG network, or both could generate gaits in a wide speed range, highlighting the high level of versatility of the neurospinal control of human locomotion.

Methods
This study used the Simulated Controller OptimizatioN Environment (SCONE) software simulation framework [21].SCONE is an open-source software performing predictive simulations of biological motion, allowing the investigation of individual models and control parameters on the motion.
To perform a predictive simulation in SCONE, it is necessary to run a scenario including the following components: • A mechanical model of the system to simulate.
• A controller that generates input for the model actuators.
• An objective that describes the target task to optimize, through a weighted combination of measures (sub-components of the cost function).• An optimizer that optimizes the free parameters in a scenario for a specific objective.
In this study, SCONE software was extended to implement and optimize the new spinal model generating gait simulations of 10 s.

Musculoskeletal model
The musculoskeletal model (figure 1) has the skeletal structure presented by Delp et al [17] with a height of 1.8 m and weight of 75.16 kg.The model is constrained in the sagittal plane and has a total of nine degrees of freedom (DoFs): a three-DoFs planar joint between the pelvis and the ground and three rotational DoFs per leg: hip flexion/extension, knee flexion/extension, and ankle dorsiflexion/plantarflexion.Three spheres per leg are used as contact models with the ground: one of radius 5 cm at the calcaneus and two of radius 2.5 cm at the toes.The musculoskeletal model is actuated by nine Hill-type muscle-tendon units per leg: gluteus maximus (GMAX), biarticular hamstrings (HAMS), biceps femoris short head (BFSH), rectus femoris (RF), iliopsoas (ILPSO), vasti muscle group (VAS), gastrocnemius (GAS), soleus (SOL), and tibialis anterior (TA).

Controller
Muscle activation is regulated by the excitation provided by the MNs.reflex controllers provide inputs to all muscles and are not regulated by any state-machine mechanism.We chose to maintain the state machine for the balance controller in order to simplify the balance control since our main goal is the simulation of locomotor movement.A physiological neuromechanical model of trunk balance control is a complex task that is outside the scope of this study.Muscle excitation is triggered by the MN output with values between 0 and 1 since MNs can only provide excitation to muscle fibers and cannot have negative outputs.To keep a reasonable level of abstraction and complexity, we will assume that the neuron's output (corresponding to its firing rate) u output follows the dynamics of a leaky integrator: where y is the neuronal response, u input is the neural input, τ the time constant (typically 0.01), u output the output of the neuron, and f the activation function.
The activation function used for MNs is the min-max operator (f(x) = min(max(0, x), 1)), and the neural input is defined as: where w j is the weight associated with the jth connection, and u output j the output of the jth neuron.The MN receives inputs from the CPGs' network (u CPGs ) and the reflex circuit (u reflexes ).These inputs are integrated according to equations (1) and (2): and generate the MN output m output .ILPSO, GMAX, and HAMS also receive inputs from the balance controller (u balance ).To avoid the activity of the balance controller from being inhibited by the other neural circuits possibly preventing the correct balance of the trunk, u balance is not integrated into the MN dynamics, and the final MN output for hip muscles moutput is defined by the following equation: where m output is the MN output resulting from the integration of u CPGs and u reflexes and u balance represents the balance controller effect on the hip muscles.The operator () + represents only the positive part of the selected signal.The amplitude of all components is regulated by the controller's parameters tuned by the optimization algorithm.The muscle activation a responds to the excitation m output (or moutput for hip muscles) as defined by Thelen [51].
The following sections will describe in detail how each neural input is computed (u balance , u CPGs , and u reflexes ).

Balance controller of the trunk
The balance controller is the one proposed by Ong et al [41], and it is the only controller part where a state-machine mechanism is present.A PD control strategy is used to activate the hip muscles balancing the forward lean angle of the trunk.ILPSO, GMAX, and HAMS receive inputs from the balance controller during the stance phase.The excitation given by the balance controller to the hip MNs is described in equation (5).(5) where k p and k v are the proportional and derivative controller's gains, and the constant θ 0 is the desired forward lean angle regulating the proportional feedback of the actual forward lean angle θ. t D represents the time delay, corresponding to t D = 5 ms for the hip muscles.The balance controller has a total of nine parameters.

CPGs
The CPGs were implemented as two coupled oscillators (one per side) composed of RG, and PF layers [36,46] inspired by the work of Aoi et al [2,3].Among the different CPGs models proposed for human locomotion, we chose to take inspiration from Aoi's model [3,4] to maintain a reasonable level of abstraction for a neural network for which the specific structure is not yet completely clear.Furthermore, the representation of five locomotor synergies is supported by past experimental studies [30].
The RG dictates a periodic command synchronized with the environment through afferents triggered by the heel-strike event.Based on Aoi's model, two coupled differential equations govern the CPG dynamics: where ϕ left/right denotes the phase of each leg, ω(t) is the basic angular frequency, and γ is the coupling constant.
The differential equation contains events that reset ϕ left/right when the leg touches the ground in order to synchronize the CPGs' phases with the environment.This is the only feedback mechanism present in the CPGs model, and it is described in equation ( 7): if the leg touches the ground In our simulations, the angular frequency ω has a constant value and represents one of the parameters under optimization.The pattern formation layer is composed of phase-dependent primitive patterns.Each pattern resembles a bell-shaped waveform with a defined width that can be centered around a specific phase value and is implemented as a raised-cosine function: where φ is the normalized gait phase, µ is the value corresponding to the peak of the bell shape, and σ is the half-width of the curve.The pattern formation layer is composed of five primitives of the same halfwidth and centered at different times of the gait phase (figure 3(a)): The choice of modeling the CPG network as the generation of five locomotor primitives has been taken to have smooth and derivable shapes of CPG circuits inputs to MN.More precisely, we took inspiration from the observations done in past experimental studies where five bell-shaped synergies active at different phases of the gait cycle were identified in human studies [30,31].However, it should be noted that those recorded synergies account for the overall MN activity and not just the CPG subcomponent.Each MN receives a weighted neural excitation or inhibition u CPGs from all primitive patterns (figure 3(b)) according to the following equation: where w m,k is the weight parameter of the pattern k to the MN m to be determined through optimization.The total number of parameters to optimize corresponds to five weights per pattern to each specific muscle and the oscillatory frequency.Therefore, the number of parameters for the CPG network is 48.The possible values assigned to w m,k are [−1:1].

Spinal reflexes
To implement a physiologically realistic model of sensory-motor control in human locomotion, we model and investigate five spinal reflexes: This reflex models the excitatory role of group II afferents [35] responding to changes in muscle length during stretch.MNs when a large activity is detected to prevent excessive muscle activation.
The spinal sensory feedback network is composed of three types of leaky integrator neurons: SNs, INs, and MNs.Each of these neurons models the activities of neural populations in the physiological spinal cord.
INs have the same properties of MNs responding to the dynamics described in equations ( 1) and ( 2) and with the same activation function f (x).SNs instead presents a rectifier function (f(x) = max(0, x)) as activation function, and also the neural input is slightly different: where r(t − ∆) is the receptor function and ∆ the delayed value of the receptor.Transmission delays are known and can be determined according to the proximity of the receptors [19].The expressions of the receptors follow the equations: where ṽm (t)), lm (t), fm (t), ff (t) are respectively the normalized quantities for contraction velocity, muscle length, muscle force, and cutaneous forces due to ground-foot contact.We choose to consider the normalized quantities to be able to easily scale for different muscles with different values of length and strength.Here, the expression for r Ia was inspired by Prochazka [43] and modified such that only lengthening ṽm > 0 triggers a response while ignoring length and activity-dependent terms.We deliberately simplified these expressions because we wanted to capture the general trend and prevent an excessive number of physiological parameters.In figure 4, we present the primitive reflex pathways that govern the connectivity within a single spinal cord segment.These rules are used to build the topological network by assuming that muscles can be categorized as agonists (A), antagonists (N), and extensors (E) or flexors (F).
The relation between agonist and antagonist muscles defines the mutual inhibitions described in figures 4(a)-(c) and (e).In addition, a muscle can be defined as extensor or flexor.In case it is an extensor muscle, the additional connections of Ib disynaptic extensor facilitation described in figure 4(d) are included.Some bi-articular muscles can be considered both extensors and flexors since they have different effects on different joints and the Ib disynaptic extensor facilitation is included also in this case.Table 1 describes the relations among agonist and antagonist muscles assigned in our models.Accounting for all the weighted connections, the sensory feedback controller has a total of 183 parameters.
Finally, figure 5 shows the whole spinal network implemented between agonist and antagonist muscles including the reflex pathways and CPGs inputs.On top of this network, ILPSO, GMAX, and HAMS also receive inputs from the balance controller.

Optimization process
In total, the controller's parameters are 256, accounting also for 16 additional parameters regulating initial positions and velocities of the model's DoFs.Because of the large size of the parameters' space and the difficulties in obtaining a stable solution when the network is in an arbitrary state, the optimization process is divided into three steps: imitation objective, optimization for stability, and optimization of metabolic energy.In the first stage, we try  determining the network's parameters such that the output of neurons is within a plausible range and MNs' activity resembles normal gait solutions.To achieve this, we begin with a previously obtained stable gait simulation generated by a simpler controller [41].Given that we know the whole state trajectories of the musculoskeletal system, we can compute the sensory afferent inputs required by the bio-inspired controller.Therefore, we can optimize for network parameters efficiently without numerically integrating the equations of the musculoskeletal system.We call this step imitation learning because we try to imitate a simulated behavior without yet producing dynamically consistent stable gaits.The optimization objective is defined as follows: minimize where e S m (t) denotes the target simulated excitation of muscle m at time t, and e N m the excitation of the network that depends on time t, the known state variables ⃗ x(t), and parameters ⃗ p.The above parameter solution does not produce stable gaits if we evaluate the model by numerically integrating the equations of motion.Our initial goal was to calibrate the network behavior within a reasonable range of operation in order to avoid neuron activities that are extreme and always make the model fall and optimization diverge.In fact, the imitation is only done to obtain a first usable solution for further optimization.
The second optimization aims at obtaining dynamically consistent stable gaits.To do so, we start integrating numerically the equation of motion producing dynamical gaits by minimizing the distance between the model's and reference's states as already expressed in equation ( 12), penalizing unstable falling solutions and solutions outside the desired range of speed, minimizing metabolic effort and joint limit torques.With this process, we aim to obtain stable solutions generated with our bio-inspired controller with physiological kinematics and muscle activation.The optimization is done using a CMA-ES algorithm with parameters λ = 40 and σ = 5 [27].The cost function for this optimization is defined as follows: The term J mimic , represents the model mimicking the reference states as expressed in equation (12), J gait penalizes the solution where the center of mass velocity is outside the [v min , v max ] range (1.10-1.25 m s −1 for healthy human gait at normal speed) and the falling solutions.The model is considered to fall when the ratio between its center of mass height (h COM ) to the initial state (h COM,i ) is smaller than a termination height threshold set to 0.8 ( hCOM h COM,i < 0.8).The term J effort defines the rate of metabolic energy expenditure [52] normalized by the product of body mass and distance traveled.J limit is associated with joint minimization of soft joint limit torques at the knee and ankle joints in order to avoid excessive joint angles [18].Finally, J head helps to maintain head stability by minimizing horizontal and vertical head accelerations outside the following ranges: [−4.90 − 4.90] m s −2 in the vertical direction, and [−2.45 − 2.45] m s −2 in the horizontal direction, as previously done by Ong et al [41].Concerning the weights, we assigned w mimic = 10, w gait = 100, w effort = 1, w limit = 0.1, and w head = 0.25 in order to promote mainly stability and mimicking.Following this optimization, we use the resulting stable solution as initial condition to find the optimal solution that minimizes metabolic energy.To do so, we remove the mimicking component of the cost function and optimize for minimize ⃗ p J gait + J effort + J limit + J head . ( J gait allows the stability of future explored solutions and J effort allows convergence toward gait efficiency.In addition, we apply external perturbations to the pelvis and randomized internal perturbations to muscle excitation to obtain more robust and stable gaits.The external perturbation is a force of 100 N applied in the forward and backward direction for a duration of 0.2 s respectively after 3 s and 4 s after the beginning of the simulation.The internal perturbations are applied to sensory receptors.For each controller timestep, a random white Gaussian noise is sampled from a normal distribution with a standard deviation of s * noise p , where noise p is the proportional standard deviation of the normal distribution, and s is the perturbed sensory signal.This three-step optimization process was used only to find a proper local optimum to replicate human gait behavior with a high number of parameters tuning the bio-inspired controller.However, once the local optimum is found, different gait behaviors can be reached starting from this solution by only optimizing according to equation ( 14) with the different gait behaviors targeted by J gait .The methodology for these simulations is explored in the following section.Furthermore, we verified that similar solutions for a specific gait speed could be achieved from different initial conditions of the optimizer as long as this initial condition represents a combination of parameters that has been obtained from the optimal solution resulting from the three optimization steps.

Gait modulation
To study the capability of the proposed bio-inspired controller to reproduce different gait behaviors in human locomotion, we focus mainly on the modulation of locomotor speed.In this way, we aim to evaluate our controller's performance, checking the maximum and minimum speeds it can achieve.Additionally, we evaluate gait analysis and muscle activation for three selected solutions far from the extremes of the achieved speed range since very slow or very fast speeds are more subject to producing artifacts in gait simulations.Therefore the three solutions selected are at 0.6 m s −1 , 1.2 m s −1 , and 1.6 m s −1 representing slow, intermediate, and fast speeds, respectively.To do so, we modulate the optimization parameters [v min , v max ] in J gait .Furthermore, we use the data acquired from our model to have possible insights into the contribution of CPGs and spinal reflexes in the neuromotor control of human locomotion and gait modulation.To do so, first, we evaluate the inputs to MNs from CPGs and reflexes and how these affect the MNs' output at different speeds.We then performed additional optimization where either CPG parameters or reflexes parameters were fixed to investigate the modulation capabilities of each controller component.The fixed values of parameters are extracted from a reference solution of the model walking at 1.17 m s −1 with 0.79 m of step length and 0.67 s of step duration.This solution is the one where the optimizer converged without imposing any restriction on the target speed.Finally, we investigate which parameters majorly contribute to gait modulation for the three controller configurations: full control, fixed reflexes, and fixed CPGs.These parameters are identified through a correlation analysis with gait speed, step length, and step duration, where parameters that have a high level of positive or negative correlation with these three gait characteristics (above 0.80 in absolute value) are considered the potential major contributors to gait modulation [18].The correlation analysis is conducted over eight solutions obtained through different target optimizations for each of the three controller configurations.
In general, the variables under investigation are: • Kinematic variables (measures in degrees).

• The correlation coefficients (unitless and bounded
between −1 and 1) between controller parameters and spatiotemporal variables.

Results
Figure 6 shows the controller's performance.The averaged gait cycle shown is obtained from 6 gait cycles.The first cycle has been removed to allow the model to reach steady-state.When minimizing the cost function without imposing restrictions on gait speed, the model converges to a gait at 1.17 m s −1 of speed, 0.79 m step length, and 0.67 s step duration.Figure 6(a) shows qualitatively the different positions of the model's joints through the gait cycle.The simulated pelvis tilt, hip flexion, knee angle, and ground reaction forces (GRFs) shown in figure 6(b) faithfully represent the experimental observations from Bovi et al [7] illustrated by the shaded gray areas.Some discrepancies can be observed for the ankle angle that tends to have excessive dorsiflexion and lacks proper plantarflexion during push-off compared to experimental observations.Indeed, the ankle angle mostly maintains its values above the zero level of plantarflexion/dorsiflexion.This likely depends on the weak activation of GAS and SOL observed in figure 6(c).
These muscles maintain a peak activation of 0.3 for the SOL and 0.2 for the GAS.However, the simulation replicates the temporal activations observed in experiments from Perry and Burnfield [42] for TA, GMAX, VAS, GAS, SOL, and HAMS.Concerning ILPSO, the muscle is active also outside its time range, having a consistent activation also in pre-swing.The model converges to different behaviors compared to experimental results for BFSH and RF that are active at the beginning and at the end of the gait cycle, respectively, rather than during swing.

Gait modulation
By optimizing the controller's parameters, the model could reproduce gaits from 0.3 to 1.9 m s −1 .Figure 7(a) shows the modulation of gait kinematics and GRFs at 0.6, 1.2, and 1.6 m s −1 representing slow, intermediate, and fast gaits, respectively.The averaged gait cycle in figure 7 is obtained through 5 gait cycles for the gait at 0.6 m s −1 , 7 gait cycles for the one at 1.2 m s −1 , and 8 gait cycles for the one at 1.6 m s −1 .Also in this case, the first cycle has been removed to allow the model to reach steady-state.As speed increases, the pelvis tilt and the lean angle of the trunk increase in the forward direction by 8 • , and the hip flexion oscillates between 35 • and −5 • at slow speed and 45 • and −3 • at fast speed.Increasing amplitudes of knee flexion are also observed at high speed, having the peak flexion in swing of 53 • at 0.6 m s −1 and 68 • at 1.6 m s −1 .Fast speed also presents a consistent increase of ankle plantarflexion to −3 • of the ankle angle during ankle push-off, whereas this value is maintained at around 10 • of ankle dorsiflexion at slow speed.Concerning GRFs, the characteristic double peak shape is very weak at 0.6 m s −1 .Double peak amplitudes increase with the increase of speed, especially the first peak that shows the reaction with the impact with the ground during heel strike.The duration of the stance phase is reduced from 65% of the gait cycle at 0.6 m s −1 to 55% at 1.6 m s −1 .The behaviors of kinematics and GRFs modulation presented resemble the ones observed experimentally by Bovi et al [7].Some differences are observed with the level of ankle dorsiflexion since the model tends to converge to a high-level of dorsiflexion that can differ from experimental data by 7 • during heel strike and by 12 • during push-off at slow and fast speeds.Additional differences are observed for the level of hip extension at slow speed and knee extension at high speed.In fact, Bovi et al observed that the maximum hip extension decreases at low speeds, and knee extension during stance increases at high speeds.In contrast, the model reproduced increased knee flexion during stance at high speeds and a similar level of maximum hip extension at 0.6 m s −1 compared to 1.2 and 1.6 m s −1 .Muscle activity is affected by gait modulation mainly through the increase of activation with the Figure 6.Gait analysis of simulated gait at 1.17 m s −1 .The average and 95% confidence intervals have been calculated over six gait cycles (a) Model position at different times of the gait cycle.(b) Kinematics and GRFs compared to experimental data [7]: gray areas report the observed experimental ranges for pelvis tilt, hip flexion, knee flexion, ankle dorsiflexion, and vertical GRFs.(c) Muscle activation analysis: muscle activity over a gait cycle for the nine muscles along the gait cycle.Blue curves represent the means of the gait signals through the gait cycles and the shaded areas the standard deviations.The activation curves are compared with the activation timing observed experimentally [42] and represented by the black horizontal bars on the top of the graphs.increase in speed.In figure 7(b), ILPSO, GMAX, HAMS, TA, and SOL are the muscles that more consistently present an increment in muscle activation.TA and HAMS pass from a maximum activation of 0.2 at slow speed to 0.5 at fast speed, whereas ILPSO has a similar maximum activation at fast speed and a higher activation of 0.3 at slow speed.GMAX has the highest increment of muscle activation, passing from a maximum activity of 0.1 to 0.7.SOL also presents a consistent increase in its activity, passing from a value smaller than 0.1 at slow speed to 0.4 at fast speed.A lower increase is present for VAS at high speed, whereas no consistent variation in muscle activity can be observed for BFSH, RF, and GAS.Therefore, the increased plantarflexion with speed mainly depends on the increased activity of the soleus.In general, the activation amplitude of all muscles increases with speed, as observed experimentally by Cappellini et al [10].

Gait modulation: CPGs and reflexes
We tested our controller to check whether the spinal connection implemented could generate rhythmic locomotion without the state-machine regulation or the presence of CPGs.With the removal of CPGs, the remaining parameters to optimize are 208.Even if the dimensional reduction could in principle simplify the convergence to a stable solution, no rhythmic gait could be reproduced, suggesting the need for the CPG networks to provide rhythm and timing in the absence of a state machine activating sensory feedback commands at specific times of the gait cycle (as in previous models).The simulation resulting from removing CPG parameters led to the human model in a standing position with the right leg in front of the left leg.The reflexes could generate the muscle activation necessary to maintain this position until the balance controller failed to stabilize the trunk, causing the model to fall.Dzeladini et al [19] suggested that the tuning of CPGs applied only to hip muscles could easily modulate human locomotion where other muscles were controlled by sensory feedback.In order to test this hypothesis in our model, we set to 0 the CPG parameters for all muscles except hip muscles and re-optimized the parameters according to the step explained in section 2.3.The resulting simulation showed very similar behavior to the one obtained without any CPG parameters, suggesting CPGs inputs may have an important contribution also for knee and ankle muscles.It should be noted that our CPG model provides only a rough waveform (made of the five primitives), while Dzeladini's CPG provided a detailed waveform replicating the sensory-driven control signals.
To investigate the contribution of each controller component in gait modulation, we investigated the inputs from CPG circuits, spinal reflexes, and balance controller provided to the MNs. Figure 8 shows how these signals contribute to generating MNs inputs and outputs following equation (4).Generally, in the model, the net effect of the reflex circuits tends mainly to inhibit the MNs providing a negative stimulation through the whole gait cycle with the exception of ILPSO and TA.Reflexes also facilitate the activation of VAS and GAS during swing for all the speed ranges and the activation of BFSH and HAMS at slow speeds.Instead, for each muscle, CPGs present specific regions of the gait cycle where they excite or inhibit the MNs.In this regard, CPGs prevent the activation of specific muscles in specific cycle phases, such as VAS and GAS in swing that were stimulated by the reflex circuits.CPGs' patterns tend to increase the amplitude of inhibition or excitation with increasing speed.This is especially the case for SOL, where the growing muscle activation with speed is primarily due to the increased excitation from CPG circuits.CPG activity also helps to have a consistent muscle activation of ILPSO in swing, but it tends to increase the activity at slow speed, and the lower muscle activity in swing is achieved by reflexes that inhibit ILPSO during swing at 0.6 m s −1 .The balance controller is applied only to ILPSO, GMAX, and HAMS, and seems to be the leading cause of HAMS activation since the CPGs excitation is entirely inhibited by spinal reflexes.
Additionally, to investigate the modulation capabilities of different neural circuits, we performed optimizations for different target speeds by keeping fixed reflexes parameters or CPGs parameters.Table 2 compares the achieved ranges of speed, step length, and step duration for the controller optimizing all parameters (full control), maintaining reflexes parameters fixed (fixed reflexes), and maintaining CPGs parameters fixed (fixed CPGs).The optimization of all parameters allows reaching wide ranges of speed from 0.30 to 1.86 m s −1 with small and large step lengths (0.23-1.08 m) and step durations (0.53-0.84 s).Removing reflexes parameters' optimization allows reaching ranges similar to the ones obtained in full control.However, the missing optimization of CPG parameters significantly limits the controller's capabilities to modulate step duration, passing from a range covering 0.53-0.84s to 0.64-0.67s.Consequently, the optimization tends to achieve slow or fast speeds, mainly modulating the step length to reach large values of 1.21 m at high speed.The achieved value of step length is higher than the one in full control because the model converges to a more energetically efficient gait reducing the step duration when all parameters are optimized.Furthermore, we verified that the modulation of CPG frequency alone is insufficient to converge to different gait behaviors.This result implies that the modulation of CPG frequency alone may be necessary but not sufficient to modulate step duration.

Correlation analysis
The correlation analysis reported in table 3 gives indications on which parameters had a correlation coefficient higher than 0.8 with the main gait characteristics and, therefore, those that could be the main responsible for gait modulation in the three controllers configurations.Specific parameters are identified as: • Pk → M.MN for the input pattern Pk weighted connections to MN M.MN, where M is the muscle name.
• M.N s A → M.N d A for parameters regulating the weighted synaptic connections between the source neuron of a specific muscle (M.N s A ) and the destination neuron of the target muscle (M.N d A ). N represents the type of neuron and can be either SN, IN, or MN, and A represents the type of afferent and can be either Ia, II, Ib, Ib+, or RC • M.N A .w0 is the activation offset of the neuron N A regulating the neuronal response.
In full control, both reflexes and patterns' connections seem to contribute to gait modulation.The first (P0) and third (P2) patterns connections to extensor muscles like GMAX and SOL positively correlate with increasing speed and step length and decreasing step duration.The CPGs' frequency (ω) has a highly consistent correlation with gait speed and step duration, suggesting again the direct influence of this parameter on gait frequency.The only reflex parameter representing an excitatory connection is the monosynaptic excitation of Ia afferents from TA (TA.SN Ia → TA.MN), having a negative correlation with speed and favoring increased dorsiflexion during slow gaits.with increasing speed, favoring the inhibition of the hamstrings muscle.Indeed, from the previous analysis, the increased activation of HAMS with speed was mainly due to the input from the balance controller.
Concerning the configuration with fixed reflexes, CPGs' frequency (ω) highly correlates with speed and step duration also in this case.Another modulator for step duration is the input from the fifth pattern to TA's MN (P4 → TA.MN), which indeed increases its activation at the end of the gait cycle with increasing speed.The first pattern (P0) tends to increase the inhibition to GAS and GMAX at the very beginning of the gait cycle with increasing speed.Then, speed modulation through modulation of step length is enhanced by tuning the excitation from the second pattern (P1) to SOL MN (SOL.MN) in order to increase propulsion in stance.
When CPG parameters are fixed, speed modulation happens mainly through changes of step length because of the controller's limited capability to modulate step duration without tuning CPGs' frequency.The controller tends to increase step length by increasing the offset to enhance the length feedback of ILPSO (ILPSO.SN II .w0).II afferents are also involved with the decreased length feedback of VAS muscle with speed rising through the excitation of VAS.IN II from VAS.SN II .The last two relevant parameters concern the inhibitory connections of RCs and Ib afferents.The RC IN of GAS (GAS.IN RC ) decreases its inhibition to the Ia IN of TA (TA.IN Ia ), decreasing the activation of GAS itself at fastest speeds.Higher speeds should, in principle, increase the activation of GAS, but in the modulation of muscle activation, we previously observed that the optimizer tends to maintain the same activation level for the gastrocnemius muscle during speed modulation.Finally, the Ib inhibitory IN of TA (TA.IN Ib ) decreases its inhibition to the Ib IN of SOL, allowing the inhibition of this muscle.Indeed, we previously observed that the increased muscle activation of soleus at higher speeds was not due to the input from spinal reflexes but primarily due to increasing excitatory inputs from CPGs.

Discussion
In this study, we aim to investigate the possibility of controlling human locomotion by relying only on spinal reflexes not regulated by a state machine mechanism and to investigate the contribution of both CPGs and spinal reflexes in generating locomotor output.To do so, we developed a bio-inspired controller composed of a balance controller, a CPG network, and a sensory feedback network based on physiological spinal reflexes maintaining a statemachine mechanism only for the balance of the trunk.The proposed controller regulated the stimulation of nine muscles per leg.The number of muscles was chosen to find a reasonable compromise between the controller performance and a limited number of parameters.However, the nine muscles do not represent the minimal number of muscles to obtain a physiologically plausible model.Indeed, a preliminary version of the controller has been tested with a musculoskeletal model actuated with seven muscles per leg reproducing faithful but less optimal human kinematics and muscle activation.The proposed controller could replicate human kinematics and GRFs with some limitations in the ankle angle in which the model converges to an excessive dorsiflexion behavior.Regarding muscle activation, the model could reproduce most of muscle activation timings observed experimentally, with the exception of BFSH, which is active outside its range in human recording.The proposed network could probably generate muscle activation closer to physiological activity with additional optimizations.However, finding this global optimal solution remains challenging because of the large number of parameters.
Many aspects of speed modulation from human recordings, such as the increased amplitudes of flexion/extension movements and the increased muscle activation with growing speed, are also matched by the model.Concerning the role of CPGs and spinal reflexes in the neural control of human movement, we investigated the possibility of finding stable solutions without relying on CPGs as suggested by previous neuromechanical studies [22,41].In our optimizations, we could not find any stable rhythmic behavior in the absence of CPGs' commands even if the number of parameters to optimize significantly decreased, suggesting the need for the CPG network to provide rhythm and timing in the absence of a state machine activating sensory feedback commands at specific times of the gait cycle.Therefore, reflex-based circuits are always active through self-regulation by the afferents, lacking any timing information without CPGs.Similarly, a pure CPG network without reflexes leads to unstable solutions.While we cannot rule out that a different network topology might give rise to highquality gaits, our model highlights the need for both types of neural mechanisms to achieve stable and natural movements.Indeed, natural locomotor behavior emerges when both CPGs and spinal reflexes are active.Our study suggests that the state-machines used in previous sensory-driven models [23,41] could, in fact, be replaced by CPGs and that one of the main roles of CPGs, in addition to simplifying speed control [19], is to serve as a gating mechanism that ensures that reflexes do not affect muscles all the time but only at specific moments of the locomotor cycle.
More specifically, the performances on gait modulation while either reflex circuits or CPG commands were fixed, and the corresponding correlation analysis highlighted the importance of CPGs' frequency in changing the step duration.Therefore, in the model, CPGs have a crucial role in determining gait timing.Additionally, the analysis of neural inputs to MNs showed that the net inputs of reflexes are mainly inhibitory through the gait cycle for the proposed model, except for ILPSO and TA, which globally receive excitatory inputs.CPG patterns excite or inhibit MNs in specific phases of the gait cycle to allow or prevent muscle activation.Therefore, CPGs seem to be important to determine activation timing other than gait frequency.Such a control strategy is similar to the one proposed by Laquaniti et al [33] where the timing and magnitude of electromyography (EMG) activity are tuned via proprioceptive feedback and CPGs that control the basic rhythms and patterns of MN activation.However, it should be highlighted that the five locomotor primitives described by Ivanenko et al [30] and Laquaniti et al [33] were not equally spaced in the gait cycle phase as they are in our controller.This is because, in these studies, the primitives were extracted with factorization of EMG activity.Yet, this activity is the result of the global input received by muscles without being able to distinguish which input was coming from spinal reflexes and which one from CPG circuits.Therefore, we decided to simplify the distribution of the five primitives and equally space the patterns through the gait cycle since the primitives measured in experiments could hardly be generated by the CPG commands alone.This choice still leads to largely reproducing the experimental activation timing.
The modulation of gait reflexes alone could still regulate muscle activation to achieve different gait behaviors, mainly through the modulation of step length.The correlation analysis highlighted the possible parameters responsible for this behavior, such as the offset of II fibers regulating the level of stretch necessary to activate length feedback for ILPSO and GMAX.Indeed, increasing these parameters allowed larger amplitude for hip flexion/extension, promoting larger step lengths.
In general, the proposed controller presents a highly redundant system where several different combinations of neural inputs can generate the same muscle activation.The correlation analysis gave possible insights into which parameters could be the most relevant in the control of gait modulation.Yet, given the high redundancy, a separate and more extensive study would be necessary.Possibly, this study should include a large dataset of optimizations and additional elements of the cost function that could guide toward the best combination of neural inputs to generate specific muscle excitation, such as the minimization of the total neuronal activation.Then, the results should ideally be validated by experimental measurements.
Some limitations of the proposed controller should be considered.Because of the large number of parameters, finding a stable solution replicating human walking with the proposed controller may be challenging since it requires the three optimization stages described in section 2.3.However, once this solution is found, it can be used as a starting point to explore different gait behaviors by only performing the last optimization stage.In this way, we could reproduce a wide range of speeds comparable to or larger than the ones previously obtained by other neuromechanical controllers [18,41,47].Yet, the used cost function is very sensitive to the changing of optimization variables.Therefore, the optimization method used might require a good initial guess.For this reason, the initial stage of the optimization requires the imitation objective from a previous solution found with a different neuromechanical controller.Without these steps, there is a low probability that a random initialization of parameters can make the convergence to a stable solution.However, the use of the imitation objective implies that any lack of performance from the imitated solution in replicating human movement will probably reflect a lack of performance of the bio-inspired controller.This could be the cause of the excessive dorsiflexion behavior performed by our model since many solutions of the reflex-based controller proposed by Ong et al [41] that we used as imitation objective presented the excessive dorsiflexion behavior.Therefore, the proper choice of the initial imitation objective is crucial for the correct optimization of our model.
Further considerations should also be made for the design of the reflex controller.In paragraph 2.2.3, we explained how we simplified the expressions for the sensory receptors to capture the general trend and prevent an excessive number of physiological parameters.In reality, the dynamics of these receptors are very complex [37,38], and there is little evidence why the same model identified in specific animal experiments can generalize to humans in the presence of dynamic movements.
Another aspect currently not captured by our model is the spatial organization of muscles neural control in the spinal cord [31].Indeed, past studies suggested that MN pools innervating muscles performing synergistic action are co-localized in the spinal cord, and nearby MN pools are likely to share the proprioceptive feedback from the muscles they innervate [50].The model we propose in this study could potentially be extended by coding specific spinal segments including the facilitatory and inhibitory mechanisms among synergistic and antagonist muscles presented in this study and having their MN pools located in that specific segment.Despite these limitations, the bio-inspired controller we propose is a promising tool for investigating spinal circuits in human locomotion.Indeed, we have already shown the insights this model could give into the relationship between CPGs and spinal reflexes.Further suggestions could be provided in investigating pathological gaits.Past studies tried reproducing neural pathologies with neuromechanical simulation by extending previous controllers, including specific connections to model the pathology in the desired degree of freedom [8].However, the controller proposed could be more suitable for studying neuropathologies like hyperreflexia considering the effects of both excessive inputs from Ia fibers and the lack of reciprocal inhibition.Furthermore, further aspects of gait modulation regulating standingto-walking transitions and acceleration and deceleration mechanisms can be investigated.Additionally, this controller could be used as a starting point to further extend the modeling of the neuromotor system by including the implementation of additional spinal neural connections like γ-MNs [20] and descending inputs from the brainstem and other supra spinal brain areas, even though this would increase even more the controller's complexity and the total number of parameters.Additionally, future implementations could include less abstract and more realistic CPG models, for instance, based on more detailed models previously proposed for mammalian circuits that could potentially be taken as a reference for modeling human locomotion [5,6,12,13,15,16].Additional connections between CPGs and spinal reflexes may be implemented, allowing SNss to interact and modulate CPGs' patterns and CPGs' patterns to interact with spinal INs other than MNs.

Conclusions
This study proposes a novel physiologically plausible neuromechanical controller maintaining a good balance between complexity and realism to investigate the spinal components governing human locomotion.The controller is composed of a balance controller from Ong et al [41], a CPG network inspired by Aoi et al [3], and a sensory feedback network that takes into account the main reflex connections in the spinal cord without being tuned by a state machine.The controller demonstrated the ability to reproduce key behaviors of human locomotion and its modulation in simulations.Results from optimizations suggested that rhythmic locomotion could not be achieved with the only contribution of spinal reflexes without accounting for a state machine mechanism.This suggests the possible need for CPG networks to generate rhythmic movements by guiding muscle activation timing in specific phases of the gait cycle.The modulation of either CPGs or reflexes parameters or both could reproduce wide ranges of gait behaviors, highlighting the high level of redundancy in human locomotor control.The modulation of CPGs' frequency appeared to be crucial for regulating gait cycle duration.The proposed controller demonstrated to be a promising tool to provide many other indications on how the spinal cord may produce locomotor outputs.Additionally, to the authors' knowledge, no previous neuromechanical controller for human locomotion modeled the human proprioceptive feedback at a low level of abstraction like in this model.Extensions of the proposed controller could potentially test if more detailed spinal networks observed in mammalian circuits [5,6,12,13,15,16] could reproduce human biped movements in closed loop.

Figure 2 .
Figure 2. Control diagram: the bio-inspired controller is composed of the balance controller of the trunk and the spinal controller divided into CPGs and spinal reflexes.The balance controller aims at keeping the balance of the trunk by stimulating the hip muscles' motoneurons, whereas CPGs and spinal reflexes generate rhythmic behavior stimulating all muscles' motoneurons.Reflexes and CPGs are integrated by motoneurons, whereas balance inputs are summed separately.Created with BioRender.com.
Ia afferents provide monosynaptic excitation to MNs innervating the same muscle and disynaptic inhibition mediated by Ia inhibitory INs to antagonistic MNs (figure 4(a)) and model the velocitydependent response to stretch [11].• II afferents provide disynaptic excitation to MNs innervating the same muscle and disynaptic inhibition to antagonistic MNs mediated by excitatory and inhibitory INs, respectively (figure 4(b)).

Figure 3 .
Figure 3. CPG structure: (a) the CPG generates five bell-shaped primitives centered at different times of the gait cycle.(b) Each k-pattern stimulates all the m-motoneurons depending on the assigned weight w m,k that can be positive or negative.(a) and (b) Created with BioRender.com.

Figure 4 .
Figure 4. Reflex pathways.Green, blue and red stand for somatosensory (SN A/N ), inter (IN A/N ), and motor neurons (MN A/N ), respectively.The connection tip o stands for inhibition while < is for excitation.Subscript letters A, N, and E denote agonist, antagonist, and extensor muscles, respectively.The rules are repeated for all antagonist muscles.(a)-(e) Created with BioRender.com.

Figure 5 .
Figure 5. Spinal network between a muscle and its antagonist.The network includes reflexes driven by Ia, II, Ib afferents and Renshaw cells and inputs from CPGs' patterns.Connections from patterns to motoneurons are represented by back arrows since these connections can be both inhibitory and excitatory.ILPSO, GMAX, and HAMS also receive inputs from the balance controller.Created with BioRender.com.

Figure 7 .
Figure 7. Gait modulation at 0.6 (calculated over 5 gait cycles), 1.2 (calculated over 7 gait cycles), and 1.6 m s −1 (calculated over 8 gait cycles).(a) Comparison of kinematics and GRFs among the three speeds.(b) Changing muscle activation at different speeds.The activation curves are compared with the activation timing observed experimentally and represented by the black horizontal bars on the top of the graphs.

Figure 8 .
Figure 8. CPGs, reflexes, and balance inputs at 0.6 (calculated over 5 gait cycles), 1.2 (calculated over 7 gait cycles), and 1.6 m s −1 (calculated over 8 gait cycles).CPGs, reflexes, and balance generate motoneurons inputs and outputs according to equation (4).The contribution of the three controller components is compared through the three speeds selected.

Table 1 .
Agonist-antagonist relationship among muscles modeled.Each antagonist relationship implies the corresponding reciprocal inhibition of Ia, II, and Ib connections, and the reciprocal excitation connections of RC.The table specifies whether the agonist is considered an extensor, which includes the disynaptic excitation from Ib+, or flexor.

Table 3 .
Correlations coefficients of controller's parameters contributing to the modulation of speed, step length, and step duration in the three controller's configurations: full control, fixed reflexes, and fixed CPGs.GAS.IN Ia → TA.IN Ia ), enhancing the activation of GAS itself because of the decreased inhibition from TA.IN Ia .Then, the RC of RF decreases its inhibition to the RC of HAMS (RF.IN RC → HAMS.IN RC )