Variational and phase response analysis for limit cycles with hard boundaries, with applications to neuromechanical control problems

In many animals, motor control involves neural oscillatory circuits that can produce rhythmic patterns of neural activity without receiving rhythmic inputs [central pattern generators (CPGs)], force generation by muscles, and interactions between the body and environment. Moreover, sensory feedback from the peripheral nervous system is known to modulate the rhythms of the electrical signals in CPGs and thereby facilitate adaptive behavior.

Motor systems show an overall robustness, but because they are highly nonlinear, understanding how they achieve robustness due to their different components is difficult. To understand how animals produce robust behavior in a variable environment, Shaw et al. (2015) and Lyttle et al. (2017) developed a neuromechanical model of triphasic motor patterns to describe the feeding behavior of the marine mollusk Aplysia californica. Like many rhythmic motor systems, feeding in Aplysia involves two distinct phases of movement: a power stroke during which the musculature engages with the substrate (the seaweed) against which it exerts a force to advance its goal (ingestion of seaweed), and a recovery stage during which the motor system disengages from the substrate to reposition itself, in preparation for beginning the next power stroke. Similarly, in legged locomotion, the stance phase corresponds to the power stroke and the swing phase corresponds to the recovery stage.

Also, like many rhythmic motor systems, feeding in Aplysia involves a closed-loop system, which integrates biomechanics and sensory feedback, and exhibits a stable limit cycle solution. It has been conjectured that sensory feedback plays a crucial role in creating robust behavior by extending or truncating specific phases of the motor pattern (Lyttle et al. 2017, §3.1). To test this hypothesis, we applied small mechanical perturbations as well as parametric perturbations to the sensory feedback pathways in the coupled neuromechanical system. It was shown in Lyttle et al. (2017) that a sustained increase in mechanical load leads to changes in both shape and timing of the limit cycle solution: the system generates stronger retractor muscle force for a longer time in response to the increased load. Qualitatively similar effects have been observed during in vivo experiments in Aplysia (Gill and Chiel 2020). In general, we expect that applying parametric changes to CPG-based motor systems leads to changes in both the shape and timing of the resulting limit cycle behavior (Fig. 1).

In Aplysia, the increased duration (timing) and increased force (shape) have opposite effects on the task-fitness, measured as seaweed consumption per unit time. Strengthening the retractor force pulls in more food with each cycle, which increases fitness, whereas lengthening the cycle time decreases fitness. Together these effects approximately cancel, making the system robust against increased load. This type of “stronger-and-longer” response may occur generically in other motor systems. Thus, in this paper, we seek to understand the roles of sensory feedback and biomechanics in enhancing robustness. To this end, we apply recently developed tools from variational analysis (Wang et al. 2021) to quantitatively study changes in shape and timing of a limit cycle under static perturbations.

A sustained change in parameter \(\Theta \) in a dynamical system \({\dot{\textbf{x}}}=\textbf{F}(\textbf{x},\Theta )\) producing a limit cycle trajectory typically causes changes in both the timing and shape of the trajectory, which may both influence the performance S of the limit cycle system

Full size image

In the first part of the present paper (cf. Sect. 2.1), we apply the classical tools of forward variational analysis to the model introduced by Shaw et al. (2015) and Lyttle et al. (2017) (denoted as the Shaw–Lyttle–Gill or SLG model for brevity) to arrive at the following insights:

• Nonlinear biomechanical properties confer resilience by immediately increasing resistance to applied loads, on timescales much faster than neural responses.

• The main effect of sensory feedback is to shift the timing of retraction neural pool deactivation; in parallel, firing rate saturation effectively censors sensory feedback during specific movement subintervals.

While the forward-in-time variational analysis is illuminating and allows us to explain in detail the robustness mechanism, it is still incomplete. Over time, the original and perturbed cycle will become increasingly out of phase due to the timing changes under sustained perturbations. Hence, the shape displacements estimated from the forward variational analysis will become less and less accurate over time. This difficulty is not limited to models of feeding in Aplysia californica. For example, if we were to compare the gaits of two subjects walking on treadmills with slightly different speeds, although the ratio of stance and swing may be the functionally important aspect, this quantity is difficult to assess directly without putting the two movements on a common footing by comparing them using a common time scale.

Thus, in order to compare perturbed and unperturbed motions with greater accuracy, in the remainder of the paper (cf. Sect. 2.2 and following) we show how to extend the local-in-time variational analysis to a global analysis by applying the infinitesimal shape response curve (ISRC) analysis and local timing response curve (LTRC) analysis developed in Wang et al. (2021). We review these methods in Sect. 3. This time-rescaled analysis accounts for both global timing sensitivity (through the infinitesimal phase response curve, IPRC), as well as local timing sensitivity (through LTRC) by rescaling time to take into account local differences in the effects of parametric variation. It yields a more accurate and self-consistent description of the oscillator trajectory’s changing shape in response to parametric perturbations and helps complete the picture by providing a complementary perspective. Specifically, our time-rescaled analysis provides additional insights, specifically that

• Increasing the applied load on the system increases the duty cycle of the neuron pool responsible for retraction, in the sense that the retraction neuron pool is activated for a larger percentage of the closed phase of the cycle. (The closed phase of the trajectory occurs while the animal’s radula-odontophore, or grasper, is closed on the seaweed, and encompasses the power stroke.) This effect ultimately results in more seaweed being consumed, despite increased force opposing ingestion.

• We are able to derive the multidimensional infinitesimal phase response curve (IPRC) despite the presence of nonsmooth dynamics in the system; we identify the mechanical component of the IPRC as the one that contributes most to robustness, and note that its contribution arises from the “power stroke” segment of the motor cycle.

• We derive an analytical expression for the robustness to the mechanical perturbation that decomposes naturally into a sum of two terms, one capturing the effect of the perturbation on the shape of the trajectory, and the other capturing the effect on the timing; this result provides a quantitative analysis of robustness that confirms the qualitative insights described previously in the literature.

• In addition to sensory feedback and intrinsic biomechanical properties, robustness against changes in applied load can arise from coordinated changes of multiple parameters such as the gain of sensory feedback and muscle stiffness.

The dynamics of the SLG model (Lyttle et al. 2017) is given by

This system incorporates firing rates of three neuron populations, corresponding to the “protraction-open” (\(a_0\)), “protraction-closed” (\(a_1\)), and “retraction” phase (\(a_2\)). Note that when a nerve cell ceases firing because of inhibition, its firing rate will be held at zero until the balance of inhibition and excitation allow firing activity to resume. Hence, we supplement model (1) with three hard boundaries introduced by the requirement that the firing rates \(a_0,a_1, a_2\) must be nonnegative:

During the limit cycle, when a neural variable \(a_i\) changes from positive to 0, we call that the \(a_i\) landing point; when it changes from 0 to positive, we call that the \(a_i\) liftoff point. The fact that the trajectory is non-smooth at the landing and liftoff points will play an important role in the analysis to follow. See Sect. 4 for further discussion of the biological basis for our modeling assumptions.

This model also consists of a simplified version of the mechanics of the feeding apparatus: a grasper that can open or close (\(x_r\)), a muscle that can protract the grasper to reach the food (\(u_0\)) and another muscle that can retract the grasper to pull the food back into its mouth (\(u_1\)). The net force exerted by the muscles is given by the sum of the two muscle forces

where

is the effective length-tension curve for muscle forces, \(c_i\), \(w_i\) and \(k_i\) denote the mechanical properties of each muscle.

\(F_{\text {sw}}\) represents the external force applied to the seaweed, which can only be felt by the grasper when it is closed on the food (\(a_1+a_2>0.5\)), during which \(r=1\). When the grasper is open (\(a_1+a_2\le 0.5\)), \(r=0\); that is, the grasper moves independently of the seaweed. This condition leads to a transversal crossing boundary; that is, the open/closing boundary given by

Values for model parameters and initial conditions are given in Tables 2 and 3 in Appendix. For additional details on the biological assumptions motivating the model, see Shaw et al. (2015) and Lyttle et al. (2017).

In this paper, we are interested in the so-called heteroclinic mode in which the neural dynamics temporarily slow down when the sensory feedback overcomes the endogenous neural excitation and forces the neural trajectory to slide along the hard boundary \(a_i=0\). Such temporary slowing down of neural variables allows the muscles, which evolve on slower timescales, to “catch up”; hence the seaweed can be swallowed and ingested successfully. Following the terminology from Wang et al. (2021), we call attracting periodic trajectories that experience sliding motions limit cycles with sliding components (LCSCs).

Using classical sensitivity analysis and our recently developed tools from variational analysis (Wang et al. 2021), we show here that biomechanics and sensory feedback cooperatively support strong robustness by changing the timing and shape of the neuromechanical trajectory. While both sensory feedback and biomechanics respond immediately to the increased load, we find that the sensory feedback effect is initially censored while the neural activity is pinned against a hard boundary of neuronal firing. Thus the effect of the sensory feedback signal is significantly delayed relative to onset of the increased load. Our analysis suggests that sensory feedback mediates robustness mainly through shifting the timing of neural activation and specifically increasing the duty cycle of the retraction neural pool. This response allows the system to digest more seaweed despite the increased force opposing ingestion and hence achieve strong robustness. In addition to uncovering the mechanisms for robust motor control, our methods allow us to quantify analytically the robustness of the model system to the mechanical perturbation. Finally, although we focus on the Aplysia californica feeding system as our working example, our methods should extend naturally to a broad range of motor systems.

