1 Introduction

Animals as diverse as tardigrades (Nirody et al., 2021) and tapirs (Catavitello et al., 2018) use interlimb coordination to walk through complex terrain. When a walking animal encounters an unexpected perturbation (e.g., it is pushed or tripped), its ability to recover and sustain locomotion can be a matter of life and death. However, the presence of significant temporal delays in animal sensorimotor systems establish fundamental limits on how quickly animals can respond to external perturbations (More and Donelan, 2018). Sources of temporal delay include neural conduction (Sterling and Laughlin, 2015), synaptic transmission (Sabatini and Regehr, 1999), and electromechanical muscle activation (Pringle, 1978; Akiyama et al., 2010).

Sensorimotor delays pose a particular challenge for systems that rely on feedback control, because they limit the availability of up-to-date information on the state of the external environment (Franklin and Wolpert, 2011). For example, time delays limit robust performance in feedback control systems (Doyle et al., 2013), with model predictive control (MPC) (Garcia et al., 1989; Camacho and Alba, 2013; Brunton and Kutz, 2022) providing one algorithmic strategy to improve performance while maintaining robustness. In biological systems, theoretical studies suggest that some animals use predictive internal models to mitigate the effect of delays (More and Donelan, 2018; Li et al., 2023; Desmurget and Grafton, 2000). Limb compliance and other biomechanical adaptations can also compensate for unexpected perturbations within a limited range and mitigate the effect of sensorimotor delays (Daley, 2018; Ashtiani et al., 2021). However, the constraints that sensorimotor delays impose on robust locomotion have been difficult to quantify, because it is experimentally challenging to manipulate sensorimotor delays and observe their effects on animal locomotion (although this has been achieved in other motor systems, for example using focal cooling in singing birds (Banerjee et al., 2021)). Computational models of locomotion do not typically include delay as a tunable parameter, and most existing models of walking cannot sustain locomotion in the presence of delays and external perturbations.

Flies are agile, robust walkers, and the availability of genetic and behavioral tools in Drosophila makes them well-suited for investigating neural mechanisms of locomotor control (Deangelis et al., 2019; Goncalves et al., 2022; Cruz and Chiappe, 2023). The fly is also the only walking animal whose nervous system is almost completely mapped at synaptic resolution (Galili et al., 2022). Flies walk rhythmically with a continuum of stepping patterns that range from tetrapod (where two of six legs lift off the ground at a time) to tripod (three of six legs alternately lift off at a time) (Deangelis et al., 2019; Mendes et al., 2013; Szczecinski et al., 2018; Wosnitza et al., 2013; Strauss and Heisenberg, 1990; Nishii, 2000; Pratt et al., 2024). Each fly leg has five joints that move through 7 mechanical degrees of freedom (Karashchuk et al., 2021; Lobato-Rios et al., 2022) and are actuated by ∼18 muscles that are innervated by ∼70 motor neurons (Azevedo et al., 2022). Each leg motor neuron is uniquely identifiable (Azevedo et al., 2022) and receives thousands of synaptic inputs from hundreds of unique premotor neurons within the fly’s ventral nerve cord (VNC, Lesser et al. (2023); Cheong et al. (2023)), a part of the invertebrate nervous system analogous to the vertebrate spinal cord. The architectural features of Drosophila locomotor control — joints, leg motor neurons, VNC — motivate the architecture of our proposed model.

Past models of insect multi-legged walking have taken three general approaches. The first approach models legs as coupled oscillators or inverted pendula without taking into account the mechanics of leg joints. Each leg is a single oscillator, and a network of oscillators is tuned to recreate oscillatory gait from measurements of foot falls (Couzin-Fuchs et al., 2015; Proctor and Holmes, 2018); alternatively, tripod gaits are approximated by spring-loaded inverted pendula (Chun et al., 2021). The second approach focuses on physical details of the legs and its joints. Schilling et al. (2013) uses a decentralized reactive controller to recreate hexapod walking patterns; Goldsmith et al. (2019) and Goldsmith et al. (2020) introduce a robotic platform; Lobato-Rios et al. (2022), Wang-Chen et al. (2023), and Vaxenburg et al. (2024) develop virtual fly simulations in a physics engine. Walking behavior in these models are also driven by coupled oscillators with tuned parameters that reproduce inter-leg coordination patterns; however, the underlying joint kinematics are typically unrealistic. The third approach uses normative learning and optimization to generate walking behaviors de novo — Ramdya et al. (2017) maximized fly walking speed by varying inter-leg coordination with a genetic algorithm, whereas (Heess et al., 2017) maximize walking speed in both bipedal and quadruped walkers with reinforcement learning. Normative models produce walking with varying degrees of realism, but they require clever selection of objectives with constraints and are computationally expensive. Further, any change in delay values would necessitate a computationally intensive retraining process. Notably, Geijtenbeek et al. (2013) include fixed sensorimotor delays in a bipedal model, but do not explore the effect of varying delay values. Overall, existing walking models offer kinematic or physiological accuracy, but few achieve both, and none consider the effect of varying sensorimotor delays.

Here, we develop a new model of fly walking with several key innovations. First, the model is trained on high-resolution 3D joint angle data from walking flies, which has only recently become possible due to new deep-learning based computer vision tools (Karashchuk et al., 2021). Second, the model’s multi-layered architecture achieves more than the sum of its parts: specifically, a neural network model recapitulates kinematic coordination of many joints, and an optimal controller allows the data-driven model to generalize to new scenarios (e.g. large delays and perturbations) without retraining. Third, the inclusion of delay as a tunable parameter allows us to systematically investigate the quantitative relationships between sensorimotor delays and robust walking. Fourth, we introduce a new method to quantitatively compare the kinematic similarity of real and simulated walking. Overall, these analyses suggest that neuromuscular delays limit how fast flies can walk while retaining robustness to unexpected perturbations. They also illustrate a general approach of using in silico experiments in virtual animals to investigate how fundamental physiological parameters constrain animal behavior.

2 Results

2.1 A kinematic model of fly walking that incorporates delays and accommodates perturbations

We designed a walking model with three functional layers (Fig. 1B), inspired by the hierarchical anatomical organization of the fly nervous system (Fig. 1A, Dallmann et al. (2021)). The three layers are an inter-leg phase coordinator, a trajectory generator, and an optimal controller that interfaces with a leg dynamics model. Each individual leg is modular and governed by its own dynamics, optimal controller, and trajectory generator. Inter-leg coordination is accomplished by the phase coordinator alone. In other words, the movement of each leg is not coupled to any other leg except through the phase of its current step cycle. This modularity is inspired by the segmental neuroanatomy of the VNC, in which each leg is controlled by distinct local premotor circuits and pools of motor neurons (Lesser et al., 2023; Cheong et al., 2023).

