Variational design of sensory feedback for powerstroke-recovery systems
Zhuojun Yu
zhuojun.yu@case.edu
Department of Mathematics, Applied Mathematics, and Statistics, Case Western Reserve University, Cleveland, OH 44106, USA
Peter J. Thomas
pjthomas@case.edu
Department of Mathematics, Applied Mathematics, and Statistics, Department of Biology, Department of Electrical, Control and Systems Engineering, Case Western Reserve University, Cleveland, OH 44106, USA
Although the raison d’etre of the brain is the survival of the body, there are relatively few theoretical studies of closed-loop rhythmic motor control systems. In this paper we provide a unified framework, based on variational analysis, for investigating the dual goals of performance and robustness in powerstroke-recovery systems. To demonstrate our variational method, we augment two previously published closed-loop motor control models by equipping each model with a performance measure based on the rate of progress of the system relative to a spatially extended external substrate – such as a long strip of seaweed for a feeding task, or progress relative to the ground for a locomotor task. The sensitivity measure quantifies the ability of the system to maintain performance in response to external perturbations, such as an applied load. Motivated by a search for optimal design principles for feedback control achieving the complementary requirements of efficiency and robustness, we discuss the performance-sensitivity patterns of the systems featuring different sensory feedback architectures. In a paradigmatic half-center oscillator (HCO)-motor system, we observe that the excitation-inhibition property of feedback mechanisms determines the sensitivity pattern while the activation-inactivation property determines the performance pattern. Moreover, we show that the nonlinearity of the sigmoid activation of feedback signals allows the existence of optimal combinations of performance and sensitivity. In a detailed hindlimb locomotor system, we find that a force-dependent feedback can simultaneously optimize both performance and robustness, while length-dependent feedback variations result in significant performance-versus-sensitivity tradeoffs. Thus, this work provides an analytical framework for studying feedback control of oscillations in nonlinear dynamical systems, leading to several insights that have the potential to inform the design of control or rehabilitation systems.
Keywords: Sensory feedback, Closed-loop control, Central pattern generator, Power stroke, Robustness, Efficiency
1 Introduction
Physiological systems underlying vital behaviors such as breathing, walking, crawling, and feeding, must generate motor rhythms that are not only efficient, but also robust against changes in operating conditions. Although central neural circuits have been shown to be capable of producing rhythmic motor outputs in isolation from the periphery (Brown,, 1911, 1914; Harris-Warrick and Cohen,, 1985; Pearson,, 1985; Smith et al.,, 1991), the role of sensory feedback should not be underestimated. Sensory feedback can play a crucial role in stabilizing motor activity in response to unexpected conditions. For example, modeling work suggests that walking movements can be stably restored after spinal cord injury by enhancing the strengths of the afferent feedback pathways to the spinal central pattern generator (CPG) (Markin et al.,, 2010; Spardy et al.,, 2011). Feedback control can also improve the performance and efficiency of movements. For instance, in a model of feeding motor patterns in the marine mollusk Aplysia californica, seaweed intake can be increased by strengthening the gain of sensory feedback to a specific motor neural pool (Wang et al.,, 2022).
We are interested in understanding how sensory feedback contributes to control and stabilization within a specific class of rhythmic motor behaviors, namely, behaviors in which an animal (or robot) repeatedly engages and disengages with the outside world. We refer to the phase of the motion during which the animal is in contact with an external substrate as the power stroke, and the component during which the animal is disengaged as the recovery phase. The decomposition of a repetitive movement into powerstroke and recovery applies naturally to many motor control systems, including locomotion (Jahn and Votta,, 1972) and swallowing (Shaw et al.,, 2015); a similar dynamical structure also appears in mechanical stick-slip systems (Galvanetto and Bishop,, 1999) as well as abstract two-stroke relaxation oscillators (Jelbart and Wechselberger,, 2020). In the motor control context, when the animal is in contact with an external substrate or load opposing the motion, we say the animal makes “progress” (food is consumed, distance is traveled, oxygen is absorbed) relative to the outside world. During the recovery phase, the animal disconnects from the external component, and repositions relative to the substrate in order to prepare for the next power stroke. Consider, for example, the ingestive behavior of Aplysia (Shaw et al.,, 2015; Lyttle et al.,, 2017; Wang et al.,, 2022). When the animal’s grasper is closed on a stipe of seaweed, it drags the food into the buccal cavity; meanwhile, the food applies a mechanical load on the grasper. Then the grasper opens, releasing its grip on the food. The grasper moves in the absence of the force exerted by the seaweed and returns to the original position to begin the next swallowing cycle.
In this paper, we present a novel analysis of feedback control for powerstroke-recovery systems. To quantitatively evaluate the behavior of a system controlled by different feedback mechanisms, we measure the sensitivity (or robustness) and performance (or efficiency). The complementary objectives of sensitivity and performance have been studied in a variety of motor control systems, from both empirical and theoretical perspectives (Lee and Tomizuka,, 1996; Yao et al.,, 1997; Ronsse et al.,, 2008; Hutter et al.,, 2014; Lyttle et al.,, 2017; Sharbafi et al.,, 2020; Mo et al.,, 2023). There are a myriad of ways to interpret performance and robustness used by engineers, biologists, neuroscientists, and applied mathematicians. Here we define the performance of a powerstroke-recovery system to be the total progress divided by the period of the rhythm (i.e. the average rate of progress), and the sensitivity to be the ability of the system to maintain performance in response to some perturbations on the external condition (i.e. the derivative of the performance with respect to the perturbation parameter). As a step towards first-principles–based design of sensory feedback mechanisms, we aim to understand what aspects of sensory feedback contribute to the coexistence of high performance and low sensitivity.
The ubiquity of powerstroke-recovery systems, and the importance of the dual goals of robustness and efficiency, motivate us to develop analytical tools for systematically studying both quantities simultaneously. In this work we apply mathematical tools based on variational analysis to evaluate the two objectives applicable for any powerstroke-recovery system. The key quantities in our analysis are the infinitesimal shape response curve (iSRC) and local timing response curve (lTRC) recently established by Wang et al., (2021, 2022) and generalized in Yu and Thomas, (2022); Yu et al., (2023). The iSRC describes, to first order, the distortion of an oscillator trajectory – the shape response – under a sustained perturbation. In contrast, the lTRC captures the effect of the perturbation on the timing of the oscillator trajectory within any defined segments of the trajectory (such as the powerstroke and recovery phases). Both the iSRC and the lTRC complement the more widely known infinitesimal phase response curve (iPRC), which quantifies the effect of a transient perturbation on the global limit cycle timing (Brown et al.,, 2004; Izhikevich and Ermentrout,, 2008; Ermentrout and Terman,, 2010; Schultheiss et al.,, 2012; Zhang and Lewis,, 2013).
Based on the iSRC and lTRC approach, we propose a general framework to investigate robust and efficient control through diverse feedback architectures. We apply our method to two neuro-mechanical models, each possessing a natural powerstroke-recovery structure but different in their levels of details and perturbations. The first model is based on an abstract CPG-feedback-biomechanics system introduced in Yu and Thomas, (2021) which studied the relative contributions of feedforward and feedback control (an idea going back to Kuo, (2002)). We extend this model to incorporate an externally applied load, enabling us to define quantitative measures of both performance and robustness, as the system alternately grasps and releases the external substrate. The second model, due to Markin et al. (Markin et al.,, 2010; Spardy et al.,, 2011), represents a locomotor system with a single-joint limb, and features more detailed CPG circuitry as well as more realistic afferent feedback pathways. We modify the Markin model so that the limb “walks” up an “incline”; modifying the angle of the incline introduces a parametric perturbation that allows us to define performance and robustness.
The activity of sensory feedback pathways is difficult to measure in many physiological systems; for this reason sensory feedback is the “missing link” for understanding the design and function of many biological motor systems. Specifically, in many experimental biological systems, the dynamics of the isolated CPG and the form of the muscle activation in response to descending motor signals is well characterized, while the precise form of the sensory feedback remains unkown. From a practical perspective, descending motor signals are generally carried by large-diameter axons from which it is easier to record high-quality (high signal-to-noise) traces, relative to ascending sensory signals, which are generally carried by much smaller diameter axons with poor signal-to-noise properties.
Given a particular specification of central neural circuit, descending output, and biomechanical response elements, but with the precise form of sensory feedback unknown, we can think of the pursuit of performance and/or robustness as an optimization (or dual optimization) problem on the space of sensory feedback functions. However, this is an infinite dimensional space of inputs to a highly nonlinear system. That is, the mapping from sensory feedback function to system trajectories to system performance and/or robustness is highly nonlinear and possibly non-convex. It may have multiple nonequivalent optima. Therefore, to restrict the problem to a manageable scope, in the specific cases studied here we restrict attention to sensory feedback functions with a prescribed, but plausible form, such as a sigmoid specified by a threshold and slope parameters, or multiple channels with different relative gain parameters, and study the restricted optimization problem there. Moreover, for simplicity, in the models we consider here, we restrict attention to sensory feedback that either monotonically increases or monotonically decreases with variables such as muscle length, tension, and/or velocity.
Our analysis of these two systems – one more abstract and the other more realistic – illustrates a technical framework for studying performance and sensitivity of powerstroke-recovery motor systems, leading to several insights that have the potential to inform the design of control or rehabilitation systems. For example, (i) the excitation-inhibition property of feedback signals determines the sensitivity pattern while the activation-inactivation property determines the performance pattern; (ii) the strong nonlinearity of feedback activation with respect to biomechanical variables may contribute to achievable performance-sensitivity optima; (iii) force-dependent feedback can prevent the performance/robustness tradeoffs commonly occuring with the length-dependent feedback. These findings may yield important information for future work modeling biphasic rhythm generation, in that they provide insights that could guide the design of feedback systems to accomplish well-balanced efficient and robust powerstroke-recovery activities in biological and robotic experiments.
The rest of this paper is organized as follows. We present the mathematical formulation of the control problem we consider, together with the variational analysis approach we provide, in Section 2. The two neuromechanical models, their feedback architectures, and their performance-sensitivity patterns are discussed in Section 3 and Section 4, respectively. Finally we summarize the broad framework as well as the main observations and insights that we obtain from the models in Section 5, where we also discuss limitations, connections to previous literature, and possible implications of our results for biology and engineering as well as future directions.
2 Mathematical formulation
The motor systems we consider integrate central neural circuitry, biomechanics, and sensory inputs from the periphery to form a closed-loop control system. The model systems we study fall within the following general framework (Shaw et al.,, 2015; Lyttle et al.,, 2017; Yu and Thomas,, 2021):
(1) |
Here, and are vectors representing the neural activity variables and mechanical state variables, respectively; the vector field represents the intrinsic neural dynamics of the central pattern generator when isolated from the rest of the body; captures the biomechanical dynamics driven by the central inputs; carries the sensory feedback from the periphery, which modulates the neural dynamics; is an externally applied load to the mechanical variables controlled by a load parameter .
In the powerstroke-recovery systems, we assume the load interacts with the mechanics only during the powerstroke phase. The portion of the trajectory comprising the powerstroke phase is specified separately for each model. For the purposes of this paper we assume the vector fields are sufficiently smooth (e.g. twice differentiable) except at a finite number of transition surfaces – for example at points marking the powerstroke-recovery transitions.
For many naturally occurring control systems, the mechanical variables may include both the position and velocity of different body components, as well as muscle activation variables. The sensory feedback function may be difficult to ascertain experimentally. For example, the feedback could have an excitatory effect on the neural dynamics, or an inhibitory effect, or a mixture at different points within a single movement; it may depend not only on neural outputs but also the length, velocity, or tension of the mechanical components; it may arise from multiple channels each with different gain. Given the broad varieties of the possible feedback functions, we restrict the scope of our investigation to some biologically plausible forms for the specific models we consider in sections 3 and 4.
2.1 Performance and sensitivity
Suppose for , system (1) has an asymptotically stable limit cycle solution with period . Let represent the rate at which the system advances relative to the outside world, and let represent the total progress achieved over one limit cycle, i.e.,
We note that in general, the instantaneous performance measure may depend both on the system variables and on the control parameter .
We consider the task performance of the system, denoted by , to be the progress divided by the limit-cycle period, or equivalently, the mean value of the rate of progress averaged around the limit cycle, defined as follows
(2) |
For a powerstroke-recovery system, we choose the time coordinate so that coincides with the beginning of the powerstroke phase, and write for the duration of the power stroke. We adopt the convention that during the recovery phase, the position with respect to the outside world is held fixed ()111This convention is consistent with a simplified single-limb swing-stance model (Markin et al.,, 2010; Spardy et al.,, 2011) as well as with a simplified model of feeding biomechanics in Aplysia californica (Shaw et al.,, 2015; Lyttle et al.,, 2017; Wang et al.,, 2022), and denote as the recovery phase duration (). In such a system, we write the performance (2) as
(3) |
Assume that , which represents the unperturbed load. When the system is subjected to a small static perturbation on the load, , the original limit cycle trajectory is shifted to a new trajectory , and its ability to resist the external change to maintain the performance is considered as a measure of robustness for the system. Since the sensory feedback pathways regulate the system dynamics, it would be desirable to obtain a feedback function so that the system is most robust against the load change, i.e.,
(4) |
Suppose we can expand around :
Then the minimization problem (4) to the first order reduces to
We quantify the sensitivity of the original system to be , which describes the (infinitesimal) response of the task performance to the external perturbation on the load. When with a certain feedback function, , which implies a strong ability of the system to maintain performance homeostasis222Homeostasis refers to the situation when a quantity remains approximately constant as a parameter varies over some range (Golubitsky and Stewart,, 2017, 2018; Golubitsky and Wang,, 2020). For limit cycle systems, Yu and Thomas, (2022) defined a homeostasis criterion in terms of the zero derivative(s) of the averaged quantity with respect to the perturbation parameter.. Define a (linear) functional , and the problem falls into a functional minimization problem, , which attains the minimum with some feedback function when the functional derivative . Finding when the underlying system has a limit cycle (as opposed to the more often studied case of fixed-point homeostasis) is an open problem that we do not attempt to solve here. Instead, we focus on a limited range of functions with practical significance and investigate the constrained optimization problem for the optimal feedback structure within that range.
2.2 Variational analysis
Variational analysis using the infinitesimal shape response curve (iSRC) and the local timing response curve (lTRC) is the key tool in our derivation of the performance sensitivity . In this section, we present a brief review of the theory and then provide two analytical methods to calculate the sensitivity. More mathematical details are given in Appendix C.
Consider system (1) written in the form
(5) |
where and is the corresonding vector field parameterized by the load . For convenience, we write trajectories as and as ; we write other quantities similarly. Expanding the perturbed trajectory yields
(6) |
The timing function rescales the perturbed trajectory to match the unperturbed trajectory, so that the series (6) is uniform with respect to the time coordinate. The linear shift in the shape of the unperturbed trajectory, , is referred to as the iSRC. It satisfies a nonhomogeneous variational equation (Wang et al.,, 2021; Yu and Thomas,, 2022)
(7) |
where measures the local timing sensitivity to the perturbation. The initial condition for the iSRC equation (7) is
(8) |
where and represent the intersection points of the trajectories with an arbitrary Poincaré section transverse to both the perturbed and unperturbed limit cycles, so that indicates the linear displacement of the unperturbed intersection point. For more details about the iSRC, see Appendix C.2, Wang et al., (2021), and Yu and Thomas, (2022).
To solve the iSRC equation (7), the lTRC is built to yield the timing sensitivity local to each phase of the motion. Approximate the perturbed phase durations by
Wang et al., (2021) and Yu et al., (2023) developed a formula to calculate the first-order approximation for the duration change in phase , given by
(9) |
Here, and denote the unperturbed entry point to, and exit point from, the specific phase, respectively. Vector , defined to be the gradient of the remaining time of the trajectory until exiting phase , is referred to as the lTRC for phase . It satisfies the adjoint equation
with a boundary condition
where is a normal vector of the exit boundary surface at . When the vector field changes discontinuously across the surface defining the boundary between two regions, the Jacobian should be evaluated as a one-sided limit, taken from the interior of the local region. With a linear time scaling for (6) (and setting to be the start of the power stroke), i.e.,
the local timing sensitivity function in (7) reduces to
which can be obtained by using equation (9) for each phase. See Appendix C.1, Wang et al., (2021), and Yu et al., (2023) for more details about the lTRC formulation.
The variational analysis above allows us to analyze the performance sensitivity of powerstroke-recovery systems. In Yu and Thomas, (2022) we provided a formula for the sensitivity of any averaged quantity with respect to an arbitrary control parameter, as long as the quantity of interest does not have explicit dependence on the control parameter. We generalize the approach of Yu and Thomas, (2022) to allow for dependence of the instantaneous performance on both the state and the parameter , and obtain
(10) |
The first term in the integral arises from the impact of the perturbation on the shape of the trajectory () as well as directly on the quantity of interest (. Here denotes the proportion of the powerstroke duration within the period (). The second term indicates the impact of the perturbation on the timing of the trajectory, in that represents the linear shift in in response to the perturbation, which can be analytically evaluated by
Given the special structure of powerstroke-recovery systems, we can derive a more succinct expression for . For any value of the second definition in (2) gives
(11) |
Recall at , we write and as shorthand for and . We can expand the perturbed progress and period around as
where is approximately given by the net change of the mechanical component of the iSRC (cf. equations (7) and (8)) within the powerstroke phase, and is the linear shift in the total period, readily given by (cf. equation (9)). Therefore, equation (11) at becomes
(12) |
which, like (10), incorporates both the shape and timing effects of the perturbation in two distinct terms. Equation (12) also suggests that the sensitivity can be directly given by the absolute difference between the first-order timing change and shape change induced by the perturbation. When the two effects completely offset each other, the system achieves “perfect” robustness.
The two expressions given by (10) and (12) for calculating the sensitivity of the task performance for powerstroke-recovery systems allow us to compare different sensory feedback mechanisms in pursuit of an efficient and robust motor pattern.
In the following sections we develop two illustrative examples: an abstract CPG-motor model introduced in Yu and Thomas, (2021), and an unrelated realistic locomotor model studied in Markin et al., (2010) and Spardy et al., (2011).
The two examples show a variety of differences in their model construction, but our analytic framework is broad enough to address both and give useful insights.
Simulation codes required to produce each figure are available at
https://github.com/zhuojunyu-appliedmath/Powerstroke-recovery.
3 Application: HCO model with external load
In Yu and Thomas, (2021), we studied a simple closed-loop model combining neural dynamics and biomechanics, as sketched in Fig. 1. The CPG sytem comprises a half-center oscillator (HCO) with two conductance-based Morris–Lecar neurons (Morris and Lecar,, 1981; Skinner et al.,, 1994). Outputs from the HCO drive a simple biomechanical system, which follows a Hill-type kinetic model based on experimental data from the marine mollusk Aplysia californica (Yu et al.,, 1999). Sensory feedback from the periphery couples the body and brain dynamics, allowing the system to interact with the changing outside world, and to modulate the central neural activities adaptively. However, the previous study of this model did not explore the performance with respect to a physical task; rather, the CPG followed an autonomous clocklike pattern. To perform a more meaningful analysis for understanding principles of closed-loop motor control, here we augment the model from Yu and Thomas, (2021) by incorporating a mechanical load exerted in a specific direction with recurrent engagement and disengagement with the system, which enables us to apply our quantitative measures of progress and sensitivity of the system.
3.1 The equations of the HCO model
The model equations we consider are as follows. For and ,
(13) |
Variable denotes the membrane voltage for HCO neuron cell , and is the gating variable for the potassium current in cell . The two neuron cells are coupled by fast inhibitory synapses, and the coupling function is given by
which closely approximates a Heaviside step function with denoting the synaptic threshold.
In the third equation of (13), represents the activation of the th muscle. The neural outputs from the HCO drive the associated muscle, modeled as
The biomechanics is represented by the movement of an object (nominally, a pendulum or limb), with each side connected to one of the two muscles. The object position relative to the center of mass of the organism, denoted as , is controlled by the muscle forces and acting on it. An external load is exerted on the object only during the powerstroke phase, as specified by the indicator variable defined to be
We assume that the powerstroke phase is at work when is in the active state (), whereas the load is absent from the system when is inhibited (). Parameter describes the strength of the load, which is considered as the perturbation parameter for this model.
The system completes an intact closed loop through the sensory feedback induced by the biomechanics on the CPG in the form of feedback currents, e.g.,
where is the length of muscle , and the function describes the feedback synaptic activation. We assume that the feedback conductance has fast dynamics, following a sigmoid function. As discussed in Yu and Thomas, (2021), the feedback synaptic architecture affords eight variations, depending on whether the feedback is (i) inhibitory or excitatory, (ii) activated by muscle contraction or muscle stretch, and (iii) modulatory on the contralateral or ipsilateral neuron. For example, when the feedback current is inhibitory to its contralateral neuron and activated when the muscle is contracted, then for the -equation we set mV and
(14) |
Figure 1 illustrates the system controlled by the inhibitory-contralateral-decreasing feedback mechanism, and Fig. 2 shows a typical solution for the system. In contrast, setting mV for the excitatory feedback current, or setting the sigmoid function to be increasing for the muscle-stretch activated case, or changing -dependence to -dependence for the ipsilateral mechanism, would specify other possible mechanisms. We compare the performance and sensitivity of all different realizations of each of the eight variations of the feedback control scheme below. The force terms , as well as additional details about the functions in system (13), parameter values used for simulations, and simulation codes are given in Appendix A.
3.2 Analysis of the HCO model
In order to establish measures of performance and sensitivity, we assume that the system advances only during the powerstroke phase. That is, the rate at which the system makes progress is given by
as indicated in Fig. 2D, F. Therefore, by (3), the performance (i.e., the average rate of progress) is
(15) |
When a sustained small perturbation is applied to the load, , the solution trajectory shifts in both its shape and timing, and the performance of the perturbed system is consequently different from the performance of the unperturbed system. As a reference, Fig. 3A, B, E, F compare the trajectories of the unperturbed solution with for the system shown in Fig. 2 and the perturbed solution with , and Fig. 3C, D, G, H illustrate the iSRC of the unperturbed trajectory, specified by Poincaré section . Note that for visual convenience, the large perturbation () is applied here, but in our actual analysis the perturbation magnitude should be small (). Our analysis yields several observations about the role of the inhibition-contralateral-decreasing sensory feedback in regulating the system’s response to the perturbation, as discussed in detail below.
With the perturbation (increased load), the transition from the power stroke to recovery is advanced, i.e., the powerstroke phase is shorter. As indicated in Fig. 3F, the immediate effect of the perturbation is the positive displacement and slower change rate in the -variable, which occurs because the object is being pulled by the stronger load opposite to its movement direction. Correspondingly, muscle 2, whose length is , is more contracted than it is in the unperturbed case. The sensory feedback current injected to neuron 1, with synaptic activation given by (14), is therefore larger and gives more inhibition to the active neuron 1. As a result, the active crosses the synaptic threshold and terminates the powerstroke phase at an earlier time, as indicated in Fig. 3A. The iSRC in each direction is consistent with the associated trajectory comparison. The significant negative peak in the -component at the end of powerstroke phase (panel C) suggests that of the perturbed solution at the rescaled time already decreases to the synaptic threshold and jumps down to the inhibited state, indicating the transition out of the power stroke is advanced. Since directly impacts , we see a different effect on (panel H) than the other variables.
Both efficiency (high performance) and robustness (low sensitivity) are important features of motor control systems interacting with the outside world. The performance for the perturbed system in Fig. 3 is smaller than the unperturbed system (, ), due to the larger magnitude in the progress decrease relative to that in the period. To measure the ability of maintaining the performance, we quantify the sensitivity of the original system in response to an infinitesimal sustained perturbation, following (10), to be
(16) |
One can also estimate the sensitivity by (12), where corresponds to the -component of . That is,
To evaluate the joint goals of high performance and low sensitivity, and to investigate how they are affected by sensory feedback, we will simultaneously study the two measures plotted together, while manipulating the shape of the feedback activation function. Specifically, we will vary the steepness parameter and position/half-threshold parameter of the sigmoid synaptic feedback activation function . Figure 4 shows the results for all eight sensory feedback mechanisms.
The eight superficially distinct feedback architectures can be reduced to four fundamentally different mechanisms in terms of their performance and sensitivity. The contralateral mechanism of muscle-stretch activated (increasing) current with threshold (), is equivalent to the ipsilateral mechanism of muscle-contraction activated (decreasing) current with threshold . Specifically, substituting the contralateral-increasing feedback activation with to the -equation of system (13) yields
where the last equation is exactly the case for the ipsilateral-decreasing feedback with . Note that for the unloaded model in Yu and Thomas, (2021), only the inhibition-excitation property of the feedback makes a fundamental difference regarding the stability and robustness of the system. However, when we incorporate mechanical interactions with an external substrate, the activating property of the feedback must be taken into account. In the following, we will discuss the performance and sensitivity for the four contralateral feedback mechanisms, which we call inhibition-increasing (II), inhibition-decreasing (ID), excitation-increasing (EI), and excitation-decreasing (ED). The other four ipsilateral feedback mechanisms can be applied accordingly.
3.3 Performance and sensitivity of the HCO model
Our analysis of the HCO model, subject to an applied external load, leads to the following observations from Fig. 4:
-
1.
Excitatory feedback is advantageous over inhibitory feedback in terms of performance.
-
2.
The qualitative patterns of performance are reversed with respect to the steepness of synaptic activation function in the activating and inactivating feedback mechanisms.
-
3.
The qualitative patterns of sensitivity are reversed with respect to the steepness of synaptic activation function in the excitatory and inhibitory feedback mechanisms.
-
4.
When the sigmoid activation function is approximately linear over the working range of the limb, the performance-sensitivity changes approximately linearly with the slope of .
-
5.
As the working range of the limb extends beyond the linear regime of , the performance-sensitivity curve can become strongly nonlinear and even non-monotonic, leading to well-defined simultaneous optima in both performance and sensitivity.
We discuss each of these points in turn below.
Excitatory sensory feedback outperforms inhibitory feedback.
This conclusion is evident from Fig. 4, in the higher location of the performance-sensitivity pattern for each of the excitation mechanisms relative to that for all of the inhibition mechanisms. To understand the advantage of excitatory over inhibitory feedback in this model system, Fig. 5 compares the trajectories of two systems, one with excitatory constant feedback and the other with inhibitory constant feedback. These systems correspond respectively to the two green dots in Fig. 4A. In the excitatory system, the extra excitation to due to the feedback has two opposing effects. On the one hand, it advances the time at which the active neuron “jumps down” to the inhibited state and thus shortens both the powerstorke phase and the total period , relative to the system with constant inhibitory feedback. Panel F of Fig. 5 shows the projection of the trajectory on the plane with points plotted at constant time intervals. Note the rapid change in the voltage component relative to the slower gating variable. The significant timing change results from the exponential deceleration of the dynamics of the active neuron when approaching the jump-down point, as shown by the contraction of points before jumping in panel F. Although the projection of the two trajectories during the active state does not differ much spatially, the difference in time needed to cover the small spatial difference is significant. On the other hand, the constant excitatory drive also reduces the net progress of the system per cycle. In particular, the progress of the system declines due to the shorter movement time (panel E). The net performance (, equation (2)) is determined by both the timing and shape of the trajectory. For the excitatory system, the resulting performance is in fact larger than for the inhibitory system, because the relative change in period is larger than the relative decrease in progress. Therefore, regardless of the structure of the feedback pathway, any system equipped with the excitatory sensory feedback is always advantageous over the system with the inhibitory sensory feedback in terms of their performance.
The performance patterns in the decreasing and increasing feedback mechanisms are qualitatively reversed with respect to the steepness of .
For example, as changes from constant to approaching a Heaviside step function, the performance of the II mechanism with (red dots in Fig. 4A bottom) monotonically increases while the performance of the corresponding ID mechanism (red + in Fig. 4A bottom) decreases monotonically. Figure 6 illustrates two systems controlled by the ID mechanism with but different values. When is more shallow (dashed), the inhibitory feedback current to cell 1 is less intense, due to the smaller synaptic activation over the contraction regime of muscle 2 (panel E). This makes neuron 1 (the neuron driving the powerstroke) more active, switching earlier to the inhibited state, and thus giving a shorter duration power stroke (panel A). The progress of the limb over the shorter powerstroke phase is however larger, which is opposite to the excitatory situation in Fig. 5. This outcome occurs because the effect of the duration decrease is not comparable to the effect of the faster velocity around the end of power stroke. The progress velocity over the powerstroke phase is controlled by the force generated by muscle 1, which becomes stronger due to the increased muscle activation (panel C). Therefore, by decreasing the steepness of for the ID mechanism, the muscles act on the limb in a stronger and faster fashion, leading to the enhanced performance. In contrast, by applying a similar analysis, we observe that the corresponding II mechanism, for which the muscle-stretch activated feedback current gives more inhibition to the active neuron 1, has the opposite effect (not shown). This example illustrates how two mechanisms with the opposite activating property of the sensory feedback can have qualitatively distinct performance changes when the sensitivity of the feedback pathway to the biomechanics is varied.
The sensitivity profiles of the excitatory and inhibitory systems to an applied load are qualitatively reversed with respect to the steepness of the synaptic activation function.
Consider the two contralateral-decreasing mechanisms with (see the two blue dot curves in Fig. 4A). As becomes smaller, the sensitivity of the ED mechanism first decreases to almost 0 (perfect robustness) and then increases, whereas the sensitivity of the corresponding ID mechanism first increases to the maximum and then decreases. Recall that following equation (12), the sensitivity can be written as
This expression contains two terms — accounts for the effect of the perturbation on the system progress (shape), while accounts for the effect on the period (timing). To the extent that the two effects are large and of opposite signs, the system becomes more sensitive; in contrast, the cancellation of the two effects leads to a robust system. Figure 7 shows the two effects in the ED mechanism (red) and ID mechanism (blue) with and varied over . We observe that the stronger mechanical load exerts positive effects on both the timing and shape for the system with the excitatory sensory feedback mechanism, by prolonging the limit cycle period and increasing the progress of the limb (positive red curves), but the shape effect is more profound. As decreases, the improvement in the progress becomes less significant, while the timing changes more dramatically, which leads to a reduction in the sensitivity. The sensitivity minimum is attained when the two effects completely offset each other at , where the system is most robust against the perturbation. In contrast, the perturbation allows the system with the inhibitory feedback to make larger progress in a shorter time (negative blue solid and positive blue dashed); this possibility was mentioned in the discussion of Fig. 6. Although making greater progress in a shorter time in response to perturbation improves performance, it is not beneficial to the robustness of the system in terms of maintaining performance homeostasis. Moreover, as decreases, both effects become stronger, inducing the system to be even more sensitive. Apart from this example, the excitation-inhibition property of the sensory feedback in this model always offers qualitatively reversed sensitivity patterns when we consider the steepness variation of the feedback activation curve. The variational analysis serves as a tool to examine the coordinated effects of perturbation on the trajectory geometry and timing, and to identify mechanisms with superior robustness.
When is approximately linear over the working range of the muscles and limb, the performance-sensitivity curve changes approximately linearly with .
As an illustration, Figure 8 zooms in on the patterns shown in Fig. 4A with varied among , over which is approximately linear within the possible range of muscle length. On this scale, a fundamental ambiguity arises on account of the dual goals of performance and robustness. Without specifying a relative weighting between these two quantities, there is no well justified way to choose which of the three traces in the upper ensemble, all moving up and to the left, are preferred. The red, black, and blue traces all simultaneously increase performance while decreasing sensitivity, relative to the constant feedback case (green dot). Among the lower ensemble, the blue curve can be rejected as its components exhibit a tradeoff: the robustness is enhanced with decreased performance. Moreover, within the linear regime, the performance and sensitivity can both be improved indefinitely by increasing the feedback gain. The possibility of disambiguating the two effects and finding a globally optimal solution arises only when the parameters are varied enough for nonlinear effects to come into play, as shown next.
The nonlinearity of the synaptic feedback activation function allows the existence of optimal combinations of performance and sensitivity.
As the linear regime of extends beyond the working range of muscle length and the curvature of the sigmoid becomes more pronounced, some sharp turning points arise in the performance-sensitivity curve of each feedback mechanism, either in terms of performance or sensitivity or both. In general, along a continuous curve in the (sensitivity , performance ) plane, indexed by a parameter (here, the sigmoid steepness parameter ), there will be an optimal region marked at one end by the condition and at the other end by . Depending on the relative weight given to and , the optimal value of will place the system somewhere within this segment of the curve. In the case of the HCO system, the optimal curve among those tested (, top blue dots in Fig. 4A) makes a hairpin turn in the upper left corner of the plot (blue arrow), leading to relatively unambiguous identification of the optimal value of the slope parameter . In such cases, these optimal points, or narrowly identified optimal regions, indicate the possibility of well-defined simultaneous optima in both of the performance and sensitivity patterns. Thus our results suggest the possibility to realize these joint goals through adjusting the structure of the sensory feedback mechanisms more broadly.
The preceding observations have provided several insights that could support the design of sensory feedback pathways, such as the selection of excitatory versus inhibitory feedback currents, and activation versus inactivation with muscle contraction, as well as the shape of the feedback activation curve, to promote efficiency and robustness in other more realistic HCO-motor models with analogous configurations.
4 Application: Markin hindlimb model with imposed slope
The second closed-loop powerstroke-recovery example we consider is based on a neuromechanical model proposed by Markin et al., (2010), as sketched in Fig. 9. The model consists of a spinal central pattern generator controlling the movement of a single-joint limb. The CPG sends output via efferent activation of two antagonist (flexor and extensor) muscles to a mechanical limb segment, which in turn generates afferent feedback signals to the CPG. The model system performs rhythmic locomotion, comprising a stance phase, during which the limb is in contact with the ground, and a swing phase, during which the limb moves without ground contact. Thus the system falls within the class of powerstroke-recovery systems.
Everyday experience teaches us that normal walking movements are robust against gradual changes in the slope of the terrain, although the detailed shape and timing of the limb trajectory changes as one ascends (or descends) a steeper or shallower incline. As a proxy for this form of parametric perturbation of rhythmic walking movements, we extend Markin et al’s model to include an incline parameter, simulating the effects of a change in the ground slope. We use this parameter in our analysis of performance and sensitivity, to define the sensitivity of the system’s progress relative to the external substrate (the ground) in response to this environmental change.
4.1 The equations of the Markin model
Following Markin et al., (2010), we model the central pattern generator as a multi-layer circuit, consisting of a half-center rhythm generator (RG) containing flexor neurons (RG-F) and extensor neurons (RG-E). These RG neurons project to pattern formation (PF) neurons (PF-F and PF-E, respectively) and to inhibitory interneurons (In-F and In-E) which mediate reciprocal inhibition between the flexor and extensor half-centers. Given sufficient tonic supra-spinal drive, the CPG generates rhythms of alternating activation of flexor and extensor neurons, and the output of PF neurons induces alternating activity in flexor and extensor motor neurons (Mn-F and Mn-E). An additional circuit of interneurons (Int and Inab-E) provides a disynaptic pathway from the extensor side to Mn-E.
The dynamics of the RG, PF, and Mn neurons are each described by two first-order ordinary differential equations, governing each cell’s membrane potential and the slow inactivation gate of a persistent sodium current:
(17) |
Here, , , and refer to the persistent sodium current, potassium current, and leak current, respectively, described by
Excitatory and inhibitory currents to neuron are respectively represented by and , given by
(18) |
The nonlinear function describes the output activity of neuron , defined to be
(19) |
where is the synaptic threshold. Parameter defines the weight of the excitatory synaptic input from neuron to neuron , while defines the weight of the inhibitory input from to ; represents the weight of the excitatory drive to neuron ; defines the synaptic weight of afferent feedback to neuron , with the feedback strength . Details on the feedback terms will be provided at the end of this section.
The dynamics of the interneurons, In-F, In-E, Int, and Inab-E, are each described by a single first-order equation:
(20) |
where the currents are in the same form as above. As shown in Fig. 9, the source of excitatory inputs to these interneurons comes from RG, supra-spinal drive, and sensory feedback. Note that In-E and In-F in particular do not receive any inhibitory input or excitatory supra-spinal drive, so the right-hand side of their voltage equation has and .
The motor neurons, Mn-F and Mn-E, respectively activate two antagonistic muscles, the flexor (F) and extensor (E), controlling a simple single-joint limb. The limb motion is described by a second-order differential equation:
(21) |
where represents the angle of the limb with respect to the horizontal. The first term accounts for the moment of the gravitational force; and are the moments of the muscle forces; denotes the moment of the ground reaction force which is active only during the stance phase, given by
(22) |
Note that the stance (powerstroke) phase is defined when the limb angular velocity is nonnegative, while the swing (recovery) phase occurs when the velocity is negative. We expand on the original model from Markin et al., (2010) by introducing parameter to describe the slope of the ground whereon the limb stands. Thus will play the role of the load parameter subjected to perturbations for this model.
The feedback signals from the extensor and flexor muscle afferents provide excitatory inputs to the RG, PF, In, and Inab-E neurons. Muscle afferents provide both length-dependent feedback (type Ia from both muscles and type II from the flexor) and force-dependent feedback (type Ib from the extensor). Linear combinations of feedback terms , written as in (18), are fed into each side of the model — Ia-F and II-F feedback go to the flexor neurons and Ia-E and Ib-E to the extensor neurons. The feedback terms are in the form
(23) |
For more details about the mathematical formulations, functional forms, and parameter values of the entire model, see Appendix B.
4.2 Analysis of the Markin model
Figure 10 shows the time courses of the output of neurons , limb angle , and feedback activity in the default flat-ground system where . Unlike the HCO model, where the active states of neurons overlap with the powerstroke and recovery phases, here the extensor and flexor active phases are slightly shifted relative to the stance and swing phases Spardy et al., (2011). The excitatory feedback of types Ia and II increase during the silent phase of the associated neuron receiving the signal, reaching a peak just before the target neuron becomes active; these feedback signals then decreases during the active phase of the neuron. In contrast, the Ib-E feedback signal, which is solely dependent on the extensor muscle force, is active only when the extensor neurons are active. It remains low until the onset of activation in the extensor units, at which point it jumps up to a high level and then recedes. The different types of the feedback pathways induce distinct performance-sensitivity patterns, as we will discuss in §4.3.
For this model, we define the performance to be the average distance the limb moves along the ground during the stance phase, given by
(24) |
Here denotes the stance duration and is the limb length. Figure 11 compares the time series of the perturbed system subjected to a small change in the ground slope () with the unperturbed system. As indicated by the green dotted vertical line, the perturbation prolongs the stance phase and delays the transition to the swing phase, in that the system decelerates because of the steeper ground slope (panel B). Consequently, the progress becomes smaller (panel D), and the performance is further deteriorated due to a longer period (, ). Note that this model and the HCO model have different conditions for transitioning between the powerstroke and recovery phases. In the HCO model the transition is determined by the neuron activation, while in the biophysically-grounded Markin model the transition is determined by the mechanical condition . This difference accounts for the opposite responses of the transition timing to analogous environmental challenges: when increasing the difficulty of the task the HCO powerstroke-to-recovery transition moves earlier, while the Markin model’s transition is delayed. Put another way, the powerstroke of the perturbed HCO system contracts relative to the net period, while in the Markin model it expands (cf. Fig. 3).
In the following, we study how the performance-sensitivity pattern is affected by the strength of each afferent feedback pathway, represented by where in (18). Figure 12 shows the patterns as one of the four feedback strengths is varied, with the other three strengths fixed at the normal strength . Before we discuss the patterns in §4.3, note that the possible range of each strength allowing for stable progressive locomotor oscillations is remarkably different:
In particular, the system can maintain stable oscillations upon cutting off the II-F or Ib-E feedback pathway, but cannot sustain oscillations without the Ia-type feedback. We finish this subsection by examining the mechanism underlying the failure of movement as the strength parameter is out of the range.
The Markin model system features fast-slow dynamics (Rubin and Terman,, 2002), as the persistent sodium inactivation time constant in (17) is large over the relevant voltage range, so evolves on a slower timescale than . The activity of neuron can be therefore determined from the location of the intersection of its nullclines in the phase plane. To illustrate the mechanism underlying the transition between phases, Fig. 13 shows the nullcline configuration and the corresponding positions of RG and In neurons in the default system around the CPG transition. Starting from the extensor-inhibited state (panels A), because the -nullcline and -nullcline of RG-E intersect at a silent stable fixed point, RG-E cannot escape from the silent state and trade dominance on its own. However, In-E, due to the increasing excitatory inputs from Ia-E and Ib-E feedback (Fig. 10), is able to cross the synaptic threshold first and begin to inhibit RG-F (panels B). This inhibition raises the -nullcline of RG-F, followed by the decrease of RG-F voltage, weakening its excitation to the downstream In-F. This reduction of excitation results in the decrease of In-F voltage. Consequently, RG-E receives less inhibition from In-F, lowering the -nullclines of RG-E, such that the critical point moves to the middle branch of its -nullcline, and hence, RG-E can reach the left knee and jump to the right branch. This transfer of active and silent states indicates that the In cells and excitatory sensory feedback dominate the transitions in the CPG. In the case without sufficiently strong feedback inputs to the silent In cell, the ipsilateral RG neuron will become deadlocked in silence, and thus the locomotion will fail.
Disentangling the effects of changing the strength of an afferent feedback pathway is nontrivial, because of the multiplicity of pathways impacting the activities of many neurons and limb within the system. For instance, Fig. 14 shows the feedback and neuron dynamics in the cases (dashed) and (solid), which is approaching the minimal allowable strength of Ia-F. Varying the strength of the Ia-F pathway alone leads to a chain of changes in all system components including the extensor feedback, and significantly impacts the transition timing. Although the magnitude of the feedback signals does not differ much in the two systems, yet the system with smaller accumulates feedback excitation more slowly. This circumstance is due to the interplay of spatial and timing measures of distance in the fast-slow dynamical systems (Terman et al.,, 1998; Rubin and Terman,, 2002; Yu et al.,, 2023), in which a small spatial difference translates into a significant temporal extension. Decreasing the feedback strength further (below 0.59) makes it impossible for the extensor feedback to jump up; consequently the In-E neuron will not be facilitated to escape from the inhibited state (not shown) and the movement will stall. This outcome may seem counterintuitive, since reducing the flexor feedback strength should make the flexor In more difficult to escape; however, its direct effect on the mechanics in turn influences the extensor feedback to an analogous (Ia-E) or even more significant (Ib-E) extent (see equations (23)). We observe similar situations beyond the nonzero extrema of other feedback strengths, in which either In-F or In-E fails to escape from the silent state and compromises the whole locomotor oscillation.
4.3 Performance and sensitivity of the Markin model
When analyzing the abstract HCO model, we considered both ipsilateral and contralateral feedback, both excitatory and inhibitory. In contrast, each of the three sensory feedback pathways in the Markin model is ipsilateral and excitatory. Although the two models are quite different in their level of details, feedback variations, and the type of perturbation, we compare and contrast the two models as far as we can in the Discussion. In this section, by varying the strength of each feedback pathway, we study the performance-sensitivity patterns of the Markin model as shown in Fig. 12 and yield the following observations.
Trade-offs between performance and sensitivity often occur as the afferent feedback strength changes.
As the strength parameter is varied, the system shows a performance-sensitivity tradeoff when the curve either moves up and to right (performance improves while robustness decreases) or moves down and to left (robustness increases while performance declines). Specifically, the system is not able to generate simultaneously efficient and robust movement with the whole flexor-feedback (Ia-F and II-F) variations as well as with the large Ia-E and small Ib-E variations. For example, Fig. 15 compares the performance of the default system to the system with a stronger Ia-F pathway (). Similarly to Fig. 14, the increase in the Ia-F strength advances the stance-swing transition, leading to a shortened working period and reduced progress. Since the change in the shape is more significant than the change in the timing, the performance () declines. Note that in the HCO systems featuring either excitatory feedback or inhibitory feedback, the dominance between the shape and timing is reversed (cf. Fig. 5). These examples indicate the interplay of temporal and spatial effects in the performance measure. Likewise, the sensitivity measure, given by (10) or (12), also consists of two factors — the timing response to the perturbation and the shape response. Although the system with strengthened Ia-F feedback is inferior to the default system in terms of performance, it benefits from enhanced robustness. In contrast, weakening the Ia-F strength contributes to an improved performance at the cost of responding more sensitively to external perturbations. With such a performance-sensitivity tradeoff, the system cannot simultaneously achieve the dual goals of efficiency and robustness, especially when varying the flexor-feedback pathway(s) alone.
Force-dependent feedback gain can optimize both performance and robustness simultaneously.
We observe that the force-dependent feedback pathway (type Ib from the extensor) impacts performance and robustness very differently than the length-dependent feedback pathways (types Ia and II). In contrast to the length-dependent feedback, the force-dependent feedback is only active during the active phase of the associated motor neuron, reflecting the fast-slow dynamics of the CPG (see equations (23) and Fig. 10). Notably, the performance-sensitivity tradeoff disappears with variations, and an optimal value of arises, giving the best combination of the pair (purple arrow in Fig. 12). Around the optimal point, the performance saturates while the sensitivity changes dramatically, and at , the sensitivity attains zero, indicating that the system achieves infinitesimal homeostasis (Yu and Thomas,, 2022) concurrent with peak performance. We verify this perfect robustness by plotting the coordinated effects of the slope perturbation on the system trajectory’s shape and timing in Fig. 16 as varies over . For values of above the optimum, the robustness decreases dramatically. Indeed, the underlying oscillation fails when reaches about 7.7, presumably due to a global bifurcation. As the parameter approaches this point, the trajectory period and amplitude might be insensitive to the parameter variation but the inherent sensitivity to perturbations may dramatically increase. This circumstance is also observed in the HCO model when (subcritical) pitchfork bifurcations occur (Yu and Thomas,, 2021). On the contrary, with a strong Ia-E pathway (left corner of Fig. 12), the system maintains robustness while its performance significantly decreases, which also appears around a (subcritical) Hopf bifurcation as shown in Yu and Thomas, (2021). Although a detailed bifurcation analysis of oscillation termination is beyond the scope of this paper, our observations demonstrate the difficulty of achieving the joint goals of high performance and low sensitivity in a realistic rhythmic physiological system. The variational method, by incorporating both temporal and spatial factors, offers the possibility of an analytic framework that may assist to identify the optimal sensory feedback control mechanism in regulating the system’s efficiency and robustness.
5 Discussion
Efficiency and robustness are important aspects of control systems in both engineering and biological contexts. The application of conceptual and mathematical ideas from control theory in the physiological sciences has a long history (Baylis,, 1966; Grodins et al.,, 1967; Khoo,, 2018). In motor control, viewed through the lens of systems engineering, the concepts of robustness and efficiency can help quantify different aspects of system performance underlying biological fitness. The roles of feedforward, reciprocal, and feedback motor pathways in contributing to robustness and efficiency must be investigated on a case-by-case basis. In many instances, control theorists have observed tradeoffs between efficiency and robustness, although the definitions of these terms may vary from one context to another (Kuo,, 2002; Boulet and Duan,, 2007; Ronsse et al.,, 2008; Vasconcelos et al.,, 2009; Alfaro et al.,, 2010; Sariyildiz and Ohnishi,, 2013). These authors describe a fundamental tradeoff between efficiency and robustness, in the sense that the system improves efficiency at the cost of becoming increasingly fragile to a slight change in any parameter. Hence, an important aspect of controller design is to achieve an acceptable balance between the goals of increasing performance versus maintaining robustness to external changes, when these objectives are in conflict. Throughout this paper, we interpret efficiency as quantifying how well a system performs, and robustness as quantifying the stability of efficiency with respect to parametric perturbations.
A remarkable feature of rhythmic motor systems underlying a wide range of animals’ physiological behaviors (crawling, breathing, scratching, swallowing) is the closed-loop control structure, which integrates brain, body, and environment via sensory feedback (Chiel and Beer,, 1997; Diekman et al.,, 2017; Lyttle et al.,, 2017; Chen et al.,, 2021; Korkmaz et al.,, 2021). In particular, feedback control plays a significant role in regulating the robustness and efficiency of these rhythmic movements.333The significance of sensory feedback may vary depending on the length scale of the animals, e.g. in small insects versus large mammals (Sutton et al.,, 2023). Through afferent feedback pathways, sensory information can communicate current demands, thereby allowing the system to update the corresponding movement on short time scales. In this way, sensory feedback may improve performance; however, robustness remains a critical issue. For example, Kuo, (2002) and Yu and Thomas, (2021) discussed the fundamental properties of the trade-off between efficiency and robustness for a combined feedback–feedforward model for rhythm generation. These papers demonstrated that a system with pure feedback control, analogous to a chain reflex system, can compensate for unexpected disturbances, but shows poor robustness to imperfect sensors, and that the best trade-off requires a proper design incorporating both feedback and feedforward pathways.
In the present paper, we propose a unified framework to investigate how the architecture of sensory feedback influences the efficiency and robustness of rhythmic systems. In the particular context of a general powerstroke-recovery system, where the system repeatedly engages and disengages with the outside world, we measure the efficiency of the system, i.e. the task performance, in terms of the physical progress of the system over the working period (or equivalently, the mean rate of progress). This measure of the task performance is adopted from Lyttle et al.’s study of the feeding behavior of the sea slug Aplysia californica, which considered mean seaweed consumption rate (cm/sec) (Lyttle et al.,, 2017; Wang et al.,, 2022). This measure can be generalized to any powerstroke-recovery systems in which the progressive activity is important. We probe the robustness of the system by studying its ability to maintain task performance despite external perturbations. Following Yu and Thomas, (2022), we generalize the analytic formula for this robustness (or sensitivity) measure by using the infinitesimal shape response curve (iSRC) and the local timing response curve (lTRC) developed in Wang et al., (2021). These tools provide a mathematically grounded numerical quantification of the system’s response to perturbations. By simultaneously studying the two objectives, we compare the performance-sensitivity patterns obtained in different feedback architectures, and identify optimal designs achieving the joint requirements of high performance and low sensitivity.
Here, we consider the control problem with two specific neuromechanical systems that produce powerstroke-recovery oscillations. In the paradigmatic half-center oscillator (HCO)-motor model with an external load, adapted from Yu and Thomas, (2021), we conduct a comprehensive analysis of different feedback mechanisms in the hopes that we uncover insights that may apply more broadly. We study both contralateral and ipsilateral feedback architectures, both excitatory and inhibitory feedback modes, and both activating and inactivating pathways, in terms of their performance-sensitivity patterns. First, we find that exchanging excitatory and inhibitory architectures reverses the effect of increasing feedback gain on the sensitivity. Second, we observe that the performance response is reversed depending on whether the feedback is activated or inactivated by muscle contractions. More interestingly, as the sigmoid activation of the feedback signal (as a function of limb position) becomes increasingly nonlinear over the working range of the limb, a well-defined simultaneous optimum in performance and sensitivity arises. These findings may guide the selection of modeling frameworks to capture experimental observations or design aims for more realistic brain-body-environment systems with analogous architectures in future work.
In the hindlimb locomotor model (Markin et al.,, 2010), based on experimental measurements from spindle primary afferents in the cat, we modify the model by imposing a ground slope, and then manipulate the gain parameters for each of four sensory feedback pathways. We find that the force-dependent (type Ib-E) feedback can simultaneously optimize both performance and robustness, while the length-dependent (types Ia-E, Ia-F, and II-F) feedback variations result in marked performance-versus-sensitivity trade-offs. Indeed, as discussed in Yu and Thomas, (2021), the situation that the performance is increased by sacrificing the robustness is also observed in the abstract HCO-motor model when the system approaches a subcritical pitchfork bifurcation; and also, the situation that the robustness is improved at the cost of losing efficiency occurs when the system approaches a subcritical Hopf bifurcation.
However, it is difficult to pursue further comparisons between the two models. In the HCO model, we consider only one feedback pathway, solely dependent on the muscle length. In contrast, the locomotor model possesses multiple feedback channels, each of which relies on several components including muscle length, muscle velocity, motor neuron output, or muscle force. As argued in Katz, (2023), although focusing on a specific model can advance understanding of the neural basis of behavior in species with similar developmental and physical constraints, doing so may overlook insights that could be obtained by considering convergent evolution as a framework leading to more general principles. Our HCO model can provide a mechanistic understanding leading to predictions about structure-function relationships, and allows us to compare different kinds of architectures that are “equivalent” only at a more general level. On the other hand the detailed Markin model is just one specific architecture where perhaps alternative architectures could have evolved to serve the same purpose. Thus, although the analogy between the HCO-motor and the Markin-motor model is limited, both models are worth studying from a conceptual point of view. In particular, we note that the Ia- and II-type afferent feedback in Markin’s more realistic locomotor model is ipsilateral, excitatory, and increasing with muscle length; the performance-sensitivity pattern for this combination of elements nominally falls into the upper ensemble of the HCO model’s ipsilateral patterns (Fig. 4B, upper star traces). In the HCO setting, these patterns are always advantageous over other patterns in terms of performance, and give rise to the only (performance, sensitivity) optimum. This superiority could be interpreted as suggesting that our modeling analysis might be consistent with the action of natural evolution by which species adapt to environmental demands.
Apart from the restricted scope of feedback functions considered in our examples, it would be desirable to have greater insights into the way any form of sensory feedback shapes the attributes of motor control systems, both in terms of performance and in terms of robustness. However, given the complexity and difficulty to ascertain sensory feedback signals experimentally, it is natural to restrict attention to monotonically increasing or decreasing sensory feedback signals, such as an excitatory or inhibitory conductance in a central pattern generator that increases or decreases steadily as a function of limb position, muscle stretch, muscle velocity, or other biomechanical variables. A realistic example is the responses of cat spindle (group Ia and group II) afferents, which monotonically change with muscle length and velocity signals, and tendon organ (group Ib) afferents, which scales in approximate proportion to muscle forces (Prochazka and Gorassini,, 1998; Prochazka,, 1999). Another example is the chemosensory pathway in the respiratory control system, which has been modeled by a sigmoidal relationship between the arterial partial pressure of oxygen and the conductance representing external drive to the CPG (Diekman et al.,, 2017). Moreover, biological signals are typically limited in range. For instance, the firing rate of a neuron is bounded below by zero and typically bounded above by a maximal firing rate which is related, for example, to the neuron’s absolute refractory period.
The muscle dynamics of the HCO model was inspired in part by Aplysia’s I2 muscle (Yu et al.,, 1999). A nominal feeding model, featuring grasper-retraction (powerstroke) and grasper-protraction (recovery), has been well studied by Chiel and colleagues (Shaw et al.,, 2015; Lyttle et al.,, 2017; Wang et al.,, 2022). In their model, the proprioceptive feedback input is assumed to be affine linear with the grasper position — decreasing for the protraction neural pool while increasing for the retraction neural pool. Depending on the grasper position and proprioceptive neutral position, the feedback signal can switch between excitation and inhibition during a single swallowing cycle. In our HCO model, we only consider the same form of feedback functions for all neurons and the excitation-inhibition property does not change during the whole movement. These differences point out an interesting and important direction for further study — whether mixed inhibitory and excitatory feedback in a single pathway, or heterogeneous activating/inactivating feedback across different pathways — could provide more realistic models for a broader range of specific motor systems.
Finally, we apply a fast-slow decomposition analysis to the Markin model when discussing the failure of oscillations as the feedback gain parameter reaches the collapse point. Although the stance and swing phases in the locomotion system do not differ much concerning their time scale, it is nevertheless possible that in some powerstroke-recovery systems, the powerstroke, under a heavy load, could proceed on a slow timescale while the recovery could proceed on a fast timescale. In this case one might connect powerstroke-recovery systems to two-stroke relaxation oscillators as studied in depth by Jelbart and Wechselberger, (2020), an interesting direction for future work. Another avenue for future inquiry would be to consider the scaling of robustness and sensitivity structures, as studied here, in motor systems spanning a range of length and time scales (Sutton et al.,, 2023). The approaches we use here should be applicable to these systems as well.
6 Acknowledgments
This work was supported in part by (1) National Institutes of Health BRAIN Initiative grant RF1 NS118606-01; (2) the National Science Foundation under Grant No. DMS-1929284 while the authors were in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the program; and (3) the Oberlin College Department of Mathematics. The authors thank Hillel J. Chiel and Yangyang Wang for helpful discussions concerning the manuscript.
Appendix A HCO model details
For completeness, we list here the full equations of the HCO model with an externally applied load, as introduced in Section 3 and originally proposed in Yu and Thomas, (2021). For cell and ,
where
The force of muscle is given by
Here,
and
where
Note that in the -equation, we leave implicit a conversion factor from mV to Hz. The limb position is controlled by
where
Table 1 specifies the parameter values used for simulations. The simulation codes are available from https://github.com/zhuojunyu-appliedmath/Powerstroke-recovery. Instructions for reproducing each figure and table in the paper are provided (see the README file at the github site).
Parameter | Value | Unit | Parameter | Value | Unit |
---|---|---|---|---|---|
1 | F/cm2 | 15 | mV | ||
0.8 | A/cm2 | 15 | mV | ||
0.005 | S/cm2 | 2 | mV | ||
0.015 | S/cm2 | 0.0005 | msec-1 | ||
0.02 | S/cm2 | as in Figure | cm | ||
0.005 | S/cm2 | as in Figure | cm | ||
0.001 | S/cm2 | 2 | none | ||
-50 | mV | 2 | N | ||
100 | mV | 10 | N | ||
-80 | mV | 4103 | Nmsec/cm | ||
-80 | mV | 2.45 | none | ||
mV | 2 | none | |||
0 | mV | 0.165 | none | ||
15 | mV | 0.703 | none | ||
0 | mV |
Appendix B Markin model details
Based on Markin et al., (2010) and Spardy et al., (2011), the following provides modeling details not specified in Section 4, especially on the biomechanics and feedback dynamics. The limb motion is described by a second-order differential equation:
where is the moment of inertia of the limb with respect to the suspension point and is the angular viscosity in the hinge joint. The first term accounts for the moment of the gravitational force; and are the moments of the muscle forces (see below); represents the moment of the ground reaction force which is active only during the stance phase, given by
The parameter describes the slope where the limb stands on, which is considered as the load parameter subjected to perturbations for this model. This parameter is not included in the original model Markin et al., (2010); in effect in the original model .
In the muscle model, the muscle length is calculated as
and the moment arm is given by
for the flexor msucle and extensor muscle, respectively. Muscle velocities are defined as
where denotes the limb angular velocity. The total force in each muscle follows the Hill-type model:
Constant is the maximal isometric force; is the output of the corresponding motorneuron. The force–length dependence is given by444Careful investigation on the code provided by Markin revealed a typo in the function in the original papers. Here is the correct version.
where is the normalized muscle length corresponding to , i.e. . The force–velocity dependence is given by
where . The passive force is calculated as follows:
Given the above, the muscle moments are given by and .
To complete the model, the afferent feedback terms fed into the CPG neurons are in the form
Ia | |||
II | |||
Ib |
Here, denotes the normalized muscle velocity (); denotes the normalized muscle length ( if and 0 otherwise); denotes the normalized muscle force ( if and 0 otherwise). The other functions in the CPG model are
The values of the weight parameters for synaptic connections are provided in Table 2, and the other parameter values are listsed in Table 3. Simulation codes required to produce each figure are available at https://github.com/zhuojunyu-appliedmath/Powerstroke-recovery.
RG-F | RG-E | In-F | In-E | PF-F | PF-E | Mn-F | Mn-E | Int | Inab-E | |
Drive connections, | ||||||||||
Supra-spinal drive, | 0.08 | 0.08 | 0.40 | 0.40 | ||||||
External drive, | 0.18 | |||||||||
Excitatory connections, | ||||||||||
RG-F | 0.41 | 0.70 | ||||||||
RG-E | 0.41 | 0.70 | ||||||||
PF-F | 1.95 | |||||||||
PF-E | 1.30 | 0.35 | ||||||||
Inab-E | 0.82 | |||||||||
Inhibitory connections, | ||||||||||
In-F | 2.20 | 6.60 | ||||||||
In-E | 2.20 | 6.60 | 2.80 | |||||||
Int | 0.55 | |||||||||
Feedback connections, | ||||||||||
Ia-F | 0.06 | 0.27 | 0.19 | |||||||
II-F | 0.0348 | 0.1566 | 0.1102 | |||||||
Ia-E | 0.06 | 0.44 | 0.10 | 0.16 | ||||||
Ib-E | 0.066 | 0.484 | 0.11 | 0.176 |
Parameter | Value | Unit | Parameter | Value | Unit |
---|---|---|---|---|---|
20 | pF | 585 | Nmm | ||
55 | mV | 60 | mm | ||
-80 | mV | 7 | mm | ||
for RG, PF, Mn | -64 | mV | for flexor | 72.5 | N |
for others | -60 | mV | for extensor | 37.7 | N |
-10 | mV | 2.3 | none | ||
-70 | mV | 1.6 | none | ||
for RG | 3.5 | nS | 1.62 | none | |
for PF | 0.5 | nS | 68 | mm | |
for Mn | 0.3 | nS | -0.69 | none | |
4.5 | nS | 0.18 | none | ||
1.6 | nS | 0.17 | none | ||
10 | nS | 6.2 | none | ||
10 | nS | 2 | none | ||
-30 | mV | 0.06 | none | ||
-50 | mV | 0.026 | none | ||
for Mn | 3 | mV | 1.5 | none | |
for others | 8 | mV | 0.06 | none | |
1.4 | none | 0 | none | ||
as in Figure | none | 1 | none | ||
300 | g | 0.6 | none | ||
0.00981 | mm/ms2 | for Ia | 60.007 | mm | |
300 | mm | for II | 58.457 | mm | |
0.002 | gmm2/ms | 3.393 | N |
Appendix C Variational analysis
The tools in the variational analysis were developed in Wang et al., (2021) and generalized in Yu and Thomas, (2022) and Yu et al., (2023). Consider a parameterized continuous-time dynamical system defined on a domain ,
(25) |
where , , and is either smooth or piecewise smooth in both and . Suppose admits a family of hyperbolic and asymptotically attracting limit cycles for including , which we consider as the unperturbed limit cycle . We assume that the domain is partitioned into two subdomains with the subdomain boundaries transverse to the flow of the unperturbed limit cycle, and that is smooth within each subdomain.
C.1 Local timing response cuvre (lTRC)
For , , let be the time remaining until the unperturbed trajectory beginning at exits the region. By construction, along ,
Hence, by the chain rule
(26) |
where we write for . Define the local timing response curve (lTRC) for subdomain to be . Differentiating both sides of (26) with respect to yields
(27) |
where denotes the Jacobian of . Let denote the exit point of tracjectory from region and denote a normal vector of the exit boundary of at . Following (26),
which gives the boundary condition of ,
(28) |
The adjoint equation (27) together with the boundary condition (28) defines the lTRC within region .
One application of the lTRC is to calculate the duration the trajectory spent in each region (phase). For a small positive , consider and the corresponding perturbed trajectory We assume that
where and represent the duration of the unperturbed trajectory and the perturbed trajectory spent in phase , respectively, and is therefore the linear shift in the phase duration in response to the perturbation. In Yu et al., (2023), we provide a general expression for :
(29) |
Here, at time the unperturbed limit cycle enters phase and at time it exits the phase; and represent the entry point to and exit point from the phase , respectively, for the limit cycle trajectory with . The expression shows that the first-order shift in the phase duration consists of three terms: the first term accounts for the impact of the perturbation on the entry point to the region; the second term arises from the impact on the exit point; the integral term shows the impact on the vector field during the transit from ingress to egress.
C.2 Infinitesimal shape response curve (iSRC)
Consider the perturbed trajectory and the unperturbed trajectory . Within each subdomain , we expand :
where is introduced as a rescaled time coordinate so that the expansion is uniform with respect to . It satisfies
The vector function is defined as the infinitesimal shape response curve (iSRC), which is piecewise-specified with period . In region , the iSRC satisfies a nonhomogenous variational equation
(30) |
where measures the (local) timing sensitivity to the perturbation in region . In the special case of with a linear scaling, i.e., , we have independent of , where is given by (29).
The iSRC equation (30) requires an initial condition, which is
where and represent the intersection points of the limit cycle trajectories with the boundary surfaces between regions. Generally, changing the surfaces (initial conditions) results in the consequent iSRC functions related by a simple phase shift, which has no effect on the sensitivity given by (10). See Lemma 1 in Yu and Thomas, (2022) or Lemma 2.3 in Wang et al., (2021)) for proof.
C.3 Derivation of sensitivity formula (10)
This section presents the derivation of equation (10). It is a general formula to calculate the sensitivity of an averaged quantity with respect to the change in any parameter, applicable to the powerstroke-recovery systems. The derivation follows the same exposition given in Yu and Thomas, (2022) but with some modifications.
Denote the unperturbed limit cycle (with ) by and the perturbed limit cycle (with ) by , respectively. Denote the periods by
During the powerstroke, each limit cycle orbit is divided into steps of equal time or ; during the recovery, each orbit is divided into steps of equal time or . Let Mark off the points as in Figure 17. Then,
Note that at the powerstroke-to-recovery (pr) transition state,
and at the recovery-to-powerstroke (rp) transition state,
Let be the average of the quantity around the trajectory , i.e.,
(31) |
and similarly . Suppose and the limit cycles are piecewise smooth, where the discontinuity occurs at the transitions between the two phases. As , following (31), since the set of discontinuous points has measure zero, we have
with a similar expression for . Then,
Let the fraction of the powerstroke duration be , i.e.,
Expand and around :
where represents the linear shift in under the perturbation. Since the quantity may be explicitly affected by the perturbation during the powerstroke, then there is an additional term in the expansion of . Note that the gradient is not well-defined at the discontinuous points. Then,
Using the iSRC of , we obtain
Dividing both sides by yields
Taking the limits and , we obtain
For the last two terms, use the directional derivative in the direction (say ) tangent to the recovery-powerstroke boundary, and expand and to be
Then,
Therefore, we obtain the sensitivity of the average for the unperturbed system:
(32) |
Note that when represents the change rate of progress, then for the models considered here (and for Shaw et al.’s Aplysia feeding model (Shaw et al.,, 2015; Lyttle et al.,, 2017; Wang et al.,, 2022)) vanishes during the recovery, indicating that for and that the second integral in (32) reduces to zero. Formula (10) is derived.
References
- Alfaro et al., (2010) Alfaro, V., Vilanova, R., Méndez, V., and Lafuente, J. (2010). Performance/robustness tradeoff analysis of PI/PID servo and regulatory control systems. In 2010 IEEE International Conference on Industrial Technology, pages 111–116. IEEE.
- Baylis, (1966) Baylis, L. E. (1966). Living control systems. Freeman.
- Boulet and Duan, (2007) Boulet, B. and Duan, Y. (2007). The fundamental tradeoff between performance and robustness-a new perspective on loop shaping-classic control revisited part II. IEEE Control Systems Magazine, 27(3):30–44.
- Brown et al., (2004) Brown, E., Moehlis, J., and Holmes, P. (2004). On the phase reduction and response dynamics of neural oscillator populations. Neural Computation, 16(4):673–715.
- Brown, (1911) Brown, T. G. (1911). The intrinsic factors in the act of progression in the mammal. Proceedings of the Royal Society of London. Series B, containing papers of a biological character, 84(572):308–319.
- Brown, (1914) Brown, T. G. (1914). On the nature of the fundamental activity of the nervous centres; together with an analysis of the conditioning of rhythmic activity in progression, and a theory of the evolution of function in the nervous system. The Journal of Physiology, 48(1):18.
- Chen et al., (2021) Chen, J., Yin, B., Wang, C., Xie, F., Du, R., and Zhong, Y. (2021). Bioinspired closed-loop CPG-based control of a robot fish for obstacle avoidance and direction tracking. Journal of Bionic Engineering, 18:171–183.
- Chiel and Beer, (1997) Chiel, H. J. and Beer, R. D. (1997). The brain has a body: adaptive behavior emerges from interactions of nervous system, body and environment. Trends in Neurosciences, 20(12):553–557.
- Diekman et al., (2017) Diekman, C. O., Thomas, P. J., and Wilson, C. G. (2017). Eupnea, tachypnea, and autoresuscitation in a closed-loop respiratory control model. Journal of Neurophysiology, 118(4):2194–2215.
- Ermentrout and Terman, (2010) Ermentrout, B. and Terman, D. H. (2010). Mathematical foundations of neuroscience, volume 35. Springer.
- Galvanetto and Bishop, (1999) Galvanetto, U. and Bishop, S. R. (1999). Dynamics of a simple damped oscillator undergoing stick-slip vibrations. Meccanica, 34:337–347.
- Golubitsky and Stewart, (2017) Golubitsky, M. and Stewart, I. (2017). Homeostasis, singularities, and networks. Journal of Mathematical Biology, 74:387–407.
- Golubitsky and Stewart, (2018) Golubitsky, M. and Stewart, I. (2018). Homeostasis with multiple inputs. SIAM Journal on Applied Dynamical Systems, 17(2):1816–1832.
- Golubitsky and Wang, (2020) Golubitsky, M. and Wang, Y. (2020). Infinitesimal homeostasis in three-node input–output networks. Journal of Mathematical Biology, 80(4):1163–1185.
- Grodins et al., (1967) Grodins, F. S., Buell, J., and Bart, A. J. (1967). Mathematical analysis and digital simulation of the respiratory control system. Journal of Applied Physiology, 22(2):260–276.
- Harris-Warrick and Cohen, (1985) Harris-Warrick, R. M. and Cohen, A. H. (1985). Serotonin modulates the central pattern generator for locomotion in the isolated lamprey spinal cord. Journal of Experimental Biology, 116(1):27–46.
- Hutter et al., (2014) Hutter, M., Gehring, C., Höpflinger, M. A., Blösch, M., and Siegwart, R. (2014). Toward combining speed, efficiency, versatility, and robustness in an autonomous quadruped. IEEE Transactions on Robotics, 30(6):1427–1440.
- Izhikevich and Ermentrout, (2008) Izhikevich, E. M. and Ermentrout, B. (2008). Phase model. Scholarpedia, 3(10):1487.
- Jahn and Votta, (1972) Jahn, T. L. and Votta, J. J. (1972). Locomotion of protozoa. Annual Review of Fluid Mechanics, 4(1):93–116.
- Jelbart and Wechselberger, (2020) Jelbart, S. and Wechselberger, M. (2020). Two-stroke relaxation oscillators. Nonlinearity, 33(5):2364.
- Katz, (2023) Katz, P. S. (2023). Conclusion and perspectives: What convergent evolution of animal forms and functions says about the predictability of evolution. In Convergent Evolution: Animal Form and Function, pages 581–594. Springer.
- Khoo, (2018) Khoo, M. C. (2018). Physiological control systems: analysis, simulation, and estimation. John Wiley & Sons.
- Korkmaz et al., (2021) Korkmaz, D., Ozmen Koca, G., Li, G., Bal, C., Ay, M., and Akpolat, Z. H. (2021). Locomotion control of a biomimetic robotic fish based on closed loop sensory feedback CPG model. Journal of Marine Engineering & Technology, 20(2):125–137.
- Kuo, (2002) Kuo, A. D. (2002). The relative roles of feedforward and feedback in the control of rhythmic movements. Motor Control, 6(2):129–145.
- Lee and Tomizuka, (1996) Lee, H. S. and Tomizuka, M. (1996). Robust motion controller design for high-accuracy positioning systems. IEEE Transactions on Industrial Electronics, 43(1):48–55.
- Lyttle et al., (2017) Lyttle, D. N., Gill, J. P., Shaw, K. M., Thomas, P. J., and Chiel, H. J. (2017). Robustness, flexibility, and sensitivity in a multifunctional motor control model. Biological Cybernetics, 111:25–47.
- Markin et al., (2010) Markin, S. N., Klishko, A. N., Shevtsova, N. A., Lemay, M. A., Prilutsky, B. I., and Rybak, I. A. (2010). Afferent control of locomotor CPG: insights from a simple neuromechanical model. Annals of the New York Academy of Sciences, 1198(1):21–34.
- Mo et al., (2023) Mo, A., Izzi, F., Gönen, E. C., Haeufle, D., and Badri-Spröwitz, A. (2023). Slack-based tunable damping leads to a trade-off between robustness and efficiency in legged locomotion. Scientific Reports, 13(1):3290.
- Morris and Lecar, (1981) Morris, C. and Lecar, H. (1981). Voltage oscillations in the barnacle giant muscle fiber. Biophysical Journal, 35(1):193–213.
- Pearson, (1985) Pearson, K. (1985). Are there central pattern generators for walking and flight in insects? In Feedback and motor control in invertebrates and vertebrates, pages 307–315. Springer.
- Prochazka, (1999) Prochazka, A. (1999). Quantifying proprioception. Progress in Brain Research, 123:133–142.
- Prochazka and Gorassini, (1998) Prochazka, A. and Gorassini, M. (1998). Models of ensemble firing of muscle spindle afferents recorded during normal locomotion in cats. The Journal of Physiology, 507(1):277–291.
- Ronsse et al., (2008) Ronsse, R., Thonnard, J.-L., Lefevre, P., and Sepulchre, R. (2008). Control of bimanual rhythmic movements: trading efficiency for robustness depending on the context. Experimental Brain Research, 187:193–205.
- Rubin and Terman, (2002) Rubin, J. E. and Terman, D. (2002). Geometric singular perturbation analysis of neuronal dynamics. In Handbook of dynamical systems, volume 2, pages 93–146. Elsevier.
- Sariyildiz and Ohnishi, (2013) Sariyildiz, E. and Ohnishi, K. (2013). Performance and robustness trade-off in disturbance observer design. In IECON 2013-39th Annual Conference of the IEEE Industrial Electronics Society, pages 3681–3686. IEEE.
- Schultheiss et al., (2012) Schultheiss, N. W., Prinz, A. A., and Butera, R. J. (2012). Phase Response Curves in Neuroscience. Springer.
- Sharbafi et al., (2020) Sharbafi, M. A., Yazdanpanah, M. J., Ahmadabadi, M. N., and Seyfarth, A. (2020). Parallel compliance design for increasing robustness and efficiency in legged locomotion–theoretical background and applications. IEEE/ASME Transactions on Mechatronics, 26(1):335–346.
- Shaw et al., (2015) Shaw, K. M., Lyttle, D. N., Gill, J. P., Cullins, M. J., McManus, J. M., Lu, H., Thomas, P. J., and Chiel, H. J. (2015). The significance of dynamical architecture for adaptive responses to mechanical loads during rhythmic behavior. Journal of Computational Neuroscience, 38:25–51.
- Skinner et al., (1994) Skinner, F. K., Kopell, N., and Marder, E. (1994). Mechanisms for oscillation and frequency control in reciprocally inhibitory model neural networks. Journal of Computational Neuroscience, 1:69–87.
- Smith et al., (1991) Smith, J. C., Ellenberger, H. H., Ballanyi, K., Richter, D. W., and Feldman, J. L. (1991). Pre-Bötzinger complex: a brainstem region that may generate respiratory rhythm in mammals. Science, 254(5032):726–729.
- Spardy et al., (2011) Spardy, L. E., Markin, S. N., Shevtsova, N. A., Prilutsky, B. I., Rybak, I. A., and Rubin, J. E. (2011). A dynamical systems analysis of afferent control in a neuromechanical model of locomotion: I. rhythm generation. Journal of Neural Engineering, 8(6):065003.
- Sutton et al., (2023) Sutton, G., Szczecinski, N., Quinn, R., and Chiel, H. (2023). Phase shift between joint rotation and actuation reflects dominant forces and predicts muscle activation patterns. PNAS Nexus, 2(10):pgad298.
- Terman et al., (1998) Terman, D., Kopell, N., and Bose, A. (1998). Dynamics of two mutually coupled slow inhibitory neurons. Physica D: Nonlinear Phenomena, 117(1-4):241–275.
- Vasconcelos et al., (2009) Vasconcelos, J., Athans, M., Fekri, S., Silvestre, C., and Oliveira, P. (2009). Stability-and performance-robustness tradeoffs: Mimo mixed- vs complex- design. International Journal of Robust and Nonlinear Control: IFAC-Affiliated Journal, 19(3):259–294.
- Wang et al., (2021) Wang, Y., Gill, J. P., Chiel, H. J., and Thomas, P. J. (2021). Shape versus timing: linear responses of a limit cycle with hard boundaries under instantaneous and static perturbation. SIAM Journal on Applied Dynamical Systems, 20(2):701–744.
- Wang et al., (2022) Wang, Y., Gill, J. P., Chiel, H. J., and Thomas, P. J. (2022). Variational and phase response analysis for limit cycles with hard boundaries, with applications to neuromechanical control problems. Biological Cybernetics, 116(5-6):687–710.
- Yao et al., (1997) Yao, B., Al-Majed, M., and Tomizuka, M. (1997). High-performance robust motion control of machine tools: an adaptive robust control approach and comparative experiments. IEEE/ASME Transactions on Mechatronics, 2(2):63–76.
- Yu et al., (1999) Yu, S.-N., Crago, P. E., and Chiel, H. J. (1999). Biomechanical properties and a kinetic simulation model of the smooth muscle I2 in the buccal mass of Aplysia. Biological Cybernetics, 81:505–513.
- Yu et al., (2023) Yu, Z., Rubin, J. E., and Thomas, P. J. (2023). Sensitivity to control signals in triphasic rhythmic neural systems: a comparative mechanistic analysis via infinitesimal local timing response curves. Neural Computation, 35(6):1028–1085.
- Yu and Thomas, (2021) Yu, Z. and Thomas, P. J. (2021). Dynamical consequences of sensory feedback in a half-center oscillator coupled to a simple motor system. Biological Cybernetics, 115(2):135–160.
- Yu and Thomas, (2022) Yu, Z. and Thomas, P. J. (2022). A homeostasis criterion for limit cycle systems based on infinitesimal shape response curves. Journal of Mathematical Biology, 84(4):1–23.
- Zhang and Lewis, (2013) Zhang, C. and Lewis, T. J. (2013). Phase response properties of half-center oscillators. Journal of Computational Neuroscience, 35:55–74.