Our paper is organized as follows. We present our analysis and main results in Sect. 2. Methods that we use to understand the robustness in the Aplysia neuromechanical model are presented and reviewed in Sect. 3. We discuss limitations and possible extensions of our approach in Sect. 4.

Recall that a sustained (parametric) perturbation often causes changes in both shape and timing of the neuromechanical trajectory solution of (1). In this paper, we adopt methods developed in Wang et al. (2021) for analyzing the joint variation of both shape and timing of limit cycles with sliding components under parametric perturbations.

2.1 Forward variational analysis

We begin our analysis by investigating how the shape of the trajectory changes in response to a sustained increase in mechanical load \(F_{\text {sw}}\). To a first approximation, the change in shape can be captured by classical sensitivity analysis (also called forward variational analysis) which we review in Sect. 3.

We apply a small static perturbation to system (1) by increasing the model parameter \(F_{\text {sw}}\) by \(\varepsilon =0.02\): \(F_{\text {sw}}\rightarrow F_{\text {sw}}+\varepsilon \), and comparing the new, perturbed limit cycle trajectory \(\gamma _\varepsilon \) to the original, unperturbed limit cycle trajectory \(\gamma _0\), beginning at the start of the grasper-closed phase (time 0 in Fig. 2 panels A–D). That is, we plot \({\textbf{v}}\approx (\gamma _\varepsilon (t)-\gamma _0(t))/\varepsilon \); see (13)–(15) for precise definitions. Note that \(F_{\text {sw}}\) is multiplied by an indicator function that is only nonzero when the trajectory is in the grasper-closed phase. Thus the perturbation is only present when the grasper is closed on the food.

A small sustained perturbation is applied to the Aplysia model (1) over the closing phase in which \(F_{\text {sw}}\rightarrow F_{\text {sw}}+\varepsilon \) with perturbation \(\varepsilon =0.02\). A, B Time series of the trajectory components for nominal force value \(F_{\text {sw}}\) (solid) and perturbed force value (dashed) over two periods, aligned at the start of the closed phase. C, D The displacement solution, \({\textbf{v}}(t)\), to the forward variational equation over two periods. The gray curve in D denotes the displacement between the perturbed and unperturbed muscle forces \(F_{\textrm{musc}}\), shown as the gray curves in panel B. The yellow dashed line in panel D approximates the displacement in \(x_r\) if the net muscle force \(F_{\textrm{musc}}\) did not change after perturbation. The intervals during which the grasper without perturbation is closed on the food are indicated by the shaded regions. The vertical magenta lines indicate the times at which the grasper under perturbation switches from closed to open. The difference in periods and the delay in the grasper opening time both accumulate, making the comparison between the two trajectories invalid except at short times. A, C show trajectories and displacements along the neural directions, while B, D show trajectories and displacements along the mechanical directions. The lines in panel C (resp., D) approximate the difference between the dashed and solid trajectories from panel A (resp., B) per perturbation

Full size image

The neural and biomechanical components of the unperturbed trajectory \(\gamma (t)\) are shown by the solid curves in Fig. 2A, B, respectively. The perturbed trajectories \(\gamma _\varepsilon (t)\) are indicated by the dashed lines. The gray shaded regions indicate phases when the grasper in the unperturbed system is closed. With the perturbation (increased load), the transition from closing to opening is delayed; this transition is indicated by the magenta vertical line. Figure 2C, D shows the difference between the two trajectories per perturbation along the neural directions and along the biomechanical directions, respectively. These curves can be approximated by the solutions to the forward variational equation (14) defined in Sect. 3. The muscle forces \(F_{\textrm{musc}}(u_0, u_1, x_r)\) before and after the perturbation are shown by the gray curves in Fig. 2B, and the difference between them is included as the gray curve in Fig. 2D.

Figure 2 yields several insights about the roles of sensory feedback and biomechanics in robustness, which we discuss in detail below.

2.1.1 Biomechanics confer resilience by immediately increasing resistance to the increased load

Immediately following the perturbation of the mechanical load, we observed a positive displacement in the grasper position \(x_r\) relative to the unperturbed trajectory (Fig. 2D, yellow curve). This displacement simply reflects the grasper being pulled by a stronger force \(F_{\text {sw}}+\varepsilon \). If nothing other than the applied load \(F_{\text {sw}}\) changes in the system, a linearized analysis suggests that the initial displacement of \(x_r\) would approximately follow the yellow dashed line given by \(\frac{1}{b_r}t\) (Fig. 2D). Nonetheless, while this line gives a good initial approximation, the true displacement in \(x_r\) (yellow solid curve) quickly sags below the yellow dashed line over time. This difference arises due to the negative displacement occurring in the muscle force \(F_{\textrm{musc}}(u_0, u_1,x_r)\) (Fig. 2D, gray curve) which acts to overcome the increased load. However, early in the retraction cycle, all other variables, including the muscle activation \(u_0\) and \(u_1\), show almost no displacement at all (see Fig. 2C, D). This observation suggests that long before the sensory feedback effect has time to act, the biomechanics may play an essential role in generating robust motor behavior, by providing an immediate, short-term resistance to increased load.

Early in the retraction cycle, increasing the load stretches both the retractor and protractor muscles, and moves them down their length-tension curves. As a result, both forces become weaker, but the magnitude of the protractor muscle force drops more quickly than the retractor muscle force (Fig. 3). Thus, the retractor muscle force grows relative to the opposing protractor muscle force. This shift endows the system with a built-in resilience, in that increasing seaweed force opposing ingestion automatically (i.e., without changes in neural activation) engages a larger resisting muscle force long before neural pools or muscles show differential activation. This is a new insight beyond the previous “longer-stronger” hypothesis (Shaw et al. 2015; Lyttle et al. 2017).

2.1.2 Sensory feedback effects are largely delayed by the firing rate hard boundary properties

Changes in \(x_r\) due to the increased load will immediately propagate to the neural variables (\(a_0, a_1, a_2\)) through the sensory feedback \(\varepsilon _i(x_r-\xi _i)\sigma _i)/\tau _a\), and hence should affect the neural variables without any lag. However, no significant displacement of the neural variables is observed until nearly the end of the retraction cycle (Fig. 2C). In other words, while the sensory feedback itself immediately responds to the increased load, the effect of the changed sensory feedback signal is not manifest until much later in the retraction cycle, when the protraction-open neuron pool is released from inhibition along its hard boundary and starts to fire (Fig. 2). Hence, the nonsmooth hard boundary conditions on neuronal firing rates significantly delay the effect of sensory feedback, and create intervals during which neurons are insensitive to sensory feedback.

The time series of the perturbed (dashed) and unperturbed (solid) protractor muscle forces \(F_{\textrm{musc,pro}}\) (red) and retractor muscle forces \(F_{\textrm{musc,ret}}\) (blue) over two periods. The gray shaded regions and the magenta lines have the same meanings as in Fig. 2

Full size image

2.1.3 Sensory feedback contributes by shifting the timing of neural activation

Due to the hard boundary effects, the displacements of the neural variables appear near the end of the closed phase, when the protraction-open neuron pool \(a_0\) lifts off from its hard boundary, and the retraction-closed neuron pool \(a_2\) deactivates to stop firing (Fig. 2A, C). A positive (resp., negative) displacement of \(a_2\) (resp., \(a_0\)) indicates that \(a_2\) deactivates (resp., \(a_0\) activates) at a later time with the increased load, and hence the retraction-closed phase is prolonged. Consequently, the retraction muscle activity will increase, because its stimulation by the retraction motor neuron is prolonged, allowing the slow retractor muscle to generate larger forces (Fig. 3). Similarly, we also observe a decreased protractor muscle activity, as the protraction neuron pool \(a_0\) turns on at a later time. This decrease leads to a stronger retractor muscle force and a weaker protractor muscle force (Fig. 3). Hence, a more negative net muscle force results (Fig. 2D, gray curve), which corresponds to a stronger resisting force pulling the seaweed towards the jaw to swallow the food. To summarize, the main effect of sensory feedback that contributes to robustness is prolonging the retraction phase to confer on the system a resilience in that increasing seaweed force opposing ingestion engages a larger resisting muscle force.

Thus, sensory feedback contributes to robustness primarily by shifting the timing of neural activation, as opposed to the magnitude of neural activation. Biologically, this distinction corresponds to affecting the timing of motor neuron burst onset or offset, rather than burst intensity.

2.2 Variational analysis with rescaled time—ISRC

Under the forward-in-time analysis, the grasper of the perturbed system lags behind the grasper of the unperturbed system throughout the closing phase; yet the net seaweed movement was measured to be greater for the perturbed system (Lyttle et al. 2017). Fig. 4 shows an expanded view of the perturbed and unperturbed systems’ grasper position (Fig. 4A) and the linearized difference produced by the variational equation (Fig. 4B). As this detailed view shows, at the time when the unperturbed system transitions from closed to open (gray-white boundary) the unperturbed grasper position is more retracted than the perturbed grasper position at the coincident time point. Similarly, the grasper component of the variational equation is positive at the gray-white boundary. Furthermore, at the time when the perturbed system transitions from closed to open (magenta line) the perturbed system continues to be less retracted than the unperturbed system. Thus, whether we compare the systems at the perturbed or unperturbed opening time, the perturbed grasper is “further behind.” Yet the overall effect in the perturbed system is a larger net intake of seaweed per cycle.