Figure 1.

Each layer in the model is an abstraction for the layer below it, such that various elements of walking (e.g. joint control, inter-leg coordination) can be integrated through modules. Below, we describe each modular component of our model; more details on its derivation and implementation are elaborated in the Methods and Materials and Appendices. Although the model can turn and side-step, our analysis focuses on forward-walking — the model is driving the fly to walk straight in all simulations, unless otherwise stated.

2.1.1 Inter-leg phase coordinator

To coordinate multi-legged walking, we modeled the step phase of each leg as an oscillator. We refer to the left and right legs as L1–3 and R1–3, respectively, where the front legs are L1 and R1. The phase coupling of the six leg oscillators establishes realistic inter-leg coordination. We use a Kuramoto oscillator model (Acebró n et al., 2005; Strogatz, 2000) to perform this coordination, as in Proctor and Holmes (2018); equations and implementation details are in Methods and Materials. Briefly, the phase coordinator takes the instantaneous phases ϕ from all legs as input, then outputs the desired phases ϕ_d to all legs. These desired phases are synchronized across pairs of legs to maintain a tripod coordination pattern, even when subject to unpredictable perturbations. The phase coupling coefficient among legs is estimated from measured 3D joint kinematics of walking flies (Karashchuk et al., 2021). Without the phase coordinator, individual trajectory generators would generate realistic kinematics for each leg, but they would not be coordinated with each other.

2.1.2 Joint kinematics trajectory generator

The trajectory generator layer is responsible for producing realistic 3D joint kinematics of each leg. To express the relationships among all leg joint angles, we use an artificial neural network trained to generate angle trajectories. The modular design makes it possible to train trajectory generators separately from the inter-leg phase coordinator and optimal controller. Ultimately, the optimal controller integrates desired angle trajectories with proprioceptive joint angle feedback to output joint torques to the physical model of each leg.

As illustrated schematically in Fig. 1B, the inputs to the trajectory generator module are desired leg phase ϕ_d, joint angles θ, joint angular velocities , and walking speed and direction of the fly. The trajectory generator then outputs the leg phase velocity and desired angular accelerations . This output is integrated to produce the desired angle and angular velocities , which are the inputs to the optimal controller.

To train the network from data, we used joint kinematics of flies walking on a spherical treadmill, obtained from tracking 3D joint angles with Anipose (Karashchuk et al., 2021). Details on training approach and network properties are described in Methods and Materials. The walking speed and direction are signals that are not generated from other modules of the model, but are instead external inputs to the trajectory generator computed from the data; biologically, these signals are analogous to descending signals from the fly brain. This organization is motivated by the observation that walking velocity and direction have a substantial effect on joint angles, but they do not have a substantial effect on parameters and outputs of the phase coordinator (e.g., phase offsets), as substantiated in Appendix B.

After training and when assembled with the other layers of the model, the trajectory generator receives proprioceptive information on the observed, current state of the leg , as well as the target desired phase ϕ_d from the phase coordinator. Thus, it generates a time-series of desired angles and angular velocities for some future interval and sends this time-series to the controller; it also estimates the current phase of the leg, which is passed to the phase coordinator.

In the absence of external perturbations, the trajectory generator produces realistic joint angles similar to those of walking flies, as we show below. When challenged with unpredictable external perturbations, the impact on the trajectory generator is mitigated by the optimal controller layer, which attempts to return the actual state to the desired state. This control is possible because the controller operates at a higher temporal frequency than the trajectory generator, so that the controller has adequate time to reject disturbances before communicating to the trajectory generator. We emphasize that all data used to train the trajectory generator came from experimental conditions with no external perturbations.

2.1.3 Control and dynamics

The optimal controller layer maintains walking kinematics in the presence of sensory and motor delays and helps compensate for external perturbations. At regular intervals, the controller receives a time-series of desired state trajectories θ_d from the trajectory generator layer. The controller then produces the necessary torques τ to track this trajectory for the given leg dynamics. External perturbations w, when present, enter through the dynamics and affect the state ; the controller then senses the state and responds accordingly.

To design the controller, we first derived dynamical equations for each leg using link-and-joint models (Fig. 1C), then linearized these dynamics and designed a linear quadratic regulator (LQR) controller. This optimal controller senses the state of the leg via proprioceptive input, then determines the optimal motor output for walking. The controller makes use of internal predictive states to accommodate sensory and motor delays. A detailed description of the controller derivations can be found in Methods and Materials. We note that the walking and compensation capabilities of the overall model are not contingent upon any specific dynamics or controller formulation; any controller that adequately tracks the trajectory generator would suffice.

2.2 The model generates realistic walking kinematics

The layered model generates realistic 3D walking kinematics when compared to real kinematic data from walking flies (Fig. 2). Below, we provide qualitative and quantitative comparisons of joint angles, joint angular velocities, and phases of walking both within and across legs.

Figure 2.

2.2.1 Qualitative evaluation of joint angle time-series and videos of walking kinematics

We first qualitatively compared simulated and real kinematics by examining time-series data and videos. In example trajectories of femur-tibia flexion angles of the right front and femur rotation angles of the left middle legs, the simulated time-series matched the mean, frequency, and pseudo-triangular shape of the fly data (Fig. 2A). Articulated animations of simulated and real trajectories for comparison are shown in Videos 1, 2, and 3. Although the model and data were largely similar, some differences stood out. For example, the simulated amplitudes were generally smaller than real amplitudes, and the simulated trajectories tended to be more regular.

2.2.2 Comparing joint angles and angular velocities versus phase

We next sought to quantitatively compare simulated and real joint angles. Direct comparisons of time-series trajectories are inadequate, because temporal offsets between time-series produce large mismatches even if the time-series are similar. For instance, if we shift a time-series trajectory by a half-cycle and compare it with itself, this will produce a large mismatch, even though the two trajectories are identical except for a misalignment in phase.

To compare real and simulated trajectories, we computed the step-cycle phase for each time-series and used this to plot the mean angles θ as a function of phase (Fig. 2B). Here we make the distinction between generated phase, the per-leg phases produced by the phase coordinator of the model, and computed phase, which can be computed for each joint from time-series data. Since we did not have access to the generated, desired phase for real fly data, all comparisons were made between computed phases.

