{'fraction_non_alphanumeric': 0.033977888287008255, 'fraction_numerical': 0.015862268004617708, 'mean_word_length': 4.510367606581442, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 7, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Control in triphasic rhythmic neural systems: a comparative mechanistic analysis via infinitesimal local timing response curves

Similar activity patterns may arise from model neural networks with distinct coupling properties and individual unit dynamics. These similar patterns may, however, respond differently to parameter variations and, specifically, to tuning of inputs that represent control signals. In this work, we analyze the responses resulting from modulation of a localized input in each of three classes of model neural networks that have been recognized in the literature for their capacity to produce robust three-phase rhythms: coupled fast-slow oscillators, near-heteroclinic oscillators, and threshold-linear networks. Triphasic rhythms, in which each phase consists of a prolonged activation of a corresponding subgroup of neurons followed by a fast transition to another phase, represent a fundamental activity pattern observed across a range of central pattern generators underlying behaviors critical to survival, including respiration, locomotion, and feeding. To perform our analysis, we extend the recently developed local timing response curve (lTRC), which allows us to characterize the timing effects due to perturbations, and we complement our lTRC approach with model-specific dynamical systems analysis. Interestingly, we observe disparate effects of similar perturbations across distinct model classes. Thus, this work provides an analytical framework for studying control of oscillations in nonlinear dynamical systems, and may help guide model selection in future efforts to study systems exhibiting triphasic rhythmic activity.


Introduction
Many of the automatic behaviors fundamental to animal life, including respiration, digestion, and various forms of locomotion, represent multi-phase processes, in each of which a sequence of motor unit activations repeats rhythmically. In each case, the behavior is driven by the generation of sustained, multi-phase rhythmic outputs by a corresponding neural circuit that, for automaticity, must be able to produce its characteristic activation pattern in the absence of temporally-varying inputs. Nonetheless, inputs do serve important functions for these neuronal activity patterns; transient inputs may serve as triggers to turn rhythms and behaviors on and off, they may allow for interruption of rhythms (e.g., pausing of inspiration to allow for swallowing or vocalization) and switching between rhythms (e.g., transitions in swimming and walking in salamanders (Ijspeert et al., 2007)), and they may serve as control signals to modulate the frequency or relative phase durations of ongoing rhythms.
The ubiquity of these rhythms, and the relative experimental accessibility of the central pattern generator (CPG) circuits that often produce them, have made them a common subject of scientific investigation. In the computational realm, past modeling works have presented several distinct mathematical frameworks that can generate repetitive, multi-phase rhythms. In each case, the rhythms can exhibit extended activations of specific units in a circuit, occurring sequentially in each cycle. Despite its functional importance, however, the topic of control of such circuits has received far less theoretical attention. On the one hand, control of linear systems is a thoroughly developed topic (cf. Brockett (2015)). Fixed-point control for nonlinear systems is understood to the extent that one can approximate the nonlinear system under weak disturbances with a related linearized system (cf. Isidori (1985)). Despite decades of development in these areas, however, control for nonlinear dynamical systems exhibiting limit cycle dynamics remains relative underdeveloped. One reason for this lack of development may be that parametric control of limit cycles requires understanding the coordinated changes in both the shape and timing of trajectories upon adjustment of a control parameter. For smooth systems, Floquet analysis can capture the effects of small perturbations, but in many naturally occurring control systems -for instance in biological motor control systems -there are switching surfaces or fast-slow dynamics that require more subtle treatment. Such nonsmooth dynamics can occur at the point of making or breaking contact between the organism's body and an external substrate, or at hard boundaries for internal dynamics such as exclusion of negative firing rates in neural network CPG models formulated in the firing rate framework.
In this paper, we present a novel analysis of control in rhythm-generating neuronal circuits, achieved by applying and extending recently developed tools that allow linearized (sometimes called variational) analysis of the parametric variations induced in periodic orbits by changes in constant input levels. As noted in previous work (Yakovenko et al., 2005;Olypher et al., 2006;Daun et al., 2009;Rubin et al., 2009), the functional control of CPGs requires not only the capacity to alter the overall output frequency but also, in general, a mechanism to tune the durations of specific phases within a rhythmic pattern, such as prolonging an inspiration without altering expiration or prolonging the stance phase of locomotion without altering the swing phase; in practice, however, changes in the duration of activity at one node in a coupled circuit may induce cascading changes in other nodes. Therefore, in this paper we study how changing the input intensity to a single unit within a model CPG circuit impacts the durations of all phases within its rhythmic output cycle, with an eye towards the independent phase modulation property.
An important aspect of the work presented here is that we consider this control problem within several classes of mathematical models that naturally produce multi-phase oscillations and have appeared in the neuroscience literature. One such modeling framework consists of relaxation oscillators coupled with mutually inhibitory synapses. The individual relaxation oscillator model includes dynamics on distinct fast and slow timescales and produces periodic dynamics featuring alternating, extended active and silent phases, with fast transitions between them. This pattern provides a strong qualitative resemblance to the outputs of units within CPGs, and the fast-slow nature of the dynamics involved has proved convenient for a broad collection of analyses (cf. Chapter 9 of Ermentrout and Terman (2010) as well as Bertram and Rubin (2017) for relevant methods). In particular, it is natural to represent such a CPG model as a collection of units, each with its own active and quiescient states, and to focus on solutions in which at most one unit is in its active state at any given time. On the other hand, it is not necessarily the case that all components of a CPG circuit have the capacity to oscillate individually in some input range. One alternative to the coupled relaxation os-cillator CPG model is provided by models featuring heteroclinic cycles, sometimes referred to as winnerless competition models (May and Leonard, 1975;Afraimovich et al., 2004). In such models, an extended passage near a saddle point yields the prolonged duration associated with each phase of a CPG output. Here, we no longer have a fast-slow timescale splitting and decomposition of the phase space into projections for individual units. A third modeling framework for multi-phase neuronal rhythms, which also lacks a timescale splitting, is offered by the class of competitive, threshold-linear networks (Morrison et al., 2016). Inhibitory interactions between threshold-linear units can yield oscillations without stable fixed points, as is suitable for functional CPG dynamics, and this modeling framework has been proposed as being convenient for mathematical tractability.
In each of these modeling frameworks, we ask how the architecture influences the controllability of the pattern of activation. Specifically, we analyze the sensitivity of the individual activation phases comprising the overall pattern of activity in response to the tuning of an excitatory, tonic driving current external tothat is, not generated by or affected by -the activity of the circuit. In all cases, we use a recently developed mathematical tool, the local timing response curve (lTRC) (Wang et al., 2021). The lTRC provides a firstorder (linear) approximation for the effect of a parametric perturbation on the time it takes for a trajectory to pass through a given region. It decomposes the effect into a component affecting the trajectory's entry position into the region, an impact on the vector field within the region, and an influence on its exit from the region. Importantly for CPG models, this idea applies to vector fields that lack smoothness, or even feature discontinuities, at switching surfaces within the phase space. The lTRC analysis provides a mathematically grounded numerical quantification of the effects of a specific parameter variation on rhythm period and phase durations, which we can represent in terms of entrance, within-region, and exit effects. In several cases, we supplement the lTRC calculation with mathematical analysis based on timescale decomposition or other ideas, to obtain additional insights about the dynamical mechanisms underlying the observed effects. Our findings yield important information for future work modeling neuronal generation of multi-phase rhythms, in that they provide insights about the control properties of each framework that we consider, which can guide the selection of a modeling framework to match experimentally observed circuit properties.
In all cases, for concreteness, we consider three-phase, or triphasic, rhythms. This choice is motivated by the common appearance of triphasic rhythms in CPG outputs, including those for mammalian respiration, tripod insect gaits, escape swimming in mollusks, scratching in turtles, feeding in sea slugs, and digestion in crustaceans (Marder et al., 2007;Smith et al., 2007;Büschges et al., 2008;Daun and Büschges, 2011;Hao et al., 2011;Wojcik et al., 2014;Shaw et al., 2015). We organize the rest of our paper as follows. A brief description of the modeling frameworks we consider and the related phase-transition mechanisms is given in section 2. Section 3 presents the main analysis and results on the controllability of phase durations in each phase transition mechanism. Last, we summarize the mapping that we obtain between model classes and control properties in section 4, where we also discuss limitations, connections to previous literature, and possible implications of our results for biology and for future work.
2 Model systems 2.1 Relaxation oscillator circuit A vast array of planar dynamical systems, including many neuronal models, can produce relaxation oscillations. In the neuronal case, the key ingredients are typically a fast inward current and a slow negative feedback, which can be the inactivation of the inward current or the activation of an outward current. From a mathematical perspective, all such models that produce a transition from a stable hyperpolarized state to oscillations via an initial Andronov-Hopf (AH) bifurcation and a subsequent transition to a stable depolarized state via another AH bifurcation, as a single input current parameter is varied, are equivalent (up to criticality of the AH bifurcations) (Izhikevich, 2007). When such units are coupled together synaptically and produce sequential oscillations, however, their responses to parameter tuning can depend on the nature of the mathematical mechanism by which the neuronal activations alternate (Daun et al., 2009;Rubin et al., 2009). Thus, in this work, we consider a specific choice of neuronal relaxation oscillator model and we vary parameters within the context of that model to realize and compare several different transition mechanisms.
Equations For our individual neuron model, we follow past works on CPGs that used a slowly inactivating persistent sodium current (I NaP ), as has been characterized in respiratory CPG neurons (Butera Jr et al., 1999a,b), to support oscillations (e.g., Rybak et al. (2006); Daun et al. (2009); Rubin and Terman (2012)). The full circuit model includes mutual inhibitory coupling among three units. We make the assumption that the coupling strength responds instantaneously to changes in voltage, which simplifies the analysis and represents a reasonable approximation because the synaptic conductance rise and decay rates are much faster than the timescale of I NaP (de)inactivation and which has been successfully used in previous analysis (see Daun et al. (2009); Rubin and Terman (2012)).
The specific model equations that we consider are modified from a three-component respiratory CPG model (Rubin and Terman, 2012): (1) In system (1), F describes a persistent sodium current and a leak current, The coupling function S ∞ is given by a monotone increasing function taking values in [0, 1], where θ I represents the synaptic threshold. This function closely approximates a Heaviside step function when σ I is small, and it is multiplied by a coupling strength b ij for i, j ∈ {1, 2, 3} each time it appears. The final term, g E d i (v i − V E ), in each voltage equation represents a tonic excitatory drive, and we consider the strength factors d i as control parameters. In the equations for gating variables, τ h is the time constant for the sodium conductance inactivation, given by where 0 < 1. Additional details about the functions in system (1), as well as parameter values used for simulations, are given in Appendix A. Figure 1 shows a typical solution of system (1), where each cell undergoes relaxation-oscillator-like dynamics. Because the coupling between units in (1) only occurs through the functions S ∞ (v i ), which transition rapidly between values near 0 and 1, we can represent the solution as a collection of three phase plane projections, one to each coordinate pair (v i , h i ). In each phase plane, since is small, the projected trajectory will evolve along a branch of the v i -nullcline, except for rapid jumps between branches at the fold points where pairs of branches meet and end. A rapid change in an S ∞ term in a neuron's voltage equation yields a rapid change in its v-nullcline position and corresponding jump to a new branch. An example of the (v 1 , h 1 ) projection is shown in Figure 1 In the solutions we consider, there is always exactly one cell active at each time (except during jumps), which inhibits the other two cells. As parameters are varied, a number of different activation patterns can be generated, e.g. (132), (1232), (13123), etc. (see more patterns in Rubin and Terman (2012)). In the rest of the paper, we consider the symmetric (123) firing pattern of the system, also known as splay states, in which the three cells alternate activating in the activation order 1 → 2 → 3 → 1 → 2 → 3 → · · · . Because all system parameters are the same for all neurons in the simulation shown in Figure 1, a splay state occurs. Thus, the phase plane projections for the other two units are identical to those shown for (v 1 , h 1 ); although the solution time courses are phase-shifted (Figure 1 (a)), that difference does not show up in the phase plane ( Figure 1(b)), since time is not explicitly represented there. Skinner et al. (1994) recognized that in pairs of half-center oscillators with reciprocal inhibitory coupling, alternations of which cell is active can occur through four different transition mechanisms: intrinsic release, synaptic release, intrinsic escape and synaptic escape (see also Wang and Rinzel (1992); Zhang and Lewis (2013); Yu and Thomas (2021)). These concepts extend naturally to larger networks with reciprocally inhibitory coupling Sakurai and Katz, 2022). The intrinsic properties of the neurons and the properties of the synapses between them, especially the position of the synaptic threshold, determine the fundamental differences of the four mechanisms. Figure 1 illustrates the dynamics of the intrinsic release mechanism. During the network oscillation, the active cell reaches the lower right knee of its active (uninhibited) v-nullcline and jumps down. As it does so, its voltage crosses through the synaptic threshold θ I , releasing the suppressed cell with the largest h-value from inhibition. Although the other suppressed cell's voltage also transiently increases, it is quickly suppressed by the activation of the neuron with the largest h (see grey arrows in Figure 1). The formerly active cell becomes inhibited by the newly-active cell and thus completes its jump down on the left branch of its inhibited v-nullcline.

Intrinsic release
Intrinsic escape In the intrinsic escape mechanism, the oscillations are controlled by the suppressed cell escaping from inhibition when it reaches the upper left knee of its v-nullcline. Figure 2 shows the time (v 1 , h 1 ) phase plane of the solution. Black solid curve: cell 1 limit cycle trajectory; green solid curve: h 1 -nullcline; blue dashed curve: inhibited v 1 -nullcline; blue solid curve: active v 1 -nullcline; magenta dashed vertical line: synaptic activation threshold θ I . Phase planes for (v 2 , h 2 ) and (v 3 , h 3 ) are the same as (v 1 , h 1 ). Note the transient upward voltage deflections (grey arrows) during the brief interval in which the silent cell is released from inhibition but then suppressed due to the faster activation of the other silent cell (cf. Rubin and Terman (2012); Park and Rubin (2022)).
course of a splay state with all transitions by intrinsic escape and the phase plane projection of the trajectory to (v 1 , h 1 )-space. At first glance, the projection may look identical to the intrinsic release case. Note, however, that cell 1 reaches the left knee of its inhibited v-nullcline (dashed blue) and escapes, or jumps up, from there despite still being inhibited. Moreover, once cell 1 is active, it converges to a stable critical point where its active v-nullcline (solid blue) intersects its h-nullcline (green) at v > θ I . The subsequent jumpdown of cell 1 happens only because cell 2 manages to reach its left knee and escape. Different from the intrinsic release in which a transient increase occurs for the voltage of the suppressed cell, in the intrinsic escape case a transient decrease exists due to double inhibition from the other two active cells (see grey arrows in Figure 2), which we call a "double inhibition notch".
Synaptic release The synaptic mechanisms come into operation when the synaptic threshold θ I intersects one of the outer branches of the voltage nullclines, rather than the middle branch. In this case, the cells' transitions between the active and suppressed states occur due to synaptic dynamics. Figure 3 illustrates an example of the synaptic release mechanism, where θ I lies on the right branch of the active v-nullcline. The jumps occur when the active cell reaches the synaptic threshold voltage. When this happens, it no longer fully inhibits the suppressed cells, and the one with largest h can activate. When that cell's voltage jumps up, it crosses the synaptic threshold and inhibits the active cell, which in turn jumps down, starting from v ≈ θ I since all of this action happens on the fast timescale. Note in Figure 3(b) that the resulting trajectory projection does not pass through either the local maximum or minimum of the v-nullcline.
Synaptic escape In this mechanism, the synaptic threshold intersects the left branch of the inhibited nullcline, as shown in Figure 4. The jumps are determined by when an inhibited cell's voltage crosses through the synaptic threshold. When this occurs, the active cell becomes inhibited and jumps down, which turns off the inhibition to the suppressed cells and allows the inhibited cell with the newly supra-threshold voltage to activate. Shaw et al. (2015) and Lyttle et al. (2017) introduced a firing-rate model for the feeding pattern generator of the marine mollusk Aplysia californica, adapted from May and Leonard's three-population heteroclinic cycling model (May and Leonard, 1975). Here we study a three-dimensional piecewise linear system adapted from Park et al. (2018) that was introduced as a piecewise linear simplification of the (smooth) heteroclinic cycling model. This system takes the form

Heteroclinic cycling model
(2) Figure 3: Dynamics of the synaptic release mechanism. Colors as in Figure 1. The synaptic threshold θ I intersects the right branch of the active voltage nullcline. The release occurs when the active cell's voltage declines to the synaptic threshold.  (b) When a 1 = a 2 = a 3 = 0.01, the system contains a stable limit cycle. Note that in this case, the critical points are no longer exactly at the corners of the unit cube.
This system incorporates firing rates of three neural pools -the "protraction-open" pool x, the "protractionclosed" pool y, and the "retraction" pool z. Parameter a i represents intrinsic neural excitation to each pool, and ρ is a coupling constant representing inhibition between the neural pools. The domain of system (2) is divided into three equal rectangular pyramids that together comprise the unit cube, as shown in Figure 5. When all a i = 0, the intrinsic neural dynamics contain a stable heteroclinic cycle connecting saddles at (1, 0, 0), (0, 1, 0), and (0, 0, 1), see Figure 5(a). When some a i are small and positive, the heteroclinic cycle is broken and a stable limit cycle arises, see Figure 5(b).

Competitive threshold-linear network
The third framework we consider is a competitive threshold-linear three-unit network governed by the following ODE system (Morrison et al., 2016): The connectivity matrix W = [W ij ] in (3) gives directed connection strengths between pairs of nodes, the parameter θ i is an external drive to node i, and the nonlinearity [·] + is given by [y] + = max{y, 0}.
Assume the connectivity matrix W to be where δ > 0 and 0 < < 1. It is shown in Morrison et al. (2016) that in a unidirectional coupled network, if < δ 1 + δ , then the network has bounded activity and no stable fixed points. Specifically, to have the nodes oscillate in the cyclic order (123), we can take where = 0.25 and δ = 0.5 to satisfy the given inequality.
If all of the drive terms θ i are set to 0, then the system trajectories converge to the stable fixed point (0, 0, 0). When the θ i are all positive, the solution trajectory always converges to a stable limit cycle, which surrounds an unstable fixed point. This oscillation mechanism differs from that in the heteroclinic cycling model ( §2.2), where the periodic orbit spends its time in small neighborhoods of invariant manifolds associated with equilibria.
To make a clear comparison with the heteroclinic cycling model (2), we again partition the unit cube domain into three equal pyramids and define the active region for each node to be

Phase-duration sensitivity
A functional CPG should produce a periodic oscillation that is stable to perturbations and can be tuned in response to environmental conditions and metabolic demands. A natural expectation might be that each component of a CPG circuit would receive an independent control signal, which could be adjusted to modulate its active phase duration and amplitude. This view, however, overlooks the reality that when a dynamical system produces a periodic behavior, the coupling among the system components may interfere with or prevent such distributed control (e.g., Daun et al. (2009)). Thus, although there may be an independent input to each network component, as is the case with the models presented in the previous section, tuning of these inputs may affect the properties of multiple phases of the output rhythm.
Therefore, in this section, we study the responses of periodic solutions in each model class to variations in a control parameter. Since we consider the case of symmetric networks, in that the dynamics of all network components and the values of all coupling strengths are identical, we (without loss of generality) vary the control parameter to component one in each model and study the effect of this variation on all three phase durations within the perturbed oscillation. We then explain the observed responses via two approaches: a quantitative analysis using the local timing response curve (lTRC) and various model-specific qualitative analyses based on other mathematical tools.
The lTRC calculation developed in Wang et al. (2021) quantifies the timing sensitivity of a limit cycle within any given local region bounded between specific Poincaré sections in phase space, by decomposing the response to a static perturbation into components that relate to entry section crossing, exit section crossing, and alteration of the vector field within the region. Consider an ordinary differential equation model for x(t) in a phase space Ω ⊂ R n , n ≥ 2, and a parameter µ ∈ R. Without loss of generality, we take µ = 0 as a baseline or unperturbed case and let F (x(t)) denote F (x(t), 0). Suppose that in the baseline case the system supports a stable periodic solution, or limit cycle, γ. Suppose that γ passes through a region or subdomain R ⊂ Ω, with entry point x in and exit point x out . For each x ∈ γ ∩ R, let Γ(x) be the time it takes for the trajectory that originates at x to exit from R. The lTRC is defined to be the gradient of Γ(x), η(t) = ∇Γ(x(t)), which satisfies the adjoint equation where DF denotes the Jacobian of the vector field F . At the exit point x out of region R, η satisfies the boundary condition where n out is a normal vector of the exit boundary surface at x out (see more details about the lTRC in Appendix B.1).
Using the lTRC, Wang et al. (2021) provide a formula for analytically calculating the first-order or linear approximation for the change in the duration a trajectory spends in region R, induced by a small change in µ. That is, writing Γ µ (x) for the time to exit R from point x under the perturbation µ, we have Γ µ (x) = Γ(x) + µT 1 (x) + O(µ 2 ), with: where x in , x in µ represent the unperturbed entry point and the perturbed entry point to region R, respectively, and F µ (x(t)) := F (x(t), µ). The first-order duration shift T 1 given by (4) consists of two terms: the first term arises from the impact of the perturbation on the entry point to the region; the integral term shows the impact of the perturbation on the vector field.
Note that in Wang et al. (2021), boundary surfaces of region R were assumed to be fixed in the derivation of (4). When considering boundary surfaces that may themselves change position as control parameters vary (e.g., the boundaries of each region in the heteroclinic cycling model (2) shift as parameters a i are varied), a third term for T 1 arises, accounting for the impact of the perturbation on the exit point: where δ(·) is the Dirac delta function. When the exit boundary is fixed, η(x out ) = ∇Γ(x out ) is parallel to n out , the vector normal to the exit boundary. In this case, n out ⊥ ∂x out µ /∂µ; therefore, term B is identically zero and (5) naturally reduces to (4). The derivation of our generalized formula (5) for the shifting-boundary case is given in Appendix B.3.
With the linear duration shift T 1 calculated from equation (4) or equation (5), we can obtain the phase duration change to the first order, given by µ × T 1 . Note that T 1 > 0 (resp. T 1 < 0) if a positive perturbation increases (resp. decreases) the phase duration. In the following subsections, we use the lTRC calculation to study the phase duration sensitivity to a specific perturbation (with a comparison to direct numerical simulation), and we supplement the lTRC analysis with other mathematical approaches such as fast-slow decomposition and vector field analysis.

Relaxation oscillator circuit
In the relaxation oscillator circuit (1), we start from the symmetric case (d 1 = d 2 = d 3 = 1) and apply a small sustained perturbation µ to d 1 (i.e., d 1 → d 1 + µ). We immediately run into a complication in this case: past work has shown that the response of such models to tuning a constant input parameter depends on the phase transition mechanism operating in the unperturbed periodic orbit (Daun et al., 2009;Rubin et al., 2009). We consider the four transition mechanisms discussed by Skinner et al. (1994) and reviewed in the previous section: intrinsic and synaptic release and escape. We display and compare the sensitivities of the active phase durations for three of the four cases; to our surprise, we find that, unlike in a two-cell network, the symmetric periodic orbit with transitions by intrinsic escape is unstable to small perturbations Figure 7: Projections to the phase plane for cell 1 in the intrinsic release case for d 1 = 0.7 (blue), d 1 = 1 (black) and d 1 = 1.3 (red). The color of each trajectory is in line with the color of the corresponding v 1 nullcline; the h 1 -nullcline appears in green. Solid nullclines are relevant when cells 2 and 3 are inactive, while dashed nullclines arise when cell 2 or cell 3 is active. Successive dots are equally spaced in time. As d 1 increases, the v 1 nullclines are lowered, and consequently, the jump-down position (right knee) of cell 1 is lowered (see bottom rectangle). The trajectories' proximity to the h 1 -nullcline translates a small spatial dispersion into a large temporal effect, such that the prolonging of the active phase at jump-down is much more significant than the reduction occurring in the entry to the active phase (upper rectangle), far from the h 1 -nullcline.
in the three-cell network, so we conclude this section with commentary on this instability (with additional analysis of the instability in Appendix C).
The dominant effect of increasing the value of d 1 is that the active v 1 -nullcline, along which cell 1 moves during the active phase, is lowered, as illustrated in Figures 7 and 8. For visual convenience, Figures 7,8 show the trajectories with a perturbation of magnitude |µ| = 0.3 applied to d 1 = 1, but in our actual analysis and lTRC calculation we impose a smaller perturbation of magnitude |µ| = 0.05.
We define neuron i as being in the active phase when v i > θ I , the synaptic threshold. The sections {v i = θ I } partition phase space into 8 regions, based on the various possible combinations of active and inactive neurons. For each neuron, we are interested in its active phase duration, so for each i ∈ {1, 2, 3} we define the Poincaré sections and we quantify the time from the passage of the trajectory through P i,up until its passage through P i,down for each i. By equation (5), for each region we can analytically calculate the linear shift in the duration of each active phase when a small static perturbation µ is applied to d 1 .
Intrinsic release Table 1 presents the active phase duration changes for a solution featuring the intrinsic release transition mechanism, obtained both by direct numerical simulation and by the lTRC calculation. , the jump-down occurs with the v 1 -coordinate invariant with respect to d 1 at v 1 ≈ θ I , which is farther from the h 1 -nullcline, such that the change in spatial position does not translate into as large of an increase in the active phase duration. Table 1: Intrinsic release: duration changes of the three individual active phases, computed by direct numerical simulation ("Simulated difference") or lTRC calculation ("LTRC difference"). The durations in the default case (d 1 = d 2 = d 3 = 1) are t 1 = t 2 = t 3 = 29.3227. Phase 1 is prolonged in response to a positive d 1 -perturbation, while phases 2 and 3 are insensitive. Simulated difference ∆t 1 ∆t 2 ∆t 3 LTRC difference ∆t 1 ∆t 2 ∆t 3 µ = 0.05 0.1118 0.0008 0.0009 µ = 0.05 0.1033 0.0008 0.0008 µ = −0.05 -0.1107 -0.0009 -0.0008 µ = −0.05 -0.1033 -0.0008 -0.0008 The analytical results show good agreement with the numerical simulation, with both revealing that an increase in d 1 prolongs the active phase of cell 1 and leaves those of cells 2 and 3 almost unchanged.
The lTRC approach affords insight into the origin of the phase duration effects of perturbation. To capture these effects, we plot the integrand in (5) as well as the accumulated integral up to time t as functions of t. Specifically, for t ∈ [t in , t out ], the integrand depends on the three terms Since the exit boundaries are fixed, −η(x out ) · ∂x out µ ∂µ µ=0 is zero in this scenario. Figure 9(a) shows the timing responses to d 1 perturbation that occur in the three active phases over one period. We observe for phase 1 Figure 9: Timing responses for (a) intrinsic-release mechanism and (b) synaptic-release mechanism, measured by equation (5). Transitions into and out of each active phase are defined by (6), and the respective timing characteristics in each phase are calculated. Top panels: time course of membrane potentials over one period. Middle panels: integrand in T 1 at each time, given by (7). Bottom panels: accumulated integral of T 1 (the area between the integrand curve and x-axis). Black, blue and red traces correspond to when cell 1, cell 2 and cell 3, respectively, is active. Intrinsic release: When a positive perturbation is applied to d 1 , cell 1 (black) starts with a negative response upon entering the active phase at t = 83.49; the negative response continues accumulating, until t = 110.88 (where the integrand is zero); near the end of the active phase, the integrand curve rapidly increases to large positive values, which cancels out the negative response accumulated before and finally ends up with a positive response. Cells 2 and 3 react slightly to the perturbation only at their transition to the active state (at t = 53.71 and t = 23.93, respectively). Synaptic release: Although the integrand curve at the end of the phase 1 transiently increases to a value larger than for intrinsic release, its accumulated effect (bottom panels) is much smaller than for intrinsic release, leading to less variation in the active duration of cell 1.
(black trace) that the timing response just before the end of cell 1 activation is dominant while the effect of the perturbation on the entry to phase 1 is insignificant. The only effect of the d 1 change on phase 2 (blue trace) and phase 3 (red trace) is the slight modulation of their jump-up points. (We note that the local timing response curve in this setting has a qualitatively similar shape to the infinitesimal phase response curve (iPRC) (cf. Brown et al. (2004); Ermentrout and Terman (2010); Schwemmer and Lewis (2012)), and we discuss their relationship in Appendix B.2.) To qualitatively explain what the lTRC analysis reveals about the dominant effect of the jump-down of cell 1 on the response to perturbation, we consider the fast-slow decomposition of the relaxation oscillator system, which makes use of the interplay of spatial and temporal measures of distance between trajectories (Terman et al., 1998;Rubin and Terman, 2002). Recall that for the intrinsic release mechanism, the phase transition is triggered when the active cell reaches the appropriate right knee of the active v-nullcline. As d 1 increases and the right knee, or fold, of the v 1 -nullcline shifts to a lower position, cell 1 spends more time in the active phase. Specifically, as cell 1 progresses through the active state, it evolves toward the fixed point (the intersection of the free v 1 -nullcline and the h 1 -nullcline in Figure 7), which lies just beyond the knee, until the knee is reached. Thus, the dynamics of cell 1 exponentially decelerates as the projected trajectory approaches the jump-down point, as shown by the contraction of points around the right knee in Figure 7.
Hence, despite a small spatial difference for the v 1 -nullcline caused by a slight change in d 1 , the difference in time needed to cover that spatial spread, and hence the difference in the duration of active phase 1, is significant.
Note that as d 1 varies, the entry point of cell 1 into the active phase also changes (upper right corner of Figure 7). A larger d 1 results in a smaller h 1 upon entry of cell 1 to the active phase, and thus a shorter time spent in its active region. This reduction partially counters the delay around the right knee. However, since the voltage changes much more rapidly at the start of the active state than at the end, the spatial perturbation upon entry has negligible temporal effects.
Next, consider the effect of the d 1 -perturbation on the active phase of cell 2, the next cell to activate after cell 1. The change of d 1 does not affect the nullclines of v 2 or the jump-down point of cell 2 to the silent phase. The timing of the entry of cell 2 into the active phase is controlled by the jump-down of cell 1 from the right knee of the v 1 -nullcline, which releases cell 2 from inhibition. When cell 1 spends more time before the end of its active state, cell 2 correspondingly spends more time in the silent state. However, since the position of cell 2 is close to the silent phase fixed point (with extremely low rate of change), the added time in the silent phase does not make much of a difference in the spatial position at which cell 2 enters the active phase and hence its active phase duration.
Specifically, Figure 10 shows a schematic illustration of the jump-up of cell 2 when d 1 is positively perturbed. The longer time of cell 2 in the silent phase is equal to the extra time of cell 1 in the active phase, denoted by ∆t in the figure. However, when cell 2 enters the active phase, the time difference from its perturbed to its unperturbed position is highly compressed. We can in fact show that where t u and t p represent the unperturbed time and perturbed time that cell 2 spends in the active phase, respectively.
Consider the extra journey of cell 2 during the silent phase in the perturbed case. The evolution of the gating variable h 2 is approximated by for some values of h L and τ L . Hence, the extra time in the silent phase (from point h u to h p in Figure 10) is given by On the right branch, h 2 approximately follows for some value of τ R . Since the jump-down position of cell 2 is invariant with respect to the d 1 perturbation, the difference between the perturbed and unperturbed active durations of cell 2, which is also equal to the time of passage from the perturbed landing point of cell 2 in the active phase to its unperturbed landing The extra time spent by cell 1 to reach the right knee h RK after perturbation is equivalent to the extra time that cell 2 spends in the silent phase, traveling from h u to h p , denoted by ∆t. When cell 2 enters the active phase, the time difference between the red and blue landing points on the uninhibited v-nullcline is much smaller than ∆t; that is, t p − t u ∆t, where t u (resp., t p ) represents the passage time for cell 2 from its arrival in the active phase until it reaches the right knee of the v-nullcline without (resp., with) the perturbation. point, can be computed from To compare these quantities numerically, we choose h L = h ∞ (−62) and τ L = τ h (−62), based on a typical value of v 2 in the silent phase, and we assume τ R = τ h (−30) to give t p − t u its largest possible value. When (8) and (9) we obtain t p − t u = 0.0248, ∆t = 0.0999.
Note that because we chose the value of τ R that occurs when the trajectory is near the right knee, which represents the largest value that τ h (v 2 ) takes during the active phase, the difference t p − t u is in fact much smaller than 0.0248. This analysis confirms our temporal-spatial observation that the transition to the active state greatly compresses the temporal difference of cell 2 active phase durations (perturbed versus unperturbed), relative to the temporal difference of active phase durations for cell 1. The jump down of cell 2 releases cell 3 to the active phase, so cell 3 inherits the insignificant active phase duration change from cell 2. The fast-slow decomposition analysis above is in accordance with the lTRC characteristics (Figure 9(a)).
The duration of phase 1 changes monotonically with d 1 , and the durations of the other two phases are almost invariant.
Synaptic release Following a similar argument, we can uncover the sensitivity of the individual activation phases for the synaptic release mechanism. Note that the two release mechanisms have the same nullcline structure (see Figures 1, and 3) and we expect that very similar underlying mechanisms are at work in the two cases. Table 2 lists the duration changes with reference to the symmetric case. The duration of active phase 1 has a positive response to the increase of d 1 while the durations of phase 2 and phase 3 are insensitive. Compare the sensitivity of cell 1 activation between the two release mechanisms (Tables 1, 2). While both mechanisms yield increases in the phase 1 duration, ∆t 1 , the duration change for the synaptic release has a considerably smaller magnitude than ∆t 1 for the intrinsic release. This result is evident from the local timing response analysis as shown in Figure 9. The plot of the integrand of T 1 in equation (5) indicates that although the integrand at the end of phase 1 grows to a large positive value under the synaptic release mechanism, the elevation in the synaptic release case occurs much more briefly, resulting in a weaker temporal effect than results due to perturbation in the intrinsic mechanism.
In the synaptic release mechanism, the synaptic threshold θ I plays the key role in determining the transition of each cell from the active state to the silent state. As shown in Figure 8, as d 1 is adjusted, the voltage coordinate of the jump-down point of cell 1 remains almost equal to θ I while the coordinate of the gating variable h 1 changes. Because cell 1 only has to reach v 1 = θ I > v RK for the transition to begin, the smallest values of τ h from the intrinsic-release case do not occur, and hence the difference in active phase duration for cell 1 is smaller in the synaptic release case. This again translates into small differences in the active phase durations for cells 2 and 3.
Synaptic escape Unlike the release mechanisms, where the active cell itself has primary control of the transition, here the duration of the active phase is predominantly controlled by the silent cell. Specifically, each transition is initiated when one of the silent cells hits the synaptic threshold, thus inhibiting the active cell and inducing a transition. As d 1 increases, the inhibited v 1 -nullcine shifts downwards, and so does the escape position of the silent cell 1; the inhibited voltage nullclines, as well as the escape positions, of the other two cells are invariant. In the following, we first show the lTRC characteristics associated with the synaptic escape mechanism. Next, we provide a geometric argument to explain why the active phase duration of cell 1 increases while the active phase durations of cell 2 and cell 3 decrease when d 1 is increased. Finally, we strengthen our argument by giving an analytical calculation to show that under the assumption that t 1 increases (which can be justified by the geometric argument), both of t 2 and t 3 must decrease. Table 3 shows the changes in the three active phase durations when a static perturbation µ = ±0.01 is applied to d 1 . The first-order changes in the phase durations given by the lTRC formula (5) match well with direct numerical simulation. The timing response characteristics of the synaptic escape mechanism are  Figure 11: Timing responses for the synaptic escape mechanism, measured by equation (5). Colors are as in Figure 9. Here we adopt an artificial synaptic threshold θ imag to demarcate the start and end of the active phases. When d 1 is increased (by perturbation µ = 0.01), all three of the active phases have significant responses: cell 2 (blue) exhibits an initial negative deflection; cell 3 (red) starts with a positive response that is subsequently overcome by a more significant negative peak; cell 1 (black) shows a positive deflection when it enters the active phase.
shown in Figure 11. As d 1 increases, the silent cell 1 jumps up earlier, as indicated by the initial positive deflection of the black curve in the integrand plot. The earlier jump-up of cell 1 results in a lower position of the double-inhibition notch of cell 2 in its silent phase and a higher jump-down position of active cell 3, as shown in Figure 13 (see below for a more detailed discussion of Figure 13). These further lead to a later jump-up of cell 2 and an earlier jump-up of cell 3, which are in accordance with the initial negative deflection of the blue curve and positive deflection of the red curve in the integrand plot of Figure 11. Right before the end of the active phase of cell 3, a significant negative response occurs and counters the initial positive response. This negative response is due to the early escape of cell 1, which advances the termination of the active phase of cell 3. Figure 12 shows the phase portraits and trajectories obtained by numerical simulation for the synaptic escape mechanism. Note that there is a notch circled in each cell panel (see also Figure 2), arising from double inhibition from the other two cells (which plays an important role in changing the active phase durations). Figure 13 provides a schematic illustration of the phase portraits for the three cells. For visual convenience, we just show one v-nullcline per input level for each cell (i.e., we ignore the shift in v-nullclines due to inhibition). One can also refer to Figure 12 for the phase portraits and trajectories obtained by numerical simulation in this case. Red marks represent the locations of the cells resulting from an increase in d 1 and blue marks represent the corresponding locations with the default d 1 . Suppose d 1 is perturbed at the moment when cell 3 is at the jump-up position h U 3 (see panels (a)). Note that the jump-up of cell 3 causes cell 2 to jump down (from the right to the left branch of its v-nullcline). Since the v 1 -nullcline shifts downwards and the current position of cell 1, at the double-inhibition notch, denoted by h I 1 , moves to the right, the passage time of cell 1 to its escape position h U 1 (panels (b)) is shortened. Correspondingly, the active duration of cell 3, t A 3 , is shortened (or equivalently, the jump-down position of cell 3, h D 3 , is higher), and the doubleinhibition notch location of cell 2 on the inhibited v 2 -nullcline, denoted by h I 2 , is lower (panels (b)). In panels (c), since the escape position of cell 2, h U 2 , is invariant with respect to d 1 , the passage time of cell 2 to from h I 2 to h U 2 becomes longer, which results in a longer t A 1 , a lower h D 1 , and a higher h I 3 . Such higher h I 3 leads to a shorter passage time of cell 3 to its escape position, and thus, t A 2 is shortened, h D 2 is higher, and h I 1 is lower (see panels (d)). Note that in the cell 1 panel of (d), we use the blue dot to mark the location of the red dot from panel (a) that resulted from the initial perturbation of d 1 . Since h I 1 is lower at the end of the cycle, in panel (d), than it was initially, if we continue to track the cells following a similar argument, then we see that the active t A 3 in the second iterate is no longer quite as short as it was in the first iterate.
To represent this situation analytically, let t A i (k) denote the active phase duration of cell i on iterate k, where k = 0 corresponds to the unperturbed solution. From above, we have (2). Since the active duration of cell 3 is equal to the passage time of cell 2 from point h D 2 to h I 2 in the silent phase, it is sufficient to consider the h 2 -dynamics of cell 2. We approximate this dynamics with the ODE h 2 = (h L − h 2 )/τ L , from which we obtain Taking the derivative of t A 3 with respect to the perturbation d 1 yields .
As discussed above, Repeating similar arguments, we obtain Continuing in this way, we get a sequence of active phase times, bounded above or below as needed, and we can conclude that when subject to a static positive d 1 perturbation, t A 1 increases while t A 2 and t A 3 decrease.
The decrease of t A 2 and t A 3 can also be obtained by a direct calculation. Consider the dynamics of cell 1, with active (silent) duration t A 1 (t S 1 ) and with h 1 -coordinate ranging from h U 1 to h D 1 . Assume that t A 1 increases when d 1 is positively perturbed. As shown in Figure 13, the perturbation lowers h U 1 and thus lowers h D 1 as well. In the silent phase, the time of cell 1 from h D 1 to h U 1 is t S 1 = t A 2 + t A 3 , and thus, by approximating the ODE for h 1 , we have .
Due to the compression around the right knee of the active v 1 -nullcine, the magnitude of the decrease in h U 1 (d 1 ) is larger than the magnitude of the decrease in h D 1 (d 1 ). Hence, x(d 1 ) increases, and (t A 2 + t A 3 ) decreases. Now, suppose that t A 3 increases; correspondingly, t A 2 must decrease. On the one hand, this indicates that t S 2 = t A 1 + t A 3 increases. On the other hand, since t A 2 is shortened and h U 2 is independent of d 1 , then h D 2 must be higher. Therefore, t S 2 must decrease, which is a contradiction. Hence, t A 3 decreases. A similar argument shows that t A 2 also decreases. In this way, we conclude that when the active phaseduration for cell 1 increases, the active phase durations for both cell 2 and cell 3 must decrease. This geometric analysis and direct mathematical calculation show agreement with the timing response analysis shown in Figure 11.   Figure 2. The system exhibits a stable short-long alternation and two double-inhibition notches, with the nontrivial Floquet multipliers less than one in magnitude.

Intrinsic escape
The system with phase transitions via the intrinsic escape mechanism has an unstable equal-duration firing pattern ( Figure 2) and multiple stable short-long patterns (e.g., Figure 14, generated with the same parameter set as Figure 2). Numerically we find that the equal-duration solution in Figure 2 has two unit Floquet multipliers, while for the short-long solution in Figure 14, all Floquet multipliers except the trivial Floquet multiplier are strictly smaller than one in modulus. Due to the particularity of this mechanism, we do not give its sensitivity analysis. Appendix C presents analysis of the instability of the equal-duration periodic solution for the intrinsic escape mechanism in the singular limit → 0.

Heteroclinic cycling model
In this section, to study the controllability of the heteroclinic cycling model (2), we apply a small sustained perturbation µ to a 1 = 0.01 (i.e., a 1 → a 1 + µ), with the other two control parameters a 2 and a 3 fixed at the default value of 0.01. As a reference, Figure 15 shows how the dynamics of the three pools evolves when a 1 = 0.0005, 0.01, or 0.02. In what follows, we explain the sensitivity of each phase duration using the lTRC as well as direct analysis based on the vector field. Figure 16 shows the timing responses of the heteroclinic cycling system computed from the lTRC ias given in equation (5). Phase 3 (red) has a significant negative response to the increase of a 1 ; phase 2 (blue) is almost invariant; phase 1 (black) has a relatively small negative response. To understand the effect upon trajectory entry to each phase, consider the terms comprising "term A" in equation (5), which captures the effects of shifting the entry point to a given region, namely η(x in ) · ∂x in µ ∂µ µ=0 . The timing response curve vector at the entry point to phase 3 is η III (x in ) = (−48.89, 0, 0.06) , which shows that the first component should dominate the response. The lTRC at the entry points to the other phases is given by cyclic permutation, so the contribution to the response at the entry to phase 2 is dominated by the third component, and the response upon entry to phase 1 is driven by the second component. The timing sensitivity depends on the inner product between the lTRC and the derivative of the entry location with respect to µ. The lTRCs are symmetric (under threefold rotation); the derivatives of the entry location vary from phase to phase.  Table 4: Heteroclinic cycling model: duration changes in three individual phases, from direct numerical simulation or lTRC calculation. The durations in the default case (a 1 = a 2 = a 3 = 0.01) are t 1 = t 2 = t 3 = 2.9080. With a small positive perturbation to a 1 , pool z shows a significant decrease in the active duration (∆t 3 ); pool y is insensitive (∆t 2 ); pool x has a relatively small negative response (∆t 1 ).
Simulated difference ∆t 1 ∆t 2 ∆t 3 LTRC difference Specifically, the entry location on phase 3 varies most, while the entry location on phase 2 varies least (as a reference, see Figure 15(b)). This accounts for the smaller or larger initial displacement effects seen in Figure 16 (grey arrows). Table 4 lists the duration changes in the three individual phases, obtained either from direct numerical simulation or from the lTRC formula (5), when a perturbation of size |µ| = 0.0005 is applied to a 1 = 0.01.
The vector field of the model is given in (2). During phase 1 (in which pool x is active), a positive perturbation to a 1 only influences the evolution of x and the boundary surfaces of phase 1. When dx/dt > 0, the increase of a 1 decreases the speed of x, making x grow slower; when dx/dt < 0, the increase of a 1 accelerates the decline of x. Note that the exit surface of x from phase 1 is x = y + a 1 +a 2 2 (and the trajectory stays in phase 1 when x ≥ y + a 1 +a 2 2 ). Compared with the default case, the effects of increasing a 1 on dx/dt cause the x on the left hand side of this condition to shrink, while the right hand side increases with a 1 .
With perturbation a 1 → a 1 +µ, denote the perturbed (x, y, z)-solution and the unperturbed (x, y, z)-solution to be (x p (t), y p (t), z p (t)) and (x u (t), y u (t), z u (t)), respectively. Let (x 0 p , y 0 p , z 0 p ) and (x 0 u , y 0 u , z 0 u ) represent the entry points for the perturbed case and unperturbed case, respectively. Let t 1 be the unperturbed exit time from phase 1 and t * 1 = t 1 + ∆t 1 be the perturbed exit time; from the above observations, ∆t 1 < 0. We Figure 16: Timing responses in the heteroclinic cycling model (2). Black, blue and red traces represent pools x, y and z, respectively. When a positive perturbation is applied to a 1 , pool z has a significant negative timing response throughout its active phase. The initial displacement at the beginning of each phase (grey arrows) indicates the effect of the perturbation on the entry point to each phase (given by term A in equation (5)).
analytically derive that the first-order approximation for the duration change in phase 1 can be estimated by (see Appendix D for the derivation) When µ = 0.0005, from a direct numerical simulation, we have t 1 = 2.9080, x u (t 1 ) = 0.3874, x 0 u = 0.3773, x 0 p = 0.3769, y 0 u = 0.0111, y 0 p = 0.0112, and (10) gives ∆t 1 ≈ −0.0071. Compared with Table 4, this result confirms that equation (10) provides an accurate estimate to the change of phase 1 duration. We note that based on numerical calculations, the fidelity of this estimate depends on the inclusion of the impact of the perturbation on the positions of the entry points into the phase (i.e., x 0 p , x 0 u , y 0 p , y 0 u ), which is consistent with the initial negative deflection of the black curve in Figure 16 (grey arrow).
For phase 2, in which pool y is active, the entry points for different values of a 1 are extremely close (see Figure 15 when a very large perturbation is applied). The change in a 1 influences neither dy/dt nor dz/dt, which also do not depend on x. Further, the exit surface, given by y = z + a 2 +a 3 2 , is also independent of a 1 . Therefore, the duration of phase 2 is insensitive to the a 1 -perturbation.
When z is active in phase 3, the increase of a 1 leads to a faster increase of x, which then results in a faster decay of z. The exit surface is z = x + a 1 +a 3 2 and the trajectory stays in phase 3 when z > x + a 1 +a 3 2 . Hence, a decrease in z on the left hand side of this condition and an increase in the right hand side due to both x and a 1 indicate that the trajectory leaves phase 3 earlier than in the default case. Similar to phase 1, we obtain From direct numerical simulation, and by (11), ∆t 3 ≈ −0.0448. Consistent with the direct numerical simulation and lTRC results in Table 4, equation (11) gives a good estimate of the phase 3 duration change. Numerical calculations suggest that similarly to phase 1, the perturbation of the entry position to phase 3 has a significant influence on the phase 3 duration (as suggested by the large initial displacement of the red curve in Figure 16 (grey arrow)).
Note that such vector field analysis is of limited general utility because it requires the solutions of the model ODEs. Since the results of this analysis give a good agreement with the timing response calculation, however, this example provides a reassuring demonstration of the relevance of the lTRC analytical formula (5), which can be used to compute responses to perturbations in general cases.

Competitive threshold-linear network
In the competitive threshold-linear system (3), we impose a small static perturbation µ on the drive term θ 1 to the first node x 1 . As an illustration, Figure 17 shows the trajectories in the default case where θ 1 = θ 2 = θ 3 = 1 and in the perturbed case where θ 1 is positively perturbed at the beginning of the first active phase. (Note that for convenience of intuitive understanding, the perturbation here is applied since the first active phase, but in the following calculation we still impose a static perturbation along the limit cycle.) We can understand the effects of perturbation as shown in Figure 17 by referring to equation (3). The increase in θ 1 implies that W 12 x 2 + W 13 x 3 must become more negative (recall that W 12 , W 13 < 0) in order to causė x 1 = 0 than was needed in the unperturbed case. Although x 2 begins to increase while x 1 is still active, even if x 2 increased at its unperturbed rate, it would take longer for x 2 to grow enough for W 12 x 2 + W 13 x 3 to reach this more negative value; moreover, x 2 actually becomes slowed once x 1 exceeds its unperturbed maximum. The longer time until x 1 starts decreasing and the larger maximal amplitude of x 1 translate into a longer active phase for cell 1. Meanwhile, again with reference to equation (3) with W 31 = −1 − δ, we see that the larger amplitude of x 1 makesẋ 3 stay negative for a longer time than before and causes x 3 to reach a lower minimum value. From this lower minimum, x 3 inhibits x 2 only weakly, allowing x 2 to stay larger than before, and it takes longer for x 3 to overtake x 2 than in the unperturbed case, resulting in a longer active phase for cell 2. Meanwhile,ẋ 1 becomes positive earlier relative to the active phase of x 3 than before, due to the increase in θ 1 . The resulting earlier increase in x 1 causes an earlier sign switch ofẋ 3 at a smaller x 3 value, resulting in a shorter active phase for cell 3.
Using the lTRC calculation, Table 5 lists the duration change for each individual active phase in response to the perturbation of magnitude |µ| = 0.01 on θ 1 , and Figure 18 shows the local timing responses. Consistent with direct analysis above, all three phases are sensitive to the perturbation: the first two phases are prolonged under a positive perturbation, while the third phase is shortened.
. With a small positive perturbation to θ 1 , nodes x 1 and x 2 display a positive response in their active duration; node x 3 displays a large negative response. With a small negative perturbation to θ 1 , the opposite responses occur. Simulated difference ∆t 1 ∆t 2 ∆t 3 LTRC difference We note that the phase duration responses are quite different from those in the heteroclinic cycling case (compare Table 4). The heteroclinic cycling case allows independent phase modulation -selectively increasing the active duration of i-th unit can be attained by decreasing the drive input to the (i+1)-st unit. In contrast, in the competitive threshold-linear case, implementing perturbation on one phase gives rise to simultaneous changes in the other phases. A key difference between these two cases is that the heteroclinic cycling model features the passage near a saddle point in each phase, which the threshold-linear model lacks. Thus, the difference in responses to perturbations between these cases suggests that the slowing of trajectories near saddle points may make an important contribution to the timing responses to perturbation in the heteroclinic cycling case.

Discussion
We have considered three different dynamical systems frameworks for the production of triphasic oscillations in which three components take turns becoming active in a fixed order. In each framework, we have presented a computational and mathematical analysis of the effects induced by changes in a parameter, which represents the strength of tonic input to one component, on the active phase durations of all of the components in the network. We find that distinct phase duration responses arise when phase transitions occur through different mechanisms. Specifically, we arrive at the following observations.
Relaxation oscillator model Each of the four mechanisms underlying the relaxation oscillator model features its own phase modulation characteristics.
• Intrinsic release: an increase in the excitatory gain factor d 1 prolongs the active phase of cell 1 and leaves those of cells 2 and 3 almost unchanged. The mechanism involves the fast-slow dynamics of the model. Although the increased d 1 advances the entry of cell 1 to the active phase, the small spatial change around the jump-down (exit) position of cell 1 from the active phase translates into a significant temporal extension. As cell 1 progresses through the active state, the dynamics exponentially decelerates towards the exit position.
• Synaptic release: similarly to intrinsic release, an increase in d 1 increases the duration of cell 1's active phase, but it is a quantitatively smaller effect in the synaptic release case than in the intrinsic release case. The mechanism by which increasing d 1 affects the synaptic release case is qualitatively similar to that of the intrinsic release case. The smaller effect size is due to the different jump-down positions in the two release mechanisms; in the synaptic case, the jump-down of cell 1 is farther from the fixed point, with less slowing down, resulting in the less significant temporal increase.
• Synaptic escape: an increase in d 1 prolongs the active phase of cell 1 and shortens the active phase durations of cells 2 and 3. In contrast to the release mechanisms, in this case an increase in d 1 allows cell 1 to enter its active phase earlier in the cycle than it otherwise would, leading to several changes in the points at which the cells in the circuit initiate jumps between phases and to the compound effects on active phase durations of all three cells.
• Intrinsic escape: although the two-unit half-center oscillator (HCO) with intrinsic escape produces stable oscillations, the 3-cell circuit produces an unstable equal-duration pattern and multiple stable short-long alternating patterns. In this study, we focus on patterns with equal phase durations, and thus the response sensitivity of the stable patterns in this case are outside of the scope of our analysis.
Overall, we conclude that in the release cases, the system can selectively increase the active duration of a single unit with minimal change in the active phases of the other units, just by increasing the drive to that specific unit. However, in the synaptic escape case, an increase in the active duration of one unit is accompanied by non-negligible decreases in the active durations of the other units. In other words, independent active phase modulation is possible with release transitions but not escape transitions.
Heteroclinic cycle model An increase in the excitatory neural drive a 1 leads to a slight decrease in the active phase duration of the first cell, a negligible change in the active phase duration of the second cell, and a large decrease in the active phase duration of the third cell. The mechanism appears mainly to involve shifting the entry point of the trajectory to the domain governed by the vector field around the saddle point corresponding to elevated activity in the cell receiving the increased drive. Nearly independent active phase modulation can be achieved in this model; in contrast to the relaxation oscillator model, selectively increasing the active duration of the i-th unit requires decreasing the input to the succeeding, (i+1)-st unit.
Competitive threshold-linear network model Increasing the external drive θ 1 to unit 1 prolongs the duration of the active phases of both cell 1 and cell 2, while decreasing the active phase duration of cell 3. This mechanism increases the active duration of the i-th unit by decreasing the drive to the (i+1)-st unit, with the side effect of increases in the active durations of the other uints. Hence, a key difference of this model from the heteroclinic cycling case is the lack of a simple mechanism for independent phase modulation. A perturbation can be implemented here such that one cell's active duration exhibits the predominant effect, but the timing of the other cells is still impacted. This contrast in phase modulation properties is likely due to the lack of saddle points, which make an important contribution to the timing of trajectories in the heteroclinic cycle model.
In the relaxation oscillator cases, phase duration changes manifest through fast threshold modulation and associated timing effects. Here, we can think of the units that are not directly perturbed as being unaffected (i.e., no change in flow or transition surfaces) except in the timing at which their transitions occur. Consistent with past work (Daun et al., 2009;Rubin et al., 2009), however, the specific transition mechanisms do matter for what outcomes result from perturbations. In the piecewise continuous heteroclinic case with state-dependent switching boundary surfaces (i.e., surfaces defined in terms of the dependent variables), the effects could arise both through changes in switching manifolds and through direct effects on the vector field. To distinguish these factors, we also consider: • The piecewise-continuous heteroclinic case with fixed boundaries. Here, we obtain almost the same results (not shown) as in the switching boundary case. This indicates that the dominant contribution to the duration changes comes from the vector field, especially the effect of changes in the vector field on the positions where the trajectory reaches the surfaces of vector field discontinuity, not from changes in these surfaces themselves.
• The continuous heteroclinic case with fixed boundaries. This model has the fewest possible factors affecting the timing. We find that it has qualitatively the same responses (not shown) to perturbation as the piecewise-continuous cases. Thus, the impact of a perturbation does not arise predominantly from the jumps in the vector field, but rather from the impact of the perturbed vector field on the trajectory position, irrespective of the vector field continuity.
Unlike the relaxation oscillator cases, in the heteroclinic cycle and threshold-linear models, we cannot decompose the flow into lower-dimensional systems for separate units, and impacts on the flow occur throughout the path of the trajectory. Interestingly, it is not obvious (at least to us) why these combinations of factors result in the diverse responses to perturbations that arise across the cases that we have considered.
A previous study performed a thorough analysis comparing effects of perturbations to single units within a fast-slow half-center oscillator model with either escape transitions, release transitions, or transitions via an intermediate "adaptation" mechanism combining elements of both escape and release (Daun et al., 2009). Interestingly, that study found that escape transitions restricted the impact of such perturbations to the silent duration of the perturbed unit or equivalently the active duration of the unperturbed unit, while release transitions led to compensatory effects and produced only weak changes in durations, and in the adaptation regime, one-sided perturbations affected durations of both units. In other words, independent silent phase modulation was possible with escape transitions but not release transitions. In fact, our results are consistent with these earlier findings. In this study, in the escape regime, we do have independent modulation of silent phase durations. Specifically, as shown in Table 3, the duration change in the silent phase 1 is given by ∆t 2 + ∆t 3 , which is significantly negative when a positive perturbation on d 1 is imposed, while the duration change in the silent phase 2, ∆t 1 + ∆t 3 , and the change in the silent phase 3, ∆t 1 + ∆t 2 , are negligible. Moreover, in the intrinsic release regime, the larger duration change that we found is linked to the slowing down near the right knee of the voltage nullcline, such that the difference between the two papers in the release case comes not from the number of cells but rather from the position of the slow nullcline relative to the right knee.
Applying similar arguments, we can consider how well each model analyzed in this study achieves inde- pendent active phase modulation as could be important for downstream effectors driven by excitation from CPG units, independent silent phase modulation as could be important for downstream effects dependent on pauses in CPG unit activity, and either invariance or sensitivity of period to perturbations. The information in Table 6 allows us to discuss these properties heuristically and empirically, where for the latter, we define a circulant matrix C for each case with first row given by the elements of columns 1, 2, 3 of the table. We see that the most pronounced independent active phase modulation comes from the intrinsic release model, in which an input perturbation to cell 1 induces an increase in the cell 1 active phase duration that is strong relative to the negligible impacts on the other two cells' active phases. As discussed above, the synaptic escape model achieves independent silent phase modulation. The threshold-linear network also exhibits this property to some extent, because changes in input to unit 1 will produce much stronger impacts on the silent phase duration of unit 3 (comprising the active phase durations of units 1 and 2) than on those of the other units. Finally, the ranking of the models in order of sensitivity of period to input modulation goes as follows: synaptic escape, heteroclinic cycle, threshold-linear, intrinsic release, synaptic release. This ranking aligns with C −1 in Table 6. Indeed, if a system wants to induce a change ρ in its three active phase durations, then that corresponds to finding a vector µ of input strengths such that ρ = Cµ. The size of the needed input is therefore bounded as µ ≤ C −1 ρ . The threshold-linear model is interesting in that it allows for significant phase perturbations with relatively small changes in overall period, as evidenced by its large condition number (see cond(C −1 ) in Table 6). On the other hand, since all of the rhythm features in the synaptic release model show little response to input changes, that framework seems to be effective for robustness but not amenable to simple control; more generally, the relatively weak responses of the release models to input perturbations is consistent with past observations in HCOs (Daun et al., 2009).
On the other hand, a new result from our work that distinguishes the triphasic case from the HCO case is the instability to perturbations for the triphasic rhythm with intrinsic escape transitions. This finding does not rule out intrinsic escape as a transition mechanism for biological CPGs, which may feature asymmetries or differences in transition mechanisms for different phases as discussed below, but it does highlight the importance of checking robustness of rhythmic behavior in models, of paying attention to transition mechanisms in model design, and of taking care in extending small circuit models to larger ones.
A paradigmatic model for the brain stem neuronal network underlying respiratory rhythm generation includes three populations coupled with mutual inhibition, along with a fourth excitatory population. An investigation of that model based on fast-slow decomposition and bifurcation analysis revealed that the expiration-to-inspiration transition occurs via escape while the inspiration-to-expiration transition results from release. Simulations yielded the predictions that varying a tonic control drive to the escaping inspira-tory neurons would shorten the expiratory phase, while varying the drive to the released expiratory neurons would have relatively little effect on phase durations . These results are consistent with the two-cell, half-center oscillator analysis (Daun et al., 2009), but it is important to note that the active phases of two of the three inhibitory respiratory populations overlap, so it is not surprising that this system acts more like a two-cell network than the type of three-cell network that we studied. Although the results from the respiratory model, which combines escape and release, agreed qualitatively with those obtained from analysis of pure escape and pure release scenarios, it is not obvious that a similarly straightforward generalization of our analysis would yield correct predictions for three-cell networks that combine escape and release.
In contrast to these studies of effects of constant input levels on phase durations, Zhang and Lewis investigated differences in the infinitesimal phase response curves (iPRC) for HCO systems (Zhang and Lewis, 2013). They found that the timing sensitivity of release systems to transient pertubations was greatest when applied to the active cell, near the end of the active phase. Similarly, they found that escape systems showed greatest sensitivity to perturbations applied towards the end of an escaping cell's silent phase. Here, we find behavior that is consistent with these observations. In Figure 9, for example, under the release mechanisms, the integrand of the local timing response curve integral was largest (greatest sensitivity) towards the end of the active phase of the cell (cell 1) subject to the perturbation (changing d 1 ). In contrast, in Figure 11, under the escape mechanism, the largest values of the lTRC integrand occurred towards the end of the inactive phase of the cell coming one step earlier in the activation sequence than the cell subject to the parametric perturbation.
Extreme sensitivity of timing for limit cycle trajectories that pass close to a saddle or are near a heteroclinic or homoclinic bifurcation point has been noted previously. For example, Shaw et al. (2012) pointed out phase-dependent divergence of infinitesimal phase response curves for a family of limit cycles that approach a heteroclinic cycle as a parameter is varied. In another timing-related analysis, Izhikevich derived the phase equations for a coupled pair of relaxation oscillators and showed that in the singular limit, the effective coupling function exhibits discontinuous jumps at the points where fast phase transitions start and stop; these jumps perturb into large slopes when the singular perturbation parameter is small but nonzero (Izhikevich, 2000).
This work offers a stepping stone towards a broader theory of control of limit cycle systems. Control theory is most fully developed for linear systems (Brockett, 2015;Stanhope et al., 2014;Bechhoefer, 2021). For a linear system such asẋ = Ax + Bu, with state vector x ∈ R n , control input u ∈ R k , and observable y = Cx ∈ R m , the notions of controllability and observability are global: one constructs the "controllable" space C ⊂ R n of all points to which the system can be driven in finite time from initial x(0) = 0 via a suitable control, and the "observable" space O ⊂ R n such that trajectories from different initial conditions within O may be distinguished by observing the trajectory y(t). For nonlinear systems these notions are necessarily more limited, and are restricted to specific subclasses of systems (Isidori, 1985;Whalen et al., 2015). Here we consider nonlinear dynamical systems supporting stable limit cycle trajectories, as typically arise in central pattern generator models. Previous work on control of limit cycle systems has focused on controlling global properties of limit cycles, for instance the amplitude or period of an ongoing oscillation, or (de)synchronization of a population of coupled oscillators Wilson and Moehlis, 2022). Here we focus on the control of specific subcomponents of a limit cycle trajectory, in motor control systems comprising three distinct phases of movement. Classic analysis for understanding timing of control systems under a weak perturbation is well established by using the infinitesimal phase response curve (iPRC). The iPRC can only capture the global change in the timing induced by the perturbation, instead of the local timing change in each subcomponent of the limit cycle.
To circumvent this limitation, Wang et al. (2021) developed the local timing response curve (lTRC), which allows one to compute local timing changes due to nonuniform timing sensitivities of the control system along the limit cycle. In this work we have conducted analysis on the duration sensitivity in each phase of limit cycle systems by applying the lTRC, which yields excellent agreement with both direct numerial simulation and other mathematical analysis. Our work provides detailed examples of how to analyze the responses to perturbations arising in multi-phase oscillations of naturally occurring biological motor control systems, by applying tools that allow linearized/variational analysis of the parametric variations.
From a biological point of view, a general important point evident in this and past work is that we cannot assume that effects of perturbation of one unit in a coupled circuit will be localized to that unit or even dominant in that unit. Moreover, such effects may not materialize instantly; instead, we may need to observe one or more full cycles to see the persistent effects. In addition, the lessons we extract from this work can provide insights into the selection of a CPG model for robotics or for representing biological systems. Specifically, the divergent properties of the heteroclinic cycle CPG model and the relaxation oscillator CPG models can guide the selection of a suitable modeling framework to capture experimental observations or design aims for triphasic rhythmic systems, while similar analysis can be applied to other models producing other forms of multi-phase rhythms.
Note that circuits of interest for biological applications would typically not have the same degree of symmetry as those studied here. For example, in the respiratory network, in addition to the possibility that both escape and release transitions arise during each respiratory cycle, it is unlikely that the effective strengths of all synaptic connections within the network are equal; indeed, these differ in models that meet experimental benchmarks (Golowasch et al., 1999;Latorre et al., 2002;Rubin et al., 2009Rubin et al., , 2011Park and Rubin, 2022). It remains an interesting open question to apply the methods developed here to more realistic circuits. The work presented here serves as a template, although it remains to be seen to what extent the lessons learned here carry over to systems with heterogeneous switching mechanisms.
In this paper we have considered deterministic CPG models. In biological control circuits, stochastic fluctuations can also play an important role (Carroll and Ramirez, 2013;Yu and Thomas, 2021;Rubin et al., 2022). For example, stochastic heteroclinic channels have been studied intensively (Stone and Holmes, 1990;Armbruster et al., 2003;Bakhtin, 2011). Horchler and colleagues used the amplitude of an additive noise signal as a control parameter to dynamically regulate the mean period of a three-population May-Leonard model (Horchler et al., 2015). It would be interesting to investigate how changing the amplitude of noise impinging on individual units affects the (mean) activation durations of all three units, in analogy to the single-input perturbations studied here, but this question goes beyond the scope of the current paper. In his doctoral thesis, Shaw observed that one could manipulate the mean period of a stochastic 2D heteroclinic-channel based oscillator by changing the underlying vector field, by varying the additive noise level, or through a combination of both (Shaw, 2014). Results in Barendregt and Thomas (2022) put this analysis on a rigorous footing and extend it to several variants of the May-Leonard three-population model with discrete (non-Gaussian, non-additive) noise arising from demographic stochasticity.
Finally, we have presented three dynamical architectures that produce triphasic rhythms but our examples are certainly not exhaustive. For instance, the repressilator is a genetic regulatory network consisting of three transcriptional repressor genes (Elowitz and Leibler, 2000). The three units are connected in a negative feedback loop, allowing sustained limit-cycle oscillations. Another example is the stomatogastric ganglion of decapod crustaceans. Its constituent neurons encompass two central pattern generators, one for the gastric as in text None mill rhythm and one for the pyloric rhythm. For the pyloric CPG, three nerves, containing all of the pyloric units, generate a triphasic rhythm. This pattern controls striated muscles that dilate and constrict the pyloric region of the stomach, which moves the food caudally through filter combs inside the pylorus (Weimann et al., 1991;Selverston et al., 1998). The approaches we use here can be applied to these triphasic systems as well.