Enlarged views of Fig. 2B, D near the first transition from closed to open. Note the displacement of \(x_r\) (i.e., the \(x_r\) component of \(\textbf{v}(t)\) shown as the yellow solid line in B) is positive both at the end of the closing time of the unperturbed system (see the green arrow near the grey/white boundary) and the perturbed system (see the green arrow near the magenta line)

Full size image

This apparent contradiction underscores the need to extend the perturbation analysis beyond the standard forward-in-time variational analysis. In particular, if one cycle is slower than another, then while the local perturbation analysis can explain the cause-and-effect relations a short time into the future, they cannot account for the net effect around a cycle in a self-consistent way. Over time, the displacements between the two trajectories grow, and the linearized approximation becomes invalid except at short times (cf. Jordan et al. 2007). Hence, unless time is rescaled to take into account the difference in cycle period, comparing the components of the original and perturbed cycles will become less and less meaningful.

To overcome this difficulty, we extend the local-in-time variational analysis to a global analysis by rescaling time so the unperturbed closing and opening events coincide with those after perturbations, respectively. We do so by applying the infinitesimal shape response curve (ISRC) analysis and the local timing response curve (LTRC) (Wang et al. 2021), which we review in Sect. 3. This method yields a more accurate and self-consistent description of the oscillator trajectory’s changing shape in response to parametric perturbations (see Fig. 5). We show that the combination of the ISRC and the LTRC gives a sensitivity analysis of an oscillator to sustained perturbations within any given region (e.g., protraction or retraction cycle, opening or closing phase) and provides a self-contained framework for analytically quantifying and understanding robustness to perturbations.

Variational analysis with piecewise uniform time rescaling. The same sustained perturbation as in Fig. 2 is applied to the Aplysia model (1). A, B Time series of the perturbed (dashed) and unperturbed solutions (solid). Here piecewise uniform rescaling is applied so the closing and opening events coincide. C, D The ISRC with piecewise uniform rescaling \(\gamma (t)\) over two periods. Shaded regions have the same meanings as in Fig. 2. Note the \(x_r\) component of the ISRC is negative at the time of opening (see green arrow). With piecewise uniform rescaling, the variational approximation is consistent across multiple periods (c.f., Fig. 2)

Full size image

We write \(\gamma _1\) for the linear shift in the limit cycle shape in response to the static perturbation \(F_{\text {sw}}\rightarrow F_{\text {sw}}+\varepsilon \), that is:

uniformly in time. Note that the time for the perturbed trajectory is rescaled to be \(\tau _\varepsilon (t)\) to match the unperturbed time points. The linear shift \(\gamma _1(t)\) is the so-called ISRC curve and satisfies a nonhomogeneous variational equation (see Sect. 3). Compared with the forward variational equation, the ISRC equation has one additional nonhomogeneous term \(\nu _1 F_0(\gamma _0(t))\) that arises from the time rescaling. In this term, \(\nu _1\) is determined by the choice of time rescaling \(\tau _\varepsilon (t)\) and \(F_0(\gamma _0(t))\) is the unperturbed vector field evaluated along the unperturbed limit cycle \(\gamma _0(t)\) (see Sect. 3 for details).

Since the perturbation is applied to the seaweed, it can only be felt by the system when the grasper is closed on the seaweed. It is natural to expect that the segment at the closing phase has a different timing sensitivity than the segment at the opening phase. We hence choose to rescale time differently in the two phases, using piecewise uniform rescaling when computing the ISRC. This leads to a piecewise ISRC equation, where \(\nu _1\) is piecewise constant. It was shown in Wang et al. (2021) that \(\nu _1\) can be estimated from the LTRC analysis (see Sect. 3).

In Fig. 5A, B, the time traces of variables along the unperturbed limit cycle are shown by the solid curves, whereas the perturbed limit cycle whose time has been rescaled to match the unperturbed time points as described above are indicated by the dashed curves. With the piecewise rescaling, the transitions between the closing and opening events of the perturbed and unperturbed systems now coincide. The relative displacements between the perturbed and unperturbed trajectories are approximately given by the piecewise ISRC \(\gamma _1\) shown in Fig. 5C, D. In contrast to the forward variational analysis, in which the displacements grow over time, the piecewise ISRC curve is periodic, meaning we have achieved a self-consistent global description of the response of the limit cycle to increased load.

We now show that the apparent contradiction that we obtained from the forward variational analysis, i.e., that the grasper displacement at the end of the closing phase is positive (cf. Fig.4), can now be resolved in the time-rescaled picture. In response to the perturbation, the relative displacement of the grasper position (the \(x_r\) component of \(\gamma _{1}\), denoted as \(\gamma _{1,x_r}\)) initially increases (i.e., the grasper becomes more and more protracted due to the increased load) and reaches its peak at about \(t=1.4\) (see Fig. 5D, yellow curve). Then, it starts to decrease and becomes negative at the time when the grasper opens. This means later in the retraction cycle, the perturbed grasper is less and less protracted than the unperturbed version and eventually become more retracted by the end of the closing phase (Fig. 5D, green arrow). In summary, the grasper perturbed by larger force begins “behind” the unperturbed version, but catches up around 60% of the way through the retraction phase (in relative time) and comes out “ahead” by the time both graspers open, consistent with having a larger net seaweed intake (Lyttle et al. 2017).

To understand what causes \(\gamma _{1,x_r}\) to be negative despite its initial big rise, we consider the effect of the perturbation on the neural pool through sensory feedback. In Fig. 5C, we observe positive displacements in \(\gamma _{1,a_2}\) (yellow curve) occurring both when the retraction neuron pool \(a_2\) activates and when it deactivates. These displacements indicate that with the increased load, the retraction neuron \(a_2\) activates earlier and turns off later relative to the unperturbed \(a_2\). In other words, increasing the applied load on the system increases the duty cycle of the neuron pool involved in retraction, i.e., the retraction neuron pool is activated for a larger percentage of the total cycle. As a result, the motor system recruits a larger retractor muscle force, as indicated by the positive displacement of the retractor muscle activation \(u_1\) during the closing phase (Fig. 5D, red curve). A similar increase in motor recruitment in response to increased external load has been observed in vivo (Gill and Chiel 2020). In the model, the stronger retraction force acts to impede the protraction of the grasper, and eventually pulls the grasper to a more retracted state. Thus the grasper displacement crosses zero and becomes negative at the end of the closing phase (Fig. 5D, green arrow).

Note that there is no perturbation during the opening phase (Fig. 5, white space). During this phase of the cycle, displacements slowly decay and are nearly zero by the time the grasper closes on the food again.

2.3 Timing responses to sustained perturbations of \(F_{\text {sw}}\).

2.3.1 Infinitesimal phase response curve

To understand the timing response of system (1) to increased load, we perform an IPRC analysis. Figure 6 shows the time traces of the IPRC curve over one cycle. As before, the shaded region indicates the phase when the grasper is closed.

IPRC for the Aplysia model. Grey shaded region indicates the period when the radula/odontophore is closed. On the bottom, the red, yellow and blue rectangles denote the protraction-closed, retraction-closed and protraction-open phases, respectively. The blue spike in the IPRC in the top panel occurs when the \(a_0\) variable “lifts off” from the \(a_0=0\) boundary; the red spike occurs when the \(a_1\) variable lifts off; the liftoff point for \(a_2\) is indicated with an arrow

Full size image

The IPRC curves associated with biomechanical variables are shown in Fig. 6, lower panel. In particular, the timing sensitivity of system (1) to the increased load on the grasper (\(F_{\text {sw}}\rightarrow F_{\text {sw}}+\varepsilon \)) can be estimated using the IPRC along the \(x_r\) direction, i.e., the yellow curve in the lower panel of Fig. 6. Since the perturbation only has effect during the closing phase, only the portion of \(z_{x_r}\) over the shaded region is relevant. This portion is strictly negative. Therefore, in response to the increased load considered above, the system undergoes phase delay, and hence the total period is prolonged. This finding is consistent with earlier results on the sensory feedback effect obtained from the variational analysis (see Sect. 2.2).

The linear shift in period can be estimated by evaluating the integral

where \(T_0\), \(T_\varepsilon \) are the periods before and after perturbation \(\varepsilon \) (see Sect. 3). For the perturbation on \(F_{\text {sw}}\), the derivative \(\frac{\partial F_\varepsilon (\gamma _0(t))}{\partial \varepsilon }\) equals \((\textbf{0},\frac{1}{b_r})^\intercal \) over the grasper-closed region, and equals \(\textbf{0}\) during the grasper-open region, where the first \(\textbf{0}\) is a \(5\times 1\) zero vector and the second \(\textbf{0}\) is a \(6\times 1\) vector. It then follows that

where \(\Theta _{\textrm{close}}\) denotes the grasper-closed phase.