When we averaged joint angles for all legs over 500 distinct walking bouts of 0.5 to 2 seconds in duration, we found that the mean differences between real and simulated joint angles were less than 6 degrees (Fig. 2C). This difference is comparable to the uncertainty associated with markerless tracking of 3D fly walking kinematics that we used as training data (5.56 degrees, from Karashchuk et al. (2021)). Errors for angular velocity were higher. Aggregate differences as a function of walking and turning velocity are shown in Appendix E.

The similarity between real and simulated data as a function of phase was consistent across the natural range of forward walking speeds (Fig. 2B). Plots for all legs and joints are shown in Appendices C and D. We further demonstrate the model’s capacity to simulate leg kinematics by comparing the phase coupling of simulated joint kinematics with real walking flies, both within and across legs, in Appendix F.

2.3 Model maintains walking under unpredictable external perturbations

When walking in natural environments, animals frequently navigate uneven or slippery terrain. Thus, robust sensorimotor control systems must detect and respond to such unexpected perturbations in order to maintain stable locomotion.

Here, we show that our model maintains realistic walking in the presence of external dynamic perturbations, despite being trained only on data of walking without perturbations. This performance is made possible by the combination of the trajectory generator and optimal controller in the model. By taking advantage of proprioceptive feedback, the optimal controller compensates for external perturbations, allowing the trajectory generator to sustain realistic joint angle trajectories.

We considered two types of perturbations: impulse and persistent stochastic. Impulse perturbations are analogous to when an animal experiences a brief, unexpected force (e.g., legs slipping on an unstable surface). We simulated impulse perturbations as a velocity that is added to all joints at a single time (i.e., a single time step in our simulation). The magnitude of this velocity is drawn from a normal distribution for each joint, where increasing perturbation strength increases the mean and variance of this distribution. Persistent stochastic perturbations displace all the legs, but not at the same time. We applied perturbations following a Poisson process with a mean rate of 10 Hz; in other words, over a period of time (1 second), for each leg at each timestep, we randomly select whether a impulse perturbation should be applied based on a Poisson distribution. Ranges of perturbation strengths used in our simulation correspond to estimates of biologically plausible values (as derived in Appendix G).

To evaluate the model’s ability to walk in the presence of perturbations, we compared time-series data before, during, and after perturbations across a range of walking speeds (Fig. 3B-C). For both types of perturbations, joint angle trajectories were different during perturbations, but recovered after the perturbation period ended. In the case of persistent stochastic perturbations, kinematics appeared approximately oscillatory even during perturbations, indicating that a semblance of walking was maintained. Example animations of simulated walking bouts with perturbations are shown in Videos 4–7.

Figure 3.

To quantify the extent to which perturbed kinematics resemble normal, unperturbed walking, we introduce a new quantitative metric termed kinematic similarity (KS). For a given window of a kinematic trajectory, KS is computed by the log likelihood that it occurred in the real fly walking data (illustrated in Fig. 4A and detailed in Methods and Materials). Briefly, we reduce the experimental data to 2 dimensions using principal components analysis (PCA), then fit a kernel density estimator (KDE) to the resulting distribution. We then project each bout of simulated walking onto this subspace and evaluate the KDE model to obtain a log probability density function estimate, which corresponds to the kinematic similarity of the simulated bout to real walking bouts in the data. Lower values of KS mean lower similarity to data. When we applied this method to bouts from experimental data, we found that the average KS of experimental data was -1.627. Thus, we use KS > −1.6 as a general threshold for evaluating the realism of joint angle trajectories. As expected, larger perturbations led to walking behaviors with lower KS, and KS was not strongly dependent on forward walking speed (Fig. 4B-C). For persistent stochastic perturbations, the model produced realistic joint angle trajectories (KS > −1.6) for all simulated perturbation strengths. Impulse perturbations appeared to be more challenging for the model, as they resulted in lower KS. From the time-series, we observed that impulse perturbations resulted in greater instantaneous deviation from the standard waveform than persistent stochastic perturbations of the same magnitude. This is due to the fact that impulse perturbations produce simultaneous changes to all legs and joints, whereas persistent stochastic perturbations are spread out in time and typically non-simultaneous.

Figure 4.

2.4 Effect of sensory and motor delays on walking

Temporal delays are inherent properties of sensorimotor control systems, but they are difficult to manipulate experimentally. Therefore, we used our model to investigate how changing sensory and motor delays affects locomotor robustness. We used measurements from the literature to estimate physiological delays for leg sensory and motor neurons in the Drosophila leg. Our estimate of sensory delay (5–15 ms) was based on the measured delay from spike initiation in a mechanosensory neuron in the Drosophila femur to the peak of an excitatory postsynaptic potential in a postsynaptic VNC neuron (Tuthill and Wilson, 2016a). Our estimate of motor delay (20–40 ms) was based on the time between spike initiation in a tibia motor neuron cell body to the onset of muscle force production, measured with a force probe (Azevedo et al., 2020).

Without external perturbations, the model produced realistic walking with arbitrary delays, since the controller can effectively compensate for large delays by using predictions of joint angles in the future. However, in the presence of external perturbations and high delay values, the model was unable to maintain realistic walking, since it could not respond rapidly enough to unexpected perturbations. Here we consider the composite effects of persistent stochastic perturbations with motor and sensor delays. Similar results for impulse perturbations are included in Appendix I.

We first fixed the sensory delay at 10 ms and measured the effect of varying motor delays. Examining the time-series data, we found that for low values of motor delay (10 ms, 20 ms), even stronger perturbations had almost no effect (Fig. 5A). However, at higher values of motor delay, the effects of the same perturbation became more pronounced, though the model still managed to recover after the perturbations ended. Over a range of perturbation strengths and walking speeds, the model maintained realistic walking (KS > −1.6) up to about 30 ms of motor delay (Fig. 5C). This value is consistent with motor delays measured from the fruit fly leg, where the time between motor neuron spiking to the onset of muscle force production is about 30 ms (Azevedo et al., 2020).

Figure 5.