A Model details
In the relaxation oscillator system (1), the functions x ∞ (v) for x = h, m p are given by . Table 7 lists the parameter values for the intrinsic release transition mechanism, which produces Figure 1. Based on Table 7, change θ I to -25 mV to generate the synaptic release mechanism ( Figure 3); change θ I to -62 mV and σ h to 5 to generate the synaptic escape mechanism (Figure 4). For the intrinsic escape mechanism shown in Figure 2 and Figure 14, make the following changes: θ I = −36, θ h = −39, σ h = 9, g I = 0.24, g E = 0.14.
The simulation codes are available at https://github.com/zhuojunyu-appliedmath/Triphasic-control. Instructions for reproducing each figure and table in the paper are provided.
B Local timing response curve (lTRC)

B.1 Introduction to lTRCs
Consider a parametrized continuous-time dynamical system defined on a domain Ω ⊂ R n , where x ∈ Ω, µ ∈ R and F µ (·) is either smooth or piecewise smooth on Ω. We are interested in the case where F µ (·) admits a family of hyperbolically stable limit cycles for some range of µ including µ = 0, which we consider the "unperturbed" limit cycle γ(t). We assume the basins of attraction are contained in Ω and that the domain is partitioned into two or more subdomains, Ω = R I ∪ R II ∪ · · · with successive domain boundaries transverse to the flow of the unperturbed limit cycle. Moreover, we suppose that F µ (·) is smooth (in x) within each domain. (See Park et al. (2018) for a fully detailed construction of such a system of domains.) In the following we will write F (x) for F 0 (x). The local timing response curve (lTRC) introduced by Wang et al. (2021) is defined to measure the timing sensitivity within each region in response to a static perturbation µ = 0.
For x ∈ R i , i = I,II, · · · , let Γ i (x) be the time remaining until the trajectory beginning at x exits the region. At any point x ∈ R i , by construction, Γ i satisfies along γ(t). Thus, by the chain rule The associated lTRC for region R i is defined to be the gradient of Γ i evaluated along the limit cycle: For convenience, we drop the superscript i but note that every region is characterized by its own lTRC. Differentiating both sides of (14), within any particular subdomain, with respect to t gives where DF denotes the Jacobian of F . Let x out denote the exit point from region R i and n out denote a normal vector of the exit boundary of R i at x out . Then, following (14), the lTRC η satisfies the boundary (normalization) condition The adjoint equation (15) together with the boundary (normalization) condition (16) defines the local timing response curve in each subdomain.