Other IPRC curves in Fig. 6 indicate the timing sensitivity of the model to other perturbations and lead to several useful insights as well as testable predictions. For example,

• The IPRC curves are continuous except at the liftoff points (Fig. 6 top panel, blue and red spikes). While all three neural variables go through liftoff points, there is no large spike in \(z_{a_2}\) (yellow curve). The absence of a yellow spike and the fact that the red spike is larger than the blue spike, imply that the system has the highest timing sensitivity to perturbing \(a_1\) and intermediate timing sensitivity to \(a_0\), both of which are significantly higher than the sensitivity to \(a_2\) perturbations.

• Excitatory inputs to neural populations lead to phase advance and hence shorten the total period, because the IPRC curves associated with neural variables are mostly positive (Fig. 6 top panel).

• Most of the time the system is not sensitive to neural perturbations, but there also exist sensitive regions when the trajectory is not restricted to the hard boundaries (e.g., Fig. 6 top panel, blue and red spikes). For instance, the system has high timing sensitivity to perturbations of \(a_0\) late in the closing phase and to perturbations of \(a_1\) late in the opening phase, whereas sensory inputs are largely ignored early in the opening phase. This effect is a concrete example of differential penetrance, a striking feature of many biological systems in which some neural activity can vary greatly, with little effect on behavior, whereas in other circumstances, a small change in neural activity may have a very large impact on behavior (Chiel et al. 1999; Beer et al. 1999; Ye et al. 2006; Cullins et al. 2015).

• Increasing the protractor muscle activation \(u_0\) causes a phase delay early in the closing phase and late in the opening phase, and a phase advance otherwise. In contrast, increasing the retraction muscle activation \(u_1\) causes a phase advance early in the closing phase and late in the opening phase, and a phase delay otherwise. “Appendix B” discusses why the system has different timing sensitivities to muscle perturbations.

Although all three neural variables go through liftoff points, there is no large yellow spike in \(z_{a_2}\) (see Fig. 6). To understand this, we note that before \(a_0\) (resp., \(a_1\)) lifts off its hard boundary, there exists no inhibition from other neurons except for inhibitory sensory feedback. However, when \(a_2\) lifts off at around \(t\approx 3.2\), it still experiences inhibition from \(a_0\) (see Fig. 5A). In other words, there are two inhibitory effects pressing neurons \(a_2\) down to the hard boundary, but only one inhibitory effect acting on the other two neuron populations. As a result, while there is a discontinuous jump of the IPRC curve corresponding to \(a_2\) at the liftoff point, it remains small as the other inhibition is still present.

2.3.2 Local timing response curve

While the IPRC is a powerful tool for understanding the global timing sensitivity of an oscillator to sustained perturbations, it does not give local timing sensitivities, which, however, are needed for computing the ISRC curve as discussed above. We hence adopt the local timing response curve (LTRC) method developed in Wang et al. (2021) and reviewed in Sect. 3. To illustrate this method, we show the LTRC associated with the closing phase and denote it as \(\eta ^{\textrm{close}}\) (see Fig. 7). Although the LTRC \(\eta ^{\textrm{close}}\) is defined throughout the full domain, estimating the effect of the perturbation within the closing region only requires evaluating the LTRC in this region. Figure 7 shows the time series of \(\eta ^{\textrm{close}}\) for the model in the closing region, obtained by numerically integrating the adjoint equation backward in time with the initial condition of \(\eta ^{\textrm{close}}\) given by its value when the grasper switches from closing to opening. Note that \(\eta _{x_r}\), the yellow curve in Fig. 7 lower panel, remains positive over the closing phase. This implies that the increased load on seaweed prolongs the time remaining in the closing region; that is, the increased load prolongs the total closing time. The relative shift in the closing time caused by the increased load can also be estimated by integrating the LTRC (see Section 3).

In addition, Fig. 7 implies that strengthening the protractor muscle activation \(u_0\) during the closing phase prolongs the total closing time, whereas increasing the retraction muscle activation \(u_1\) decreases the total closing time. Similarly, we can compute the LTRC over other phases, such as the retraction phase, in order to estimate local timing sensitivities of the system in other regions.

Time series of the LTRC \(\eta ^{\textrm{close}}\) over the closing phase. The liftoff point on \(a_0=0\) coincides with the spike in \(\eta _{a_0}\) (blue curve, top panel). The cusp where \(\eta _{x_r}\) changes from increasing to decreasing (intersection of yellow and vertical black dashed curves in the bottom panel) also occurs at the liftoff point for \(a_0\)

Full size image

Finally we note an interesting feature in \(\eta ^{\textrm{close}}\): there is an abrupt change in \(\eta _{x_r}\) at the \(a_0=0\) liftoff point (Fig. 7 bottom panel, dashed vertical line). To understand this behavior, note that an instantaneous perturbation of \(x_r\) directly propagates to neural pools through sensory feedback. While all three neural pools are affected by this mechanical perturbation, the neural components of \(\eta ^{\textrm{close}}\) are zero most of the time except when the trajectory lifts off from the \(a_0=0\) constraint (Fig. 7 top panel, blue spike). This observation implies that the system has a high local sensitivity to \(a_0\) during the blue spike, whereas the sensitivity to \(a_1\) and \(a_2\) are significantly smaller than unity at all times. Thus, to understand the effect of perturbing \(x_r\) on the local timing, it is sufficient to focus on \(\eta _{a_0}\) and examine how \(a_0\) reacts to perturbing \(x_r\).

Similar to the forward variational analysis, perturbing \(x_r\) delays the activation of \(a_0\), i.e., \(a_0\) lifts off from \(a_0=0\) at a later time. That is, the displacement in \(a_0\) near the \(a_0=0\) liftoff point is negative. Since \(\eta _{a_0}\) is negative near the liftoff point, perturbing \(x_r\) prolongs the total closed time (i.e., \(\eta _{x_r}\) is positive during the closed phase).

Next we address the cusp phenomena observed in \(\eta _{x_r}\) (Fig. 7, bottom panel, yellow curve). Note that perturbations arriving before the trajectory lifts off from \(a_0=0\) delay the activation of \(a_0\) by increasing the inhibition from its sensory feedback. Moreover, the closer the time of perturbation to the time of liftoff, the larger the delay on the activation of \(a_0\). Such a larger delay leads to a greater increase in the total closed time due to perturbing \(x_r\). Hence, before the liftoff time (Fig. 7, bottom panel, vertical black dashed line), \(\eta _{x_r}\) gradually increases. Once the trajectory has passed the liftoff point, perturbing \(x_r\) delays the activation of \(a_0\) by decreasing its sensory feedback, the effect of which now becomes excitatory. The size of this effect decays exponentially as the trajectory gradually leaves the boundary \(a_0=0\). Thus, there is a cusp in the \(\eta _{x_r}\) curve at the liftoff point, after which \(\eta _{x_r}\) rapidly decreases.

2.4 Robustness to static perturbations

In this section, we show how the robustness of the Aplysia model (1), the ability of the system to maintain its performance despite perturbations, can be quantified using the ISRC, IPRC and LTRC analysis.

Following Lyttle et al. (2017), we quantify the performance or task fitness via the average seaweed intake rate

where \(\Delta x_{r,\varepsilon }\) is the net change in perturbed grasper position \(x_{r,\varepsilon }\) during the grasper-closed phase and \(T_\varepsilon \) is the perturbed period. Note that we assume the seaweed is moving together with the grasper when it is closed and not moving at all during the grasper-open component of the trajectory. Hence, \(-\Delta x_{r,\varepsilon }\) denotes the total amount of seaweed consumed per cycle.

Since the vector field \(F_{\varepsilon }(\textbf{x})\) in system (1) is piecewise smooth in the coordinates \(\textbf{x}\) and smooth in the perturbation \(\varepsilon \), it follows that the following expansion holds:

where \(\Delta x_{r,0}\) is the net change in the unperturbed grasper position during the grasper-closed component of the trajectory. Here, \(\Delta x_{r,1}\) is approximately given by the net change of the \(x_r\) component of the ISRC \(\gamma _1\), which is denoted as \(\gamma _{1,{x_r}}\) (see Sect. 2.2), over the grasper-closed phase. Suppose the grasper closes at \(t^{\textrm{close}}\) and opens at \(t^{\textrm{open}}\) over one cycle. It follows that \(\Delta x_{r,1} =\gamma _{1,{x_r}}(t^{\textrm{open}}) - \gamma _{1,{x_r}}(t^{\textrm{close}})\).

Lyttle et al. (2017) show that the robustness, i.e., the relative shift in the task fitness per relative change in perturbation, for small \(\varepsilon \), can be written as

as \(\varepsilon \rightarrow 0\). Recall that \(T_0\) is the period of the unperturbed limit cycle and \(T_1\) denotes the linear shift in period, which can be estimated using the IPRC (see Sect. 2.3.1).

In summary, the robustness formula can be decomposed into two parts, one involving changes in shape (in particular, the grasper position \(x_r\)) and the other involving the timing change. As discussed before, changes in shape can be estimated using the ISRC and the LTRC analysis, whereas the latter can be quantified using the IPRC. Below, we illustrate the quantification of the robustness by considering the perturbation to be the increase in the constant applied load \(F_{\text {sw}} \rightarrow F_{\text {sw}}+\varepsilon \).