Next, we fixed the motor delay at 30 ms and observed the effect of varying sensory delays. From time-series data, we found that even at low values of sensory delay (0 ms, 5 ms), the effects of perturbations on model output were significant (Fig. 5B). When sensory delay increased, the effects of the perturbation became more pronounced. Over a range of perturbation strengths and walking speeds, the model maintained realistic walking (KS > −1.6) up to about 10 ms of sensory delay (Fig. 5D). This range is consistent with experimental estimates of delays from mechanosensory neurons in the fly’s fly femur (Tuthill and Wilson, 2016a), although delays could be even longer for more distal sensors. For example, a sensory neuron at the tip of the tarsus would take ∼20 ms to travel 2 mm at an estimated conduction velocity of ∼0.3 m/s (Tuthill and Wilson, 2016a).

For any fixed set of values of motor and sensory delays, the KS of perturbed walking was not dependent on forward walking speed. In other words, slower walking did not improve the model’s ability to sustain walking during perturbations (Fig. 5C and D). We expand on this observation in the Discussion.

Last, we explored combinations of values for motor and sensory delays for fixed walking speed and perturbation strength. We observed that robust walking was not contingent on the specific values of motor and sensory delay, but rather the sum of these two values (Fig. 5E). That is, in order for the model to overcome a perturbation, the key parameter is the total delay between the perturbation onset and the motor response.

In summary, we used a virtual fly model of walking to perform in silico manipulations of sensory and motor delays, including delays exceeding physiologically-realistic values. When we tested the model across a range of delay values, we found that the model was able to maintain robust walking only within the physiological range. Outside of this range, the model could not maintain walking that was kinematically similar to real flies. This finding suggests that the fly locomotor control system operates close to the temporal limit at which it can detect and respond to external perturbations. Thus, we propose that sensory and motor delays establish a fundamental tradeoff between the speed and robustness of fly walking.

3 Discussion

In this paper, we developed a computational model that realistically imitates 3D joint kinematics of walking Drosophila and incorporates sensorimotor delay as a tunable parameter. We then used the model to establish a quantitative relationship between sensorimotor delays and locomotor robustness. We found that the model’s ability to maintain walking following external perturbations significantly degrades for delay values that exceed known physiological values, suggesting that these parameters are fundamental constraints on fly locomotion. The formulation of a modular, multi-layered model for locomotor control makes new experimentally-testable hypotheses about fly motor control and can also be applied to investigate limbed locomotion in other organisms. Future extensions of the model include inclusion of premotor neural circuits from the fly connectome (Lesser et al., 2023; Cheong et al., 2023) and biomechanical interactions between the limb and the environment (Lobato-Rios et al., 2022; Wang-Chen et al., 2023; Vaxenburg et al., 2024).

3.1 Fundamental constraints on locomotion imposed by sensory and motor delays

Sensory and motor delays are inextricable properties of animal locomotor systems. To study the impact of delays on locomotor control, we developed a hierarchical walking model with explicit inclusion of physiological delays as tunable parameters. Importantly, the model incorporates delay while preserving behavioral realism, thus laying the groundwork for future studies on the effect of delay on other aspects of locomotion and sensorimotor control.

Our model predicts that at the same perturbation magnitude, walking robustness decreases as delays increase. This could be experimentally tested by altering conduction velocities in the fly, for example by increasing or decreasing the ambient temperature (Banerjee et al., 2021). If a warmer ambient temperature decreases delays in the fly, but fly walking robustness remains the same in response to a fixed perturbation, this would indicate a stronger role for central control in walking than our modeling results suggest.

Tradeoffs between energy, performance, and robustness establish fundamental constraints on animal locomotion. It is energetically expensive to build and maintain muscles and neurons that operate with low delay (Sterling and Laughlin, 2015). From a performance perspective, long delays limit reaction speed and limit control robustness, which can ultimately impact animal survival (More and Donelan, 2018). A third consideration is locomotor robustness, or the capacity to detect and respond to external perturbations. In flies, the axons of leg motor neurons and some proprioceptor axons are among the largest diameter cables in the leg nerve (Phelps et al., 2021), suggesting that the speed of these systems are under selective pressure.

Our model suggests that fly walking operates in a middle ground with respect to speed and robustness. Specifically, our model maintains normal walking for up to about 30 ms of motor delay and 10 ms of sensory delay. These values are strikingly close to measured values from Drosophila: 30 ms from motor neuron spike to peak force production in a femur muscle (Azevedo et al., 2020) and 6 ms from a femoral mechanosensory neuron spike to the onset of a postsynaptic response in the VNC (Tuthill and Wilson, 2016a).

In our model, robust locomotion was constrained by the cumulative sensorimotor delay. This result could be experimentally validated by comparing how animals with different ratios of sensory to motor delays respond to perturbations. Alternatively, it may be possible to manipulate sensory vs. motor delays in a single animal, perhaps by altering the development of specific neurons or ensheathing glia (Kottmeier et al., 2020). If sensory and motor delays have significantly different effects on walking quality, then additional compensatory mechanisms for delays could play a larger role than we expect, such as prediction through sensory integration, mechanical feedback, or compensation through central control. A rich and related topic for future exploration is the interaction of delays with body size and behavioral ecology, as longer limbs, heavier bodies, and/or more mechanically complex tasks may alter the tradeoffs between energetics, speed of force generation, noise, and robust sensorimotor control (Labonte et al., 2024).

3.2 Role of proprioceptive feedback in fly walking

Our model provides insight into the role of proprioceptive feedback in fly walking, which remains an active area of research (Dallmann et al., 2021). Models of fly walking have generally ignored the role of feedback, relying instead on central pattern generators (Lobato-Rios et al., 2022; Szczecinski et al., 2018; Aminzare et al., 2018) or metachondral waves (Deangelis et al., 2019) to model kinematics. This assumption is in line with experimental work in cockroaches, which suggests that these fast-running insects rely on central control and mechanical feedback, particularly at high-speeds (Couzin-Fuchs et al., 2015; Ayali et al., 2015b). In contrast, studies in stick insects have shown that these slow-walking animals are highly dependent on proprioceptive feedback for leg coordination during walking (Bässler, 1977; Ayali et al., 2015a; Schilling et al., 2013).

Silencing mechanosensory chordotonal neurons alters step kinematics in walking Drosophila (Mendes et al., 2013; Pratt et al., 2024). However, the role of leg proprioception in overcoming external perturbations has not previously been studied in flies. In our model, which includes no limb compliance or other biomechanical adaptations, the fly effectively overcomes perturbations using only proprioceptive feedback. The need for proprioception to compensate for perturbations may be reduced when biomechanical mechanisms are present. Nonetheless, we hypothesize that removing proprioceptive feedback would impair an insect’s ability to sustain locomotion following external perturbations.