B.2 Relationship between iPRC and lTRC
The formulation of the lTRC given by (15) and (16) is similar to the formulation of the infinitesimal phase response curve (iPRC). The iPRC, denoted by z(t), satisfies the same adjoint equation with a different normalization condition F (γ(t)) · z(t) = 1. Figure 19 shows the lTRC and iPRC during the first active phase for the relaxation-oscillator system with each of the two release mechanisms. (Due to symmetry, the curves for phase 2 (or phase 3) are the same if the same perturbation is applied to v 2 (or v 3 ).) As expected, the lTRC and iPRC have an identical geometric shape and are mutually symmetric about the time axis.
B.3 Derivation of generalized T 1 formula (5) This section presents the derivation of equation (5), which is a generalization of (4) from the fixed-boundary case to the shifting-boundary case. It follows the same generally exposition given in Wang et al. (2021) but with some modifications. Like the lTRC, the linear time shift T i 1 is considered separately in each subdomain R i , but the superscript i is dropped for the time being.
Suppose limit cycle γ µ is a solution to the perturbed system where τ = τ (µ, t) is the perturbed time coordinate satisfying τ (0, t) = t, τ (µ, t + T 0 ) = τ (µ, t) + T µ , where T 0 and T µ represent the durations spent by the unperturbed trajectory and the perturbed trajectory, respectively, in region R i . For simplicity, consider the linear scaling Suppose the unperturbed trajectory γ(t) enters R i through boundary Σ in at x in when t = t in , and it exits the region through Σ out at x out when t = t out ; we use similar notation for the perturbed trajectory γ µ (τ ). Assume the boundaries are not fixed, but perturb to Σ in µ and Σ out µ under the static perturbation µ, as illustrated in Figure 20. Without loss of generality, let t out = 0, and thus, the time duration the limit cycles spend in R i is T 0 = −t in , T µ = t out µ − t in µ . Therefore, at entry and exit points to region R i , the following expressions give the times that the unperturbed and perturbed limit cycles have left to spend in region R i : As in (13) and (14), for any µ, − 1 = dΓ µ (γ µ (τ )) dτ = F µ (γ µ (τ )) · η µ (γ µ (τ )).
Integrating (19) from τ = t out µ to t in µ gives On the other hand, we can expand T µ as where T 1 is the first-order change in the duration the trajectory spends in region R i resulting from perturbation µ. Our goal is to obtain the expression for T 1 using Taylor expansions for all of the terms in (20).
Suppose we can expand F µ , Γ µ , η µ and γ µ as follows: Figure 20: A schematic illustration of the shifting entry and exit surfaces for a system with three regions. The unperturbed trajectory γ(t) (black circle) enters region R I (grey shaded) through boundary Σ in and exits R I through boundary Σ out ; the perturbed trajectory γ µ (t) (red circle) enters region R I µ (magenta shaded) through boundary Σ in µ and exits R I µ through boundary Σ out µ . The grey dot on Σ in (resp., Σ out ) represents the entry position x in (resp., exit position x out ) of the unperturbed trajectory into (resp., from) region R I . The magenta dot on Σ in µ (resp., Σ out µ ) represents the entry position x in µ (resp., exit position x out µ ) of the perturbed trajectory into (resp., from) region R I µ .
Thus, we finally obtain the expression for the first-order change in the duration through region R i : which is exactly equation (5).

C Instability for intrinsic escape
This section gives an analysis of the instability of the equal-duration firing pattern for the intrinsic-escape transition mechanism. We consider the case in the singular limit → 0, as sketched in Figure 21. (Note that compared with Figure 2, there is no double-inhibition notch.) In this activity pattern, all cells jump up from the left knee of the inhibited nullcline. Suppose cell 1 (black) is ready to escape and jump up. Let x denote the passage time of cell 3 (red) from its jump-down position to its escape position (left knee of the inhibited v 3 -nullcline). Let y denote the passage time of cell 2 (blue) from its initial position to the escape position. Correspondingly, the time cell 1 spends in the active phase is equal to y. Since all cells jump up from a specific position, we can define the map f (y) that translates the time y that a cell spends in the active phase to the time that it spends in its next silent phase.
Each phase within a full evolution cycle starts with an analogous configuration, with a "lower silent" cell that just jumped down, an "upper silent" cell that will be next to escape, and an escaping cell. Table 8 gives the passage time for the two cells in the silent phase as time passes; the left column denotes the total time  x y cycle 2 f (y) . . . . . . . . .
Substituting the first two equations into the third equation and using Taylor expansion, with cancellation