The ISRC with or without timing rescaling corresponding to the perturbation on the applied load have already been computed and discussed in Sects. 2.1 and 2.2. Note that \(\Delta x_{r,1}\) in (5) is the net change in the ISRC during the grasper-closed phase. Choosing the ISRC with rescaling based on the timing of the closing and opening events provides a more accurate estimate of \(\Delta x_{r,1}\). Hence, we use the ISRC with piecewise rescaling to estimate \(\Delta x_{r,1}\), which is the net change in \(\gamma _{1,{x_r}}\) over the closing region per cycle (see the yellow curve over the shaded region in the lower right panel of Fig. 5 and the green arrow marking the difference at the end of the closing phase). Furthermore, the linear shift in the period \(T_1\) can be estimated by (3) using the IPRC.

From the above analysis, we obtain \(\frac{\Delta x_{r,1}}{\Delta x_{r,0}}\approx 0.4806\) and \(\frac{T_1}{T_0}\approx 1.6532\), both of which are positive and are consistent with the concept of an adaptive “stronger and longer” change in the motor pattern in response to increased load. It follows that the robustness is approximately \(-1.1726\times 10^{-2}\). (Note that the smaller this number is in magnitude, the more robust the system is.) To the first order in \(\varepsilon \), the relative change in the performance is then given by

which is illustrated by the red circle as the perturbation size \(\varepsilon \) varies (see Fig. 8, top panel). To see what this means, we take a data point on the line indicated by the arrow, i.e., \((0.42,-0.005)\). Here \(\varepsilon /F_{\textrm{sw}}=0.42\) indicates a \(42\%\) increase in load \(F_{\text {sw}}\), which only causes a \(0.5\%\) decrease in the task fitness, corresponding to a highly robust response. Here the “stronger” effect (i.e., the first term in the robustness formula (5) being positive) contributes to the robustness, whereas the “longer” effect (i.e., the second term in the robustness) reduces it. However, these two effects are not independent from each other: it is the longer retraction-closed time that allows the muscle to build up a stronger force, thereby contributing to a robust response.

Relative change in task fitness \((S(\varepsilon )-S_0)/S_0\) computed numerically (blue stars) versus those obtained analytically from the ISRC and the IPRC according to formula (5) (red circles), as the perturbation \(\varepsilon \) on the seaweed load \(F_{\text {sw}}\) varies. Without perturbation, the nominal applied load is \(F_{\text {sw}}=0.01\). The approximation using the ISRC with different timing rescalings during the grasper-closed (\(\nu _{1,\mathrm close}\)) versus grasper-open phases (\(\nu _{1,\mathrm open}\)) estimated from the LTRC analysis matches the actual simulation (top panel), whereas the ISRC with uniform rescaling \(\nu _1=T_1/T_0\) estimated from the IPRC no longer gives a good approximation (bottom panel)

Full size image

We also compute the relative change in S with respect to \(\varepsilon \) using direct numerical simulations (see Fig. 8, blue stars), which show good agreement with our analytical results. In contrast, if we estimate \(\Delta x_{r,1}\) using the ISRC with a uniform timing rescaling (see Fig. 11), the resulting estimated robustness becomes more negative and no longer gives an accurate approximation to the actual robustness (see Fig. 8, bottom panel). That is, the ISRC using different rescaling factors over the grasper-closed phase (\(\nu _{1,\mathrm close})\) versus the grasper-open phase (\(\nu _{1,\mathrm open}\)), gives a much better approximation to the robustness than the ISRC based on a global timing rescaling \(\nu _1=T_1/T_0\). The fact that the \(\nu _{1,\mathrm {open/close}}\) are obtained via the LTRC analysis highlights the contribution of this novel analytical tool. This observation demonstrates that for systems under certain circumstances [e.g., non-uniform perturbation as considered in system (1)], the ISRC together with the LTRC greatly improves the accuracy of the robustness, compared to the ISRC with global timing analysis given by the IPRC.

2.5 Sensitivity of robustness to other parameters

In general, the performance of motor control systems may be affected not only by external parameters, such as an applied load, but also be internal parameters, for instance describing the physical properties of the biomechanics or neural controllers. The variational tools used in the previous section to understand mechanisms of robustness to increases in applied load—the IPRC, ISRC and LTRC—can also give insights into the effects of changing internal model parameters. For instance, in the SLG model, appropriately varying strengths of protractor or retractor muscles can overcome effects of the increased mechanical load \(F_{\text {sw}}\rightarrow F_{\text {sw}}+\varepsilon \). Because of the SLG model’s relative simplicity, we can relate many of these changes to specific components of the fitness equation in detail.

Below, we first consider how varying sensory feedback strengths can help restore the reduced seaweed intake rate due to increased applied load. Then, we examine how changing the strengths of the protractor and retractor muscles affects robustness to applied loads.

2.5.1 Varying sensory feedback strengths

Figure 9 shows the seaweed intake rate and robustness to the increased load \(F_{\text {sw}}\) with respect to changes in sensory feedback strengths \(\varepsilon _i\), \(i\in \{0,1,2\}\). The performance \(S_0\) becomes negative when \(\varepsilon _0\) or \(\varepsilon _1\) is relatively small (e.g., smaller than \(10^{-5}\)) or when \(\varepsilon _2\) is relatively big (e.g., larger than \(10^{-3}\)), during which the system is in a fast limit cycle/biting mode and hence cannot swallow seaweed.

When the system is in the heteroclinic/swallowing mode, as one might expect, increasing the sensory feedback (e.g., \(\varepsilon _2\)) improves the performance. Surprisingly, our results show that increasing sensory feedback strengths to the two protraction-related neural pools leads to opposite results by decreasing the performance. These results seem to suggest that to restore the deficit caused by the increased load and achieve an increased robustness, we can either increase \(\varepsilon _2\) or decrease \(\varepsilon _0\) and/or \(\varepsilon _1\). However, this is not true. As shown in Fig. 9, a decrease in the robustness can be induced by either decreasing \(\varepsilon _0\) or increasing \(\varepsilon _2\). Moreover, the robustness is largely insensitive to changes in \(\varepsilon _1\), despite the fact that it influences the performance. Understanding these effects on the robustness would require analysis of a second-order variational problem and represents a future direction for understanding neuromodulation.

Effects of varying \(\varepsilon _0\) (top row), \(\varepsilon _1\) (center row) and \(\varepsilon _2\) (bottom row) on the robustness (5) to \(F_{\text {sw}}\) and the unperturbed seaweed intake rate \(S_0\), when \(F_{\text {sw}}=0.01\). Blue curves: Fitness \(S_0\). Black curves: Robustness

Full size image

2.5.2 Varying muscle strengths

Next we investigate how variations of \(k_0\) and \(|k_1|\), the strengths of the protraction and retraction muscles, affect the robustness to changes in seaweed load.

Figure 10 shows that performance improves with the increased protractor muscle strength \(k_0\) or the increased retractor muscle strength \(|k_1|\). This suggests that increasing \(k_0\) or \(|k_1|\) can help restore the deficit in the performance due to the increased mechanical load and hence boost the robustness, which agrees with our numerical simulations (see Fig. 10, top panel, black curve).

Recall that the robustness can be approximated as \(F_\text {sw}\left( \frac{\Delta x_{r,1}}{\Delta x_{r,0}}-\frac{T_1}{T_0}\right) \) [see Eq. (5)]. Understanding the underlying mechanisms of the robustness requires one to investigate how the two quantities involving shifts in shape and timing change with respect to \(k_0\) or \(|k_1|\) (see Fig. 10, lower three panels). We find that increasing \(k_0\) or \(|k_1|\) reduces \(T_1\) and \(-\Delta x_{r,1}\) while \(T_0\) and \(-\Delta x_{r,0}\) are almost unaffected. Hence, both \(\Delta x_{r,1}/\Delta x_{r,0}\) (the “stronger” effect in response to perturbations on the seaweed load) and \(T_1/T_0\) (the “longer” effect ) are decreased as we increase the muscle strengths. However, the reduction in the “stronger” effect is smaller than the reduction in the “longer” effect. As a result, the robustness approximated by \(F_\text {sw}\left( \frac{\Delta x_{r,1}}{\Delta x_{r,0}} - \frac{T_1}{T_0}\right) \) increases as \(k_0\) or \(|k_1|\) increases.

Together, our analytical tools suggest ways in which coordinated changes in intrinsic parameters could maintain fitness and thus enhance robustness.

In this section, we review the classical variational theory for limit cycles (e.g., Filippov 1988; Bernardo et al. 2008; Leine and Nijmeijer 2013; Park et al. 2018), and new tools that we recently developed in Wang et al. (2021) for linear approximation of the effects of small sustained perturbations on the timing and shape of a limit cycle trajectory in both smooth and nonsmooth systems.