3.3 Predictive control is critical for responding to perturbations due to motor delay

Sustaining realistic walking kinematics in the presence of perturbations and motor delays is a challenging task for any model. In our model, this design criterion motivated us to develop a compensatory prediction in the optimal controller formulation, based on Stenberg et al. (2022). The controller compensates for the known motor delay by predicting future dynamics. We found that the model was quite sensitive to the prediction horizon, i.e., how far into the future the controller predicts. The model works best when the controller’s prediction horizon matches the motor delay — we thus used matched values in all simulations. We experimented with altering the prediction horizon to be less than the motor delay, with catastrophic consequences: the model ceased walking and mostly produced noise (not shown). From this observation, we deduce that future predictions are crucial in compensating for motor delays in our walking model, in agreement with previous theoretical work on predictive or “forward” models for sensorimotor control (Desmurget and Grafton, 2000; Li et al., 2023). We propose that fly motor circuits may encode predictions of future joint positions, so the fly may generate motor commands that account for motor neuron and muscle delays. Consistent with this hypothesis, Dallmann et al. (2023) recently found that descending motor commands from the brain excite GABAergic interneurons in the VNC that inhibit velocity-encoding proprioceptors. Thus, some proprioceptive feedback signals from the fly leg are predictively suppressed during self-generated movements.

3.4 Layered model produces robust walking and facilitates local control

One key finding of our model is that robust walking in the presence of external perturbations can emerge from a local controller in combination with a trajectory generator that is trained only on perturbation-free walking. We did not model adaptation because compensation emerges as a property of a tuned controller with feedback. The success of this simple, layered model suggests that fast, robust locomotion could be maintained locally, requiring minimal plasticity within central circuits.

Layering is a familiar concept in network architecture and hardware/software stacks (Chiang et al., 2007) and has recently found applications as a modeling framework in biology and neuroscience (Doyle and Csete, 2011; Nakahira et al., 2021). In our model, the complex problem of multi-legged locomotion is broken down into three sub-problems (delayed control, joint kinematic generation, inter-leg coordination), each of which is delegated to a different layer.

Each layer provides an abstracted interface for the layer below it, which can also be thought of as a special case of modularity. The separation of function between layers reinforces that robust locomotion can be produced by local (i.e. per-leg) control signals. Per-leg modules could also be extended to control more legs, generalizing to locomotor control in isopods and millipedes. Instead of redesigning the overall coordination and locomotor control circuits, basic duplication and subsequent fine tuning may be sufficient.

Beyond walking, our layered model framework could be applied to other animals with distinct locomotor strategies, including flying, swimming, slithering, and digging. The key ingredients required for the model are: (1) a functional inter-limb coordinator, (2) sufficient 3D kinematics data to train a trajectory-generator layer, and (3) a controller that adequately tracks the trajectory for some dynamical model of the animal. The dynamical model may be linearized (as we did here) and controlled with a standard controller; other approaches include successive (per-timestep) linearization with linear model predictive control (MPC) (Berberich et al., 2022), or fully nonlinear control techniques (Slotine and Li, 1991).

The layered model approach also has potential applications for bio-mimetic robotic locomotion. For simplicity, we used a linearized link-and-joint dynamical model. However, one could also replace this dynamical model to control a hexapod robot such as in Goldsmith et al. (2019). Due to its modular nature, the other modules (trajectory generator, inter-leg coordinator) would remain unchanged and the resulting model could, in theory, generate 3D kinematics for bio-mimetic walking on a robotic platform. Overall, the model can be thought of as performing a layered implementation of imitation learning, which is a popular technique in robotics.

3.5 Towards biomechanical and neural realism

The goal of our model was to produce realistic 3D joint kinematics while incorporating sensory and motor delays. To achieve this, the model contains several physiological simplifications. First, our dynamics model did not allow dynamical coupling between legs through the mechanics of the body, as the legs are only coupled neurally through the phase coordinator. However, in real bodies, the legs are also dynamically coupled through the body and its weight distribution over the legs (Dallmann et al., 2017). Our model also did not consider explicit leg-ground contact interactions. Rather, interactions with the ground were implicitly taken into account by the trajectories learned by the neural network, though they were not made explicit in the dynamics. Our goal was to mimic the kinematic trajectory, a problem known in robotics as motion control. However, including ground contact interactions would require computing ground contact forces, which are currently unavailable in the kinematics dataset we used.

In order to model ground contact forces and joint torques, force-based learning would need be incorporated into the trajectory generator, requiring a new dataset or extrapolating force values (e.g. via a physics-based model). Further, the controller would need to be reformulated to use impedance control or hybrid control techniques (Buss et al., 2002; Arevalo and Garcia, 2012; Sciavicco and Siciliano, 2012). These approaches are commonplace in robotics to deal with dynamics control problems, which concern both kinematic trajectories and external contact forces. The inclusion of explicit leg-ground contact interactions would also make it harder for the model to recover when perturbed, because perturbations during walking often occur upon contact with the ground (e.g. the ground is slippery or bumpy).

A promising avenue for future investigation is integration of our controller architecture with a virtual physics model (Lobato-Rios et al., 2022; Wang-Chen et al., 2023; Vaxenburg et al., 2024), which would facilitate incorporation of dynamical coupling between legs, as well as leg-ground contact interactions. The inclusion of these features may require additional coordination between the legs, which might decrease allowable values of sensory and motor delay.

Another step toward biological realism is the incorporation of explicit proprioceptors, muscles, tendons, and other biomechanical aspects of the exoskeleton. The proprioceptive neurons in the femoral chordotonal organ of each fly leg encode angles and angular derivatives (Mamiya et al., 2018). Additional proprioceptive feedback is provided by hair plate sensory neurons (limit detectors) and campaniform sensilla (load sensors), which are distributed across each leg (Tuthill and Wilson, 2016b). Thus, our use of joint state θ, , likely underestimates the resolution of proprioceptive feedback to the fly motor system. We anticipate that the increased sensory resolution from more detailed proprioceptor models and the stability from mechanical compliance of limbs in a more detailed biomechanical model would make the system easier to control and increase the allowable range of delay parameters. Conversely, we expect that modeling the nonlinearity and noise inherent to biological sensors and actuators may decrease the allowable range of delay parameters.