Effects of varying muscle strengths \(k_0\) (left panels) and \(|k_1|\) (right panels) on the robustness to \(F_{\text {sw}}\) (top panels, black curve) and the unperturbed seaweed intake rate \(S_0\) (top panels, blue curve). Default parameters \(k_0=1,\,k_1=-1\) represent the strengths and directions of protraction and retraction muscles. The second and third rows of panels show the effects of muscle strengths on timing (\(T_0,T_1\)) and shape (\(-\Delta x_{r,0},-\Delta x_{r,1}\)), respectively. The bottom panels shows how \(T_1/T_0\) (blue) and \(\Delta x_{r,1}/\Delta x_{r,0}\) (red) change as muscle strengths vary

Full size image

In the next two sections we treat the smooth and nonsmooth cases, respectively. In each case, we consider a one-parameter family of n-dimensional dynamical systems

indexed by a parameter \(\varepsilon \) representing a static perturbation of a reference system

3.1 Timing and shape responses to static perturbations in smooth systems

Following Wang et al. (2021), we make the following assumptions:

Assumption 1

• The vector field \(F_{\varepsilon }(\textbf{x}): \Omega \times {\mathcal I}\rightarrow {\mathbb {R}}^n\) is \(C^1\) in both the coordinates \(\textbf{x}\) in some open subset \(\Omega \subset {\mathbb {R}}^n\) and the perturbation \(\varepsilon \in {\mathcal I}\subset {\mathbb {R}}\), where \({\mathcal I}\) is an open neighborhood of zero.

• For \(\varepsilon \in {\mathcal I}\), system (6) has a linearly asymptotically stable limit cycle \(\gamma _\varepsilon (t)\), with finite period \(T_\varepsilon \) depending (at least \(C^1\)) on \(\varepsilon \).

It follows from Assumption 1 that when \(\varepsilon =0\), \(F_0(\textbf{x})\) is \(C^1\) in \(\textbf{x}\in \Omega \) and the unperturbed system (7) exhibits a \(T_0\)-periodic linearly asymptotically stable limit cycle solution \(\gamma _0(t) = \gamma _0(t+T_0)\) with \(0<T_0<\infty \).

Assumption 1 also implies that the following approximations hold:

where \(T_1\) is the linear shift in the limit cycle period in response to the static perturbation of size \(\varepsilon \). This global timing sensitivity, \(T_1\), is strictly positive if increasing \(\varepsilon \) increases the period. The perturbed time \(\tau _\varepsilon (t)\) satisfies the conditions \(\tau _0(t)\equiv t\) and \(\tau _\varepsilon (t+T_0)-\tau _\varepsilon (t)=T_\varepsilon \); it allows approximation (10) to be uniform in time^{Footnote 1} and permits us to compare perturbed and unperturbed trajectories at corresponding time points. The vector function \(\gamma _1(t)\equiv \frac{\partial \gamma _\varepsilon (\tau _\varepsilon (t))}{\partial \varepsilon }|_{\varepsilon \rightarrow 0}\) is the linear (i.e., first-order) shift in the limit cycle shape.

The timing and shape aspects of limit cycles are complementary, and may be studied together by considering the infinitesimal phase response curve (IPRC) and the variational analysis of the limit cycle, respectively.

Infinitesimal Phase Response Curve (IPRC) The IPRC is a classical analytic tool that measures the timing response of an oscillator due to an infinitesimally small perturbation delivered at any given point on the limit cycle. It satisfies the adjoint equation (Schwemmer and Lewis 2012)

with the normalization condition

The linear shift in period \(T_1\) can be calculated using the IPRC as

Forward Variational Equation Classical sensitivity analysis (Wilkins et al. 2009) has been used in many applications to study the shape sensitivity or response of an oscillator to sustained perturbations.

Variational analysis with uniform rescaling. The same perturbation as in Fig. 2 is applied to the Aplysia model (1). Left: The ISRC \(\gamma _1(t)\) with a uniform rescaling over one period. Right: Time series of the difference between the perturbed and unperturbed solutions along \(a_0\)-, \(a_2\)-, and \(x_r\)-directions. The black curve denotes the numerical displacement (\(\Delta y(t) = y_\varepsilon (\tau _\varepsilon (t))-y(t)\)) computed by subtracting the unperturbed solution trajectory from the perturbed trajectory, after globally rescaling time, and aligning trajectories at the onset of closing. The red dashed curve denotes the product of the perturbation \(\varepsilon \) and the ISRC curve. The vertical blue dashed lines indicate the times at which the unperturbed grasper switches from closed to open. Shaded regions and the vertical magenta lines have the same meanings as in Fig. 2. The perturbation is the same as in Fig. 2

Full size image

The dynamics of the linear shift

at time t of the periodic orbit \(\gamma _\varepsilon (t)\) due to a sustained parametric perturbation \(\varepsilon \) initiated at time 0 satisfies the following forward variational equation:

with initial condition \(\textbf{v}(0)\) set by the difference in the perturbed and unperturbed trajectories at the point where they cross the Poincaré section defined by the beginning of the closed phase. Specifically,

Compared with the homogeneous variational equation, which studies the shape sensitivity to instantaneous perturbations, the forward variational equation (14) contains a non-homogeneous term arising directly from the parametric perturbation acting on the vector field.

However, since the perturbed limit cycle has a different period \(T_{\varepsilon }\) and hence a different perturbed time \(\tau _\varepsilon \) due to sustained perturbations, the forward variational equation which neglects such changes in timing fails to give a valid comparison between the perturbed and unperturbed trajectories for times on the order of a full cycle or longer (see Fig. 2C, D). Hence, we adopt a new tool developed in Wang et al. (2021), the infinitesimal shape response curve (ISRC), which incorporates both the shape and timing aspects and captures a more accurate first-order approximation to the change in shape of the limit cycle under a parametric perturbation.

Infinitesimal Shape Response Curve (ISRC) Suppose the rescaled perturbed time can be written as \(\tau _\varepsilon (t)=t/\nu _\varepsilon \in [0, T_\varepsilon ]\) for \(t\in [0, T_0]\). It follows that the relative change in timing denoted by \(\nu _\varepsilon =T_0/T_\varepsilon \) can be represented as \(\nu _\varepsilon =1-\varepsilon \nu _1+ O(\varepsilon ^2)\) where \(\nu _1=\frac{T_1}{T_0}\).

Wang et al. (2021) denote the linear shift in the periodic orbit, \(\gamma _1(t)\) in (10), as the ISRC and adapted Lighthill’s method of “strained coordinates” (Jordan et al. 2007) to show it satisfies the following variational equation. An equation similar to (16) can also be derived by simultaneously Taylor expanding the state variable \(\textbf{x}\) around the limit cycle and its frequency (Keener 2018) (see “Appendix C” for details).

This equation resembles the forward variational equation (14), but has one additional non-homogeneous term arising from time rescaling \(t\rightarrow \tau _\varepsilon (t)\). In contrast to the forward variational dynamics \(\frac{\partial \gamma _\varepsilon (t)}{\partial \varepsilon }\), the ISRC \(\gamma _1(t)\) is periodic with period \(T_0\) (see Fig. 11, left). To see how well the ISRC approximates the actual linear shift between the perturbed and unperturbed trajectories, we plot the linear shift approximated from the ISRC (black curve) and the actual displacement (red dashed curve). Overall, they show good agreement with each other except near the transition between the grasper-closed and grasper-open phases. Such discrepancies arise from the fact that the solution segment at the closing phase has different timing sensitivity to the parametric perturbation compared with the segment at the opening phase, as discussed before. While these small errors are nearly unnoticeable (see Fig. 11, right), they expand when the ISRC result is used to calculate the robustness (see Fig. 8, bottom panel).

Thus, in the case when a parametric perturbation leads to different timing sensitivities in different regions, we use the local timing response curve (LTRC) defined by Wang et al. (2021) to compute shifts in timing in different regions in order to improve the accuracy of the ISRC, as demonstrated when considering perturbations to the load applied to the seaweed (see Fig. 12).

Local Timing Response Curve (LTRC) The accuracy of the ISRC in approximating the linear change in the limit cycle shape evidently depends on its timing sensitivity, that is, the choice of the relative change in frequency \(\nu _1\). In (16), we chose \(\nu _1\) to be the relative change in the full period, by assuming the limit cycle has constant timing sensitivity. It is natural to expect that different choices of \(\nu _1\) will be needed for systems with varying timing sensitivities along the limit cycle. To more accurately capture timing sensitivity of such systems to static perturbations, Wang et al. (2021) defined a local timing response curve (LTRC) which is analogous to the IPRC but measures the linear shift in the time that the trajectory spends within any given region. Specifically, the LTRC is the gradient of the time remaining in a given region until exiting it through some specified Poincaré section - a local timing surface corresponding to the exit boundary of this region. Such a section could be given as a boundary where the dynamics changes between regions, or where a perturbation is applied in one region but not another. For instance, in the feeding system of Aplysia californica (Shaw et al. 2015; Lyttle et al. 2017), the open-closed switching boundary of the grasper defines a local timing surface.

Let \(\eta ^{\textrm{I}}\) denote the LTRC vector for region I. Suppose that at time \(t^{\textrm{in}}\), the trajectory \(\gamma _0(t)\) enters region I upon crossing the surface \(\Sigma ^{\textrm{in}}\) at the point \(\textbf{x}^{\textrm{in}}\); at time \(t^{\textrm{out}}\), \(\gamma _0(t)\) exits region I upon crossing the surface \(\Sigma ^{\textrm{out}}\) at the point \(\textbf{x}^{\textrm{out}}\). Similar to the IPRC, the LTRC \(\eta ^{\textrm{I}}\) satisfies the adjoint equation