A further step towards biological realism would be to constrain the trajectory generator and optimal controller using patterns of synaptic connectivity within sensorimotor circuits of the fly VNC. This is now feasible using recent connectomes of the Drosophila VNC (Azevedo et al., 2022; Takemura et al., 2023). However, many challenges remain for connectome-constrained models, because many important physiological parameters are still unknown.

Our layered approach could be a useful framework for learning principles from closed loop models, despite having incomplete information about the biological system. For instance, the trajectory generator model could be replaced by a combination of a connectome simulation and artificial neural networks, with the constraint on the artificial neural network dynamics coming from the complete circuit generating robust walking in simulation. Future work could also incorporate more detailed delay structures based on compartmental models of neurons from the connectome. With increasing realism, such closed-loop models provide a promising means to investigate how complex neural circuits interact with proprioceptive feedback to control robust locomotor behaviors.

4 Methods and Materials

4.1 Tracking joint angles of D. melanogaster walking in 3D

We obtained fruit fly D. melanogaster walking kinematics data following the procedure previously described in Karashchuk et al. (2021). Briefly, a fly was tethered to a tungsten wire and positioned on a frictionless spherical treadmill ball suspended on compressed air. Six cameras captured the movement of all of the fly’s legs. Using Anipose, we tracked 30 keypoints on the fly at 300 Hz, which are the following 5 points on each of the 6 legs: body-coxa, coxa-femur, femur-tibia, and tibia-tarsus joints, as well as the tip of the tarsus.

To fit the model described in this paper, we extracted a subset of the tracking data when the fly walking, as opposed to non-walking behaviors including standing, grooming, etc. To isolate bouts of walking, we used the behavior classifier described in Karashchuk et al. (2021). We further selected walking bouts of at least 0.5 seconds (150 video frames) in length, and where the femur-tibia flexion angle of the left front leg had a range of at least 30 degrees.

In total, our dataset consisted of 3473 walking bouts from 45 flies total. The average length of a walking bout was 0.877 seconds (263 frames), with 3049.7 seconds of walking total (914,909 frames).

4.2 Inter-leg phase coordinator

We model the coordination between legs as phase-coupled Kuramoto oscillators (Strogatz, 2000), where the frequency of each oscillator is driven by the trajectory generator described in the next section.

Specifically, the phase for a leg i is ϕ_i and evolves according to its derivative

where F_i is the trajectory generator function for leg i, α is the coupling strength, and is the steady-state phase offset between legs i and j.

We model the coupling across the legs as all-to-all coupling, with coupling strength α = 6.5. We found this coupling strength best reproduced the phase coupling distributions from the real data (as shown in Appendix F).

We estimate from the walking data by taking the circular mean over phase differences of pairs the legs during walking bouts. We find that the phase offset across legs is not modulated across walking speeds in our dataset (see Appendix B), so we model as a single constant independent of speed.

4.3 Trajectory generator

A trajectory generator model was formulated for each of the six legs and fit separately to fly walking data tracked during tethered walking without any external perturbations.

4.3.1 Model formulation

We formulate the trajectory generator as the function

where θ is a vector of joint angles, s a vector of joint angle derivatives, v is the desired walking speed and direction, and ϕ is the phase of the leg. Initially, we explored using the trajectory generator to directly output angles and angular velocities ; however, we found that more realistic (i.e. similar to data) trajectories were produced when we used the trajectory generator to output angular acceleration , which we integrated to produce the desired angle and angular velocities. Note that v is not communicated to or from the layers above and below; instead, we consider walking speed and direction to be given as commands descending from the brain.

To compute a trajectory given an input v and an initial ϕ, θ, and , we integrate the function F numerically using the midpoint method (Lotkin, 1956). Following methods from Holden et al. (2017) and Zhang et al. (2018), we represent the function F as a multilayer perceptron neural network with 2 hidden layers of with 512 units each. We use ELU (Clevert et al., 2015) as our nonlinearity. In total, the multilayer perceptron has 274,437 parameters for T1 legs and 272,388 parameters for T2 and T3 legs, with the slight difference in parameters due to the different number of joint angles (dimension of θ) modeled for a given leg.

4.3.2 Training data

To train the multilayer perceptron network used to represent F, we used the fly walking data, tracked as described in the section above. The training data consists of joint angles θ, computed and , and walking velocity v. We estimated the walking cycle phase ϕ using a Hilbert transform over the femur-tibia flexion angle for T1 legs, femur rotation angle for T2 legs, and coxa-femur angle for T3 legs. For each phase, we filtered the corresponding angle using a first-order Butterworth bandpass filter with 3 Hz and 60 Hz as critical frequencies, using the scipy library (Virtanen et al., 2020). Then, we applied a Hilbert transform to each angle to obtain a complex waveform. We estimate the walking cycle phase from each complex waveform by estimating the angle of each point in the waveform.

4.3.3 Training procedure

Training the neural network representing F from data was performed in two steps, minimizing its error in predicting one time step, then minimizing its error in predicting a short time trajectory.

In the first step, we minimized the error of F for predicting over one time step, given the corresponding from the training data. We minimized the mean squared error of the prediction, normalized by the variance for each dimension. We trained our network for 300 iterations over the full training data using a batch size of 2500 training samples, using gradient descent with the Adam algorithm (Kingma and Ba, 2017). To ensure a robust function at this step, we applied dropout to a random 5% of the hidden units (Srivastava et al., 2014). We standardized the input and output training data to the multilayer perceptron so that has a mean of 0 and standard deviation of 1.

In the second step, we minimize the error of F for predicting a trajectory θ when numerically integrated over a short time horizon in the future. Specifically, we integrate F over T = 60 steps given initial conditions to produce an estimated trajectory of . Here we minimized the loss:

using gradient descent with the Adam algorithm (Kingma and Ba, 2017). During training, we clip gradients to a norm of 10 to stabilize training.

The training was implemented using Tensorflow (Abadi et al., 2015) running on a computer with NVIDIA GeForce RTX 2070 GPU and AMD Ryzen Threadripper 1920X 12-Core Processor.

4.4 Leg dynamics and optimal controller formulation

All techniques used for dynamics formulation are standard tools from control theory. We begin with a link-and-joint model of the fly leg, as shown in Fig. 1C. For simplicity, we only model joints that exhibit large ranges of movement during naturalistic walking and turning. For instance, varying femur rotation is important to the movements of the middle legs, but the front legs exhibit near-constant femur rotation (Karashchuk et al., 2021); thus, a femur rotation joint is included for the middle and hind legs only. The joints included for each leg is shown in Table 1.

Table 1:

We write the Denavit-Hartenberg (DH) table of the leg model and use this to systematically derive the Euler-Langrange matrix equations of motion:

where τ is the vector of joint torques; θ, , and are vectors of joint angles, angular velocity, and angular acceleration; M, C, B, are the inertia, Coriolis, and friction matrices, and g is the gravity vector.

Let us define the state to be the angles and angular derivatives and input to be torques τ. We next rearrange (4) into the form, so that

To linearize this system, we choose equilibrium values and , such that . In particular, is taken to be the average joint angles for each leg joint computed from the data. It follows that and gives the desired equilibrium.

We linearized about this equilibrium point, which leads to an equation of the following form:

where A and B are the Jacobians with respect to q and τ, respectively. In other words,

In our code, we use the SymPyBotics toolbox (Sousa, 2013) to obtain symbolic equations for the quantities in (4), then numerically compute Jacobian values.

Next, we rewrite the system in discrete time using a sampling interval T, which is typically chosen to be an integer multiple of sampling interval from the data (the tracking data was acquired every T = 1/300 seconds). In our controller simulations, we use T = 1/600 seconds. The discretized dynamics are thus written as:

where A = I + A_cT and B = B_cT.

Finally, we perform a coordinate shift to error dynamics. This allows us to apply standard control techniques for trajectory tracking. We define the tracking error to be y = q − q_d, where q_d is the desired state, and this error obeys the following dynamics:

where w is the external perturbation. w_traj represents the effect of constantly changing trajectories — for example, if the desired trajectory at the current time-step is some value a, and the desired trajectory at the next time-step is some other value b, then this is equivalent to introducing a perturbation of b − a. Error y is of dimension n_y (8 for front legs, and 6 for the other legs), and input u is of dimension .

To include motor delay (known as actuation delay in controls literature) and sensory delay, we make use of an augmented state formulations as introduced in Stenberg et al. (2022). Let the motor delay be d_motor steps and let the sensory delay be d_sense steps. Our augmented state vector z(t) is then written as follows:

where y is the error from (10); f contains predicted future errors up to d_motor time steps in the future of dimension d_motor * n_y; a contains motor signals from up to d_motor steps ago of dimension d_motor * n_u; s contains sensing signals from up to d_sense steps ago of dimension d_sense * n_y; and g contains information about future trajectory effects w_traj up to d_motor time steps in the future of dimension d_motor * n_y.

Finally, we write the overall system in the form of

where F, G, and H are formulated using A and B according to techniques described in Stenberg et al. (2022), and w_aug(t) contains the perturbations from (10), appropriately rearranged and zero-padded.

We note that u represents the motor signal, which must be delayed for some amount of timesteps (via a) before affecting the system; similarly, the sensory signal is delayed for several timesteps (via s) before reaching the controller via r.

To achieve effective trajectory tracking, we seek a control law under which y remains small. This can be achieved using the Linear Quadratic Gaussian (LQG) controller (Å ström and Murray, 2021) The controller is governed by the following equations:

where is the estimate of the state, estimated via a steady-state Kalman filter; L and K are the optimal observer and controller matrices, respectively, synthesized via discrete algebraic Riccati equations.

4.5 Generate joint trajectories of the complete model with perturbations

The full model integrates all of the modules to generate trajectories of joint angles over time. The phase coordinator and trajectory generators combine to compute the desired joint angles θ_i for each leg i, and the controller implements them constrained by the dynamics. We run the full model at 600 Hz.

At each time step, we first update the phases of the legs ϕ_i based on the phase coupling equation above. Every 2 timesteps, we update the target joint angles using the trajectory generator model. We run the trajectory generator F for d_motor time steps on its own, by continually integrating its output. This future trajectory forms the basis of f, the predicted future errors used to guide the controller.

Each time step, we run the controller and dynamics model to control the torque τ so that the joint angles θ go towards the target joint angles If there is a disturbance, we apply it to the joint angles and derivatives at this point. Every 8 timesteps, we set , so that the trajectory generator predicts an intended trajectory in line with the current state.

For the perturbation numerical experiments described in the Results, we ran all simulations for 1800 timesteps at 600 Hz. For persistent stochastic perturbations, we applied the perturbation from 600th timestep to 1200th timestep. For impulse perturbations, we applied a single strong perturbation at 600th timestep.

We ran all our simulations on an Intel(R) Core(TM) i9-9940X CPU. We used GNU Parallel (Tange, 2011) to run simulations on multiple cores simultaneously.

4.6 Computing Kinematic Similarity (KS) by quantifying likelihood of walking kinematics relative to ground truth

We followed a multi-step procedure in order to quantify the likelihood of the simulated walking kinematics relative to observed distribution of walking kinematics. This procedure is schematized in Fig. 4A.

To fit our likelihood model, we first performed Principal Components Analysis (PCA) to reduce the tracked joint angles from our 3D kinematics data to 2 dimensions per frame. To account for the circular nature of rotation angles, we performed PCA on the combination of sines and cosines of each angle. Next, we used a Gaussian Kernel Density Estimation (KDE) to estimate the probability density function of the principal components. Thus, we obtained a likelihood model for joint angle kinematics at each frame.

We chose 2 dimensions for PCA for two key reasons. First, these 2 dimensions alone accounted for a large portion of the variance in the data (43.1% total, with 30.5% for first component and 12.7% for second component). There was a big drop in variance explained from the first to the second component, but no sudden drop in the next 10 components. Second, the KDE procedure only works effectively in low-dimensional spaces, and the minimal number of dimensions needed to obtain circular dynamics for walking is 2.

We run our model described above to produce simulated joint angle trajectories. To estimate a likelihood of a simulated set of angles, we first projected them onto the same principal components identified from the observed kinematics. Then, we use the KDE model to estimate the logarithmic probability density function (log PDF) for each frame during the perturbation; we refer to this as the kinematic similarity (KS). For persistent stochastic perturbations, we estimated mean KS during the perturbation, from 600th timestep to 1200th timestep. For impulse perturbations, we estimated mean KS in the transient recovery process, which we estimated to be from 610th timestep to the 800th timestep.

4.7 Visualization of joint movement trajectories

For Video 1, we visualized the simulated and real joint movements using the biomechanical fly body model from Vaxenburg et al. (2024). For each frame, we ran an inverse kinematics optimization over the model angles to match the simulated or real joint positions of the fly. We did not simulate realistic physics of the fly legs and their interactions with the ground in these visualizations. In our data, the fly thorax was fixed and the wings were removed, so in these visualizations, we also fixed every other degree-of-freedom in the fly model besides the 6 legs. For the remaining videos, we visualized the joints as ball-and-stick models using matplotlib (Hunter, 2007).