together with the boundary condition at the exit point

where \(n^{\textrm{out}}\) is a normal vector of \(\Sigma ^{\textrm{out}}\) at the unperturbed exit point \(\textbf{x}^{\textrm{out}}\). The linear shift in the total time spent in region I, \(T^{\textrm{I}}_1\), is given by

where \(\textbf{x}_\varepsilon ^{\textrm{in}}\) denotes the coordinate of the perturbed entry point into region I. It follows that the relative change in frequency local to region I is given by \(\nu _1^{\textrm{I}}=T^{\textrm{I}}_{1}/({t^{\textrm{out}}-t^{\textrm{in}}})\).

Piecewise uniform ISRC The existence of different timing sensitivities of \(\gamma (t)\) in different regions therefore leads to a piecewise-specified version of the ISRC (16) with period \(T_0\),

where \(\gamma _1^j\), \(F_0^j\), \(F_{\varepsilon }^j\) and \(\nu _1^j\) denote the ISRC, the unperturbed vector field, the perturbed vector field, and the relative change in frequency in region j, respectively, with \(j\in \{\textrm{I, II, III,} \ldots \}\). Note that in a smooth system, \(F_0^j\equiv F_0\) for all j.

Displacements between perturbed and unperturbed trajectories estimated from the ISRC \(\gamma _1\) with piecewise uniform rescaling (red dashed, \(\varepsilon \gamma _1\)) agree well with the actual displacement \(\Delta y(t) = y_\varepsilon (\tau _\varepsilon (t))-y(t)\) where \(y=\{a_0,a_1,a_2,u_0,u_1,x_r\}\) (solid black). Shaded regions, vertical magenta and blue lines have the same meanings as in Fig. 11. The perturbation is the same as in Fig. 2

Full size image

As discussed before, the piecewise-specified ISRC, where \(\nu _1\) takes different values in the closing and opening phases, nicely complements the forward variational analysis. It provides a more self-consistent global description of the shape response of the limit cycle to the mechanical perturbation (see Fig. 5). Displacements between perturbed and unperturbed trajectories estimated using the piecewise-specified ISRC agree well with the actual displacements (see Fig. 12). Moreover, it yields a much better approximation to the robustness compared with the ISRC with uniform rescaling (see Fig. 8).

3.2 Timing and shape responses to static perturbations in nonsmooth systems

As discussed before, system (1) is a piecewise smooth system with one transversal crossing boundary \(\Sigma _{o/c}\) and three hard boundaries (\(\Sigma _0,\Sigma _1,\Sigma _2\)). The study of limit cycle motions in such nonsmooth systems requires analytical tools beyond the standard arsenal of phase response curves and variational analysis, developed for systems with smooth (differentiable) right-hand sides (Spardy et al. 2011a, b; Park et al. 2017). For small instantaneous displacements, variational analysis has been extended to nonsmooth dynamics with both transversal crossing boundaries and hard boundaries for studying the linearized effect on the shape of a trajectory (Filippov 1988; Bernardo et al. 2008; Leine and Nijmeijer 2013; Dieci and Lopez 2011). Analysis in terms of infinitesimal phase response curves (IPRC) has likewise been extended to nonsmooth dynamics for studying the linear shift in the timing of a trajectory following a small perturbation, provided the flow is always transverse to any switching surfaces at which nonsmooth transitions occur (Shirasaka et al. 2017; Park et al. 2018; Chartrand et al. 2018; Wilson 2019). Recently, Wang et al. (2021) extended the IPRC method to nonsmooth systems with hard boundaries.

In nonsmooth systems with degree of smoothness one or higher (i.e., Filippov systems), the right-hand-side changes discontinuously as one or more switching surfaces are crossed. A trajectory reaching a switching surface or boundary has two behaviors: it may cross the boundary transversally or it may slide along it. Hence, there are two types of boundary crossing points: transversal crossing points, at which the trajectory crosses a boundary with finite velocity in the direction normal to the boundary, and non-transversal crossing points including the landing point at which a sliding motion along a switching boundary begins, and the liftoff point at which the sliding terminates. The time evolutions of the solutions to the variational equation (i.e., the forward variational dynamics and the ISRC) and the solutions to the adjoint equation (i.e., the IPRC and the LTRC) may experience discontinuities at a boundary crossing point (Filippov 1988; Bernardo et al. 2008; Leine and Nijmeijer 2013; Park et al. 2018; Wang et al. 2021).

The discontinuity in the variational dynamics when a trajectory meets a boundary crossing point \(\textbf{x}_p\) at crossing time \(t_p\) can be expressed with the saltation matrix \(S_p\) (see Table 1):

where \({\textbf{v}}(t)\) denotes the solution of the forward variational equation or the ISRC, \({\textbf{v}}_p^-=\lim _{t\rightarrow t_p^-} {\textbf{v}}(t)\) and \({\textbf{v}}_p^+=\lim _{t\rightarrow t_p^+} {\textbf{v}}(t)\) represent the solution just before and just after the crossing, respectively.

The discontinuity in \({\textbf{z}}(t)\), the solution to the adjoint equation, at a boundary crossing point \(\textbf{x}_p\) may be expressed with the forward jump matrix (\(J_p\))

where \({\textbf{z}}_p^-=\lim _{t\rightarrow t_p^-} {\textbf{z}}(t)\) and \({\textbf{z}}_p^+=\lim _{t\rightarrow t_p^+} {\textbf{z}}(t)\) are the IPRC or the LTRC just before and just after crossing the switching boundary at time \(t_p\) in forwards time. However, Wang et al. (2021) showed that the jump matrix is not well defined at a liftoff point and hence introduced a time-reversed version of the jump matrix, denoted as \({\mathcal {J}}_p\), defined as follows:

Table 1 summarizes the saltation and jump matrices at different types of boundary crossing points.

3.3 Simulation codes

Simulation codes written in Matlab are available at https://github.com/yangyang-wang/AplysiaModel.

Overview Motor systems are robust—they maintain their performance despite perturbations. Understanding the mechanisms of robustness is important but challenging. To unravel the contributions of different components of robustness, we adopted tools we established in Wang et al. (2021) and reviewed in the methods section (Sect. 3) for studying combined shape and timing responses of both continuous and nonsmooth limit cycle systems under small sustained perturbations. We applied these tools to understand the mechanisms of robustness in a neuromechanical model of triphasic motor patterns in the feeding apparatus of Aplysia developed in Shaw et al. (2015) and Lyttle et al. (2017). We show in the results section (Sect. 2) that this framework lets us analyze how a small sustained perturbation alters the shape and timing of a closed loop system, and thus we began to describe how the neural and biomechanical components interact to contribute to robustness.

The first perturbation we considered was a sustained increase in mechanical load (\(F_{\text {sw}}\rightarrow F_{\text {sw}}+\varepsilon \)). To our surprise, we discovered that long before sensory feedback affected the system, biomechanics played an essential role in robustness by producing an immediate force increase to resist the applied load (Figs. 2, 3 and 4). Furthermore, although the sensory feedback immediately responded to the perturbation, its effect was delayed by the hard boundary properties of the neural firing rates. Our analysis suggests that sensory feedback contributes to the robustness primarily by shifting the timing of neural activation as opposed to changing neuronal firing rate amplitude (Figs. 5, 6 and 7). Our methods can also be readily used to quantify how changes in timing and shape of trajectory affect the robustness (Fig. 8). We find that sensory feedback and biomechanics contribute to the robustness of the system by generating a stronger retractor muscle force build-up during the prolonged retraction-closed phase that resists the increased load. The increased retractor muscle force ultimately leads to more seaweed being consumed during the slightly longer cycle time despite the large opposing forces, thereby contributing to a robust response. These new insights have refined and expanded a previous hypothesis that sensory feedback is the major mechanism that plays a crucial role in creating robust behavior (Lyttle et al. 2017).

Robustness is sensitive to other model parameters. For example, in Sect. 2.5 we investigated how varying internal parameters such as strengths of sensory feedback and muscle activity can help restore the performance that was reduced by an increased applied load (Figs. 9 and 10). Again, we obtained some non-intuitive results. For example, increasing the sensory feedback strength can reduce the robustness rather than improving it (Fig. 9). Moreover, increasing sensory feedback gain has opposite effects on performance and robustness, whereas increasing the protractor or retractor muscle strength improves both performance and robustness. Understanding sensitivities of performance to mixed parameters requires us to go beyond our existing methods. This second-order sensitivity represents an interesting future direction for understanding neuromodulation - the coordinated change of multiple system parameters in order to most effectively counter the effect of an external perturbation (Cropper et al. 2018). There are multiple pathways for neuromodulation, and the simplicity of the model lends itself to detailed analysis of multifactor sensitivities. In future work, we may apply the variational tools used in the present paper for understanding how changes in multiple parameters simultaneously could impact model performance and robustness (cf. Sect. 2.5).

Nonsmooth dynamics and biological realism Our model incorporates two types of nonsmooth dynamics. Both of these features complicate the model analysis, and one might ask whether an “equivalent” smooth formulation might have been employed. We emphasize that both types of nonsmooth dynamics provide better reflection of the underlying biophysics than a “smoothed” version would do, and contribute in fundamental ways to the biological mechanisms we study. Our model assumes that neurons fire, once excited, at a nonzero rate, and maintain a rate of exactly zero (rather than “very small”) when inhibited. Thus our motor pool variables \(a_i\) have hard boundaries at zero firing rate. It is well known that the Hodgkin–Huxley model, for example, fires at very high rates when provided even relatively small currents, and slower rates are only possible with additional ionic conductances, such as the A current (see, for example, Hille (2001)). Studies on the energetics of neurons in real brains do not assume that they can fire at values much lower than a few Hertz. See, for example, Figure 2 in Laughlin and Sejnowski (2003), where the minimum firing rate for rat cortex is estimated to be about 3 Hz. Thus, the model we are using Lyttle et al. (2017) is more realistic than one that assumes that neurons can fire infinitely slowly [e.g., a sigmoid function or a hyperbolic tangent function (Ermentrout 1998)]. Thus, assuming that neurons fire once excited at a nonzero rate (and thus have a “hard boundary” at 0 firing frequency) is more biologically realistic. Moreover, we have attempted to replicate the results in this paper using an alternative formulation in which we replace the hard boundary with a “soft” boundary, implemented using a sigmoidal firing rate function (Harris and Ermentrout 2015) and found that eliminating the hard boundary drives the system from the so-called heteroclinic cycling regime (Shaw et al. 2015; Lyttle et al. 2017) to the “limit cycle regime” in which it fails to consume seaweed at a rate sufficient to support survival. In addition to the hard “sliding boundary” at zero firing rate Wang et al. (2021) our model has a nonsmooth transverse crossing of a Poincaré section at the point when the grasper transitions from “open” to “closed,” and the biomechanics switch from being free of the mechanical loading to engaging the mechanical load of the seaweed. It is important to realize that when an animal encounters a load, that interaction creates a nonlinear change in the dynamics of the system. Although some investigators attempt to finesse this aspect by trying to add a smoothly changing load function (such as a spring with a very stiff spring constant), an entire field devoted to hybrid systems has developed to study how dynamics evolves when there is a discontinuous change due to interaction with the environment [see, for example, Holmes et al. (2006) and Aihara and Suzuki (2010)].

Experimentally testable predictions The surprising result that the length-tension curves of the opposing muscles generate an instantaneous response to force perturbations could be tested, at least initially, using some of the more realistic biomechanics models that have been developed of Aplysia feeding.

For example, in a detailed kinetic model that does not have sensory feedback (Sutton et al. 2004; Novakovic et al. 2006), one could apply a step increase in force when the odontophore is closed and the retractor muscle is activated while measuring the force resistance to that change, and compare that to a purely passive response in which the retractor muscle is not activated. The results of this paper predict that there will be significant differences between these conditions.

In a model that does have sensory feedback (Webster-Wood et al. 2020), one could apply a step increase in force when the odontophore is closed and measure the change in force and the duration of the cycle to determine how that perturbation alters fitness. This paper’s results predict that the response to a sustained perturbation will be smaller in the presence of sensory feedback and will be larger if sensory feedback is removed.

The model suggests that there may be delays from the time that sensory feedback is available to the time that force changes. Using the model, sinusoidal force changes could be applied at different frequencies to determine the predicted phase lag, and this effect could be tested in the real animal.

Results shown in Fig. 9 suggest that the model is relatively insensitive to changes in the strength of sensory feedback over a wide range of gains. Thus, one experimental test might be to increase or decrease the strength of sensory feedback to show that robustness to changing mechanical loads is not significantly affected. One way to test this hypothesis would be to use the newly developed technology of carbon fiber electrode arrays, which could be used to excite, inhibit, and record from many sensory neurons simultaneously (Huan et al. 2021).

In contrast, results shown in Fig. 10 suggest that changing the relative strengths of the muscles can have larger effects on robustness. Previous studies have shown that neuromodulators can speed up and strengthen muscular contractions and thus might contribute to robustness (Taghert and Nitabach 2012; Lu et al. 2015; Cropper et al. 2018). Studies of the neuromuscular transform (Brezina et al. 2000) suggested that neuromodulation could effectively speed up and strengthen feeding responses in normal animals, and thus might contribute to robustness.

Future experimental studies could be guided by coordinated changes of parameters in this model using the analysis tools we have presented.

Caveats and limitations Tracking possible transitions into and out of constraint surfaces becomes combinatorially complex as the number of distinct constraint surfaces grows. Here we impose three hard boundaries at \(a_i\ge 0\), as discussed above, by requiring firing rates to be nonnegative. An earlier model specification given in Shaw et al. (2015) and Lyttle et al. (2017) also required firing rates to be bounded via the constraint \(a_i\le 1\). Here we relax this constraint for computational convenience, since the coexistence of multiple constraints requires encoding entry/exit conditions and vector field restrictions for all feasible combinations of constraints. In practice, comparison of simulations with and without the \(a_i\le 1\) constraint give qualitatively and quantitatively indistinguishable results under most conditions.

Our analysis is in principle limited to small perturbations. Large perturbations lead to crossing of bifurcation boundaries in which the behavior switches to a different dynamical mode. “Robustness” in a broader sense can mean the distance to a basin of attraction of another dynamical attractor. For example, if the force is increased too much, the model will collapse into a stable fixed point with overextended protraction, while the animal will engage a different response to release or sever the seaweed to avoid damage to its feeding apparatus. This aspect is not captured in the variational approach. Nonlinear and bifurcation analysis could complement the present study and is ripe for investigation in future work.

In this paper, we considered a specific perturbation, namely increasing the force opposing seaweed ingestion \(F_\text {sw}\rightarrow F_\text {sw}+\varepsilon \). Note that in this formulation, the perturbation parameter \(\varepsilon \) carries the same units (force) as \(F_\text {sw}\). Consequently, in order to use a unitless measure of robustness, expression (5) includes a factor of \(F_\text {sw}/\varepsilon \). Also, in this formulation, the timing sensitivity \(T_1\) (shift in period per increase in force) and shape sensitivity \(\gamma _1\) (shift in limit cycle shape per increase in force) have units including reciprocal force. As an alternative formulation, which might facilitate comparison of robustness to perturbations across different modalities, one could rewrite the force perturbation as \(F_\text {sw}\rightarrow F_\text {sw}(1+\varepsilon )\). In this case \(\varepsilon \) would represent a unitless measure of relative perturbation size. The subsequent variational, IPRC, ISRC and LTRC analysis would remain unchanged, except the resulting quantities Z, \(T_1\), \(\gamma _1\), and \(\eta _1\) would undergo a change in units, hence a multiplicative (fixed) change in scale. An advantage of specifying perturbations as a relative or unitless quantity would be that a similar analysis to that undertaken in this paper could be applied to other modalities in the same system or across disparate systems.

Generalizability to other systems Although we focused in the present work on the robustness of the mean rate of seaweed intake with respect to increases in the force opposing ingestion, our analysis carries over to other objective functions (e.g., calories consumed per energy expenditure) as well as other perturbations (e.g., temperature, which may alter the speed of feeding in Aplysia). The variational approach to analyzing robustness should apply to any reasonable (e.g., smoothly differentiable) objective function and any parameter represented in the system, e.g., adjustments to changes in speed, steepness, or right-left asymmetry of walking movements on a (split) treadmill system (Frigon et al. 2013; Embry et al. 2018).

Our methods might also provide insights into how rapidly a system can adjust to small modulation of forces. One could possibly conduct experiments to study the linear response of a system to modulation of applied force \(F_{\textrm{sw}}\), such as an instantaneous small change from one static force to another, or a small amplitude sinusoidal modulation of \(F_{\textrm{sw}}\). The infinitesimal shape response curve and other variational tools developed in Wang et al. (2021) might play a role in the linear response analysis. This treatment could represent an interesting future direction.

The present manuscript applies variational methods to understand the robustness in a specific Aplysia neuromechanical model (Lyttle et al. 2017). This model makes significant simplifications to the real feeding apparatus control system in order to gain mathematical tractability and analytical and biological insights. Nonetheless, the framework developed in Wang et al. (2021) applies naturally to more elaborate dynamical models of Aplysia feeding such as Webster-Wood et al. (2020) and models incorporating conductance-based network descriptions of the central pattern generator (Cataldo et al. 2006; Costa et al. 2020). Thus, what we have done here provides a framework for understanding neural control of motor behaviors like the one considered in this paper.

More broadly, motor control beyond the Aplysia feeding system is also amenable to the analysis of the sort developed in Sect. 3 (Wang et al. 2021). For example, the stability of bipedal walking movements remains a challenge in the field of mobile robotics Vukobratovic et al. (2012); Westervelt et al. (2018). Biologically inspired robotics continues to provide alternative approaches with greater robustness than conventional devices (Beer 2009; Pfeifer et al. 2007; Beer et al. 1997; Goldsmith et al. 2019). The variational framework exhibited here applies to these systems as well (Fitzpatrick et al. 2020). In the context of any motor control model, the variational analysis we present here should allow analysis of robustness of any reasonable objective function with respect to any system parameter.