4.8 Code and data availability

Code is available at https://github.com/lambdaloop/layered-walking/. The data used in this study will be released publicly upon publication and privately by request.

Contributions

LK, JSL, SLB, JCT, and BWB conceived of the study. LK and JSL developed the model, ran simulations, and performed the analyses. GMC and SW-B collected data of fly walking. LK and GMC visualized the walking fly videos. LK, JSL, JCT, and BWB interpreted the results and wrote the paper. All authors reviewed and edited the paper.

Acknowledgements

We thank J. Doyle for initial discussions on layered feedback models, N. Csomay-Shanklin for initial discussions on deriving models of legged locomotion, A. Nair for discussions around coupled oscillator models, A. Schwarze for discussions around sensorimotor feedback loops, E. Dickinson for technical assistance, and members of the Tuthill and Brunton labs for helpful discussions and feedback on this manuscript. We are grateful to R. Vaxenburg, I. Siwanowicz, and S. Turaga for help adapting their virtual fly model for visualizing real and simulated fly walking.

This work has been funded by a National Science Foundation Graduate Research Fellowship to LK; a NSERC Postgraduate Scholarship (NSERC PGSD3-557385-2021) to JSL; a Searle Scholar Award, a Klingenstein-Simons Fellowship, a Pew Biomedical Scholar Award, a McKnight Scholar Award, a Sloan Research Fellowship, the New York Stem Cell Foundation, and NIH R01NS102333 to JCT; JCT is a New York Stem Cell Foundation – Robertson Investigator; Air Force Office of Scientific Research award FA9550-19-1-0386 and the Richard & Joan Komen University Chair to BWB.

Appendix

B Velocity and phase distributions in data

Figure 6.

C Angle vs. phase plots for all legs and joints

These plots are supplementary to Fig. 2B.

Figure 7.

D Angular velocity vs. phase plots for all legs and joints

These plots are supplementary to Fig. 2D.

Figure 8.

E Differences between model and data

These plots are supplementary to Fig. 2E.

Figure 9.

F Phase coupling within and across legs

We quantify the model’s capacity to simulate leg kinematics by comparing the phase coupling of simulated joint kinematics with real walking flies, both within and across legs. A phase difference probability density peak of close to zero between two joint angles means that they are strongly coupled, whereas a broader peak means they are weakly coupled.

To assess the effectiveness of the model’s inter-leg coordination, we compared phase coupling across legs for simulated vs. real walking. We found that the model reproduced patterns in the data, with some exceptions. Here, a phase difference of zeros means the two legs are synchronized, whereas a phase difference of π means the two legs move in anti-phase. Fly data showed largely tripod coordination, so that L1-R2-L3 are synchronized with each other and also anti-phase with R1-L2-R3. The model showed coupling properties that were qualitatively similar to coupling properties of the data, with moderate variations in phase difference (i.e. peak location) and synchronization strength (i.e. peak height). Some minor differences may also be due to violations of our assumption of constant phase offset across forward walking speeds. This assumption holds true over a large range of speeds, but breaks down at the lowest walking speeds (see Appendix B). In the model, the phase coordinator is responsible for representing inter-leg coupling; these results suggest that the underlying Kuramoto oscillator reproduces naturalistic coupling patterns.

Next, we examined within-leg coupling in data, where phase synchronizations between different joints on the same leg can be both strong (e.g. L3 coxa-femur flexion to femur-tibia flexion) and weak (e.g. R2 femur-tibia flexion to femur rotation). Despite the mirror symmetry of the fly, within-leg coupling is not necessarily symmetric in the data — for instance, compare L2 and R2 coxa-femur flexion. All within-leg synchronizations (i.e. peaks) in data were reproduced by the model; however, the model generally exhibited stronger synchronization (i.e. higher peaks) than the data. In some cases (e.g. L3 femur rotation to femur-tibia flexion), the model exhibited synchronizations that are not present in the data. Overall, the model generally exhibits stronger within-leg coupling when compared to the data.

The trajectory generator is responsible for reproducing the per-leg inter-joint coupling, so these results suggest that the neural network learned stronger coupling values than are present in the data. This result is consistent with the time-series comparisons, where we observed that the model produced more regular joint angle trajectories than the fly (Fig. 2A).

Figure 10.

G Estimate of plausible stochastic perturbation values

Here we derive an estimate of the maximum plausible stochastic perturbation magnitude corresponding to a sudden gust of wind. Wind force experienced by a given leg is

where ρ, v, and A are the density of air, wind speed, and perpendicular area (to wind direction) of the leg. The resulting acceleration is

where m is the mass of the fly.

We assume, generously, that the leg experiences the same acceleration. Since the controller operates at a frequency of f, we can approximate the angular velocity (i.e. perturbation) ω generated by the slip acceleration in one timestep for one joint on the leg as

where r is the length of the leg.

Substituting in values used in our model (m = 0.7e-6 kg, f = 600 Hz, r = 0.0015 m), and let ρ = 1.2 kg/m³ (standard value), v = 6 m/s (corresponding to a strong 21 km/h wind), and A = 0.15e-6 m² (corresponding to a 1.5mm by 0.1mm leg), we obtain a maximum joint perturbation magnitude of ω = 5.1 rad/s.

H Estimate of plausible impulse perturbation values

We estimate the maximum plausible impulse perturbation corresponding to a sudden slip. Assume that a leg steps at an angle of θ (from the vertical) onto a flat, frictionless ground. The resulting horizontal acceleration would be

where g is acceleration due to gravity.

Assume the controller operates at a frequency of f. We can approximate the angular velocity (i.e. perturbation) ω generated by the slip acceleration in one timestep for a joint on the leg as

where r is the length of the leg. Plugging in values of g = 9.81 m/s² (standard value), f = 600 Hz (value used in the model), r = 0.0015 m (value used in the model), and θ = 45 degrees (quite a large slip angle), we obtain a perturbation value of n ω = 10.9 rad/s.

I Effect of sensory and motor delays on walking under impulse perturbations

The results shown in Fig. 5 correspond to persistent stochastic perturbations. Here, we present similar results for impulse perturbations here. Maximum values of motor and sensory delay values remain reasonably close to known values in physiology.

Figure 11.