A Scalar Poincaré Map for Anti-phase Bursting in Coupled Inhibitory Neurons With Synaptic Depression

Short-term synaptic plasticity is found in many areas of the central nervous system. In the inhibitory half-center central pattern generators involved in locomotion, synaptic depression is believed to act as a burst termination mechanism, allowing networks to generate anti-phase bursting patterns of varying periods. To better understand burst generation in these central pattern generators, we study a minimal network of two neurons coupled through depressing synapses. Depending on the strength of the synaptic conductance between the two neurons, this network can produce symmetric n : n anti-phase bursts, where neurons fire n spikes in alternation, with the period of such solutions increasing with the strength of the synaptic conductance. Relying on the timescale disparity in the model, we reduce the eight-dimensional network equations to a fully-explicit scalar Poincaré burst map. This map tracks the state of synaptic depression from one burst to the next and captures the complex bursting dynamics of the network. Fixed points of this map are associated with stable burst solutions of the full network model, and are created through fold bifurcations of maps. We derive conditions that predict the bifurcations between n : n and (n + 1) : (n + 1) solutions, producing a full bifurcation diagram of the burst cycle period. Predictions of the Poincaré map fit excellently with numerical simulations of the full network model and allow the study of parameter sensitivity for rhythm generation.


INTRODUCTION
Short-term synaptic plasticity may have a role in burst activity in central pattern generators (CPGs).Short-term synaptic depression is commonly found in neuronal networks involved in the generation of rhythmic movements, such as in the pyloric CPG of the spiny lobster [1,2], or in the lumbosacral cord of the chick embryo [3].Synaptic depression modulates the strength of synapses in response to changes to the presynaptic firing frequency.At a high neuronal firing frequency, depression weakens the strength of synapses and therefore reduces the magnitude of the postsynaptic response.At low firing frequency, it allows sufficient time for the synapse to recover from depression between spikes, leading to a stronger postsynaptic response.In reciprocal networks, synaptic depression has been shown to act as a "switch, " giving rise to a wide range of network dynamics such as synchronous and multi-stable rhythms, as well as fine tuning the frequency of network oscillations [4][5][6].
Brown [7] pioneered the idea that synaptic depression acts as a burst termination mechanism in CPGs composed of reciprocally inhibitory neurons and involved in rhythm generation of locomotion.When one side is firing during a burst the other, antagonistic side, is prevented from firing by synaptic inhibition.However, the weakening of inhibition as a result of synaptic depression eventually releases the antagonistic side so that it starts firing, terminating the burst on the side that had originally been firing.This rhythmogenesis hypothesis has been considered one of a handful of standard mechanisms for generating locomotion rhythms in vertebrates [8][9][10].It has been proposed as an explanation of the antiphase burst rhythm in struggling in Xenopus tadpoles [11].
Bose and Booth [6] investigated burst generation in a generic half-center CPG that consists of two identical, tonically active Morris-Lecar [12] neurons coupled through inhibitory depressing synapses.Numerical simulations showed that when the reciprocal synaptic conductance between the two neurons is varied, the network produces symmetric n : n anti-phase bursts, with stronger synaptic coupling leading to longer bursts.They used methods from geometric singular perturbation theory to separate the timescales of the fast membrane, and the slow synaptic dynamics of the network to derive one-dimensional conditions necessary for the existence of stable n : n solutions (for n ≤ 2).According to these conditions the type of firing pattern largely depends on the slow depression dynamics of the synapses between the two neurons, and can therefore be predicted by knowing the strengths of the synaptic conductances of the two synapses.Thus, the scalar conditions derived in Bose and Booth [6] provide a method to numerically identify the type of stable n : n pattern for any given value of the coupling strength and n ≤ 2. However, they do not predict the exact period of such solutions.Furthermore, while they provide good arguments for the validity of their reduction assumptions and the resulting scalar conditions, they do not verify them numerically.
Here we extend the previous analysis by providing a Poincaré map of the slow depression dynamics.This allows us not only to predict the types of stable n : n solutions the full network can produce, (for any n), but also to study how varying the coupling strength affects the period of such solutions.To do this, we build on, and numerically test, the assumptions on the fast-slow timescale disparity made in Bose and Booth [6].We reduce the two-cell model to a scalar Poincaré map that tracks the evolution of the depression from the beginning of one burst to the beginning of the next burst.Stable fixed points of our map are associated with stable n : n burst solutions.Our map construction is motivated by the burst length map of a T-type calcium current, utilized by Matveev et al. [13], which approximates the anti-phase bursting dynamics of a network of two coupled Morris-Lecar neurons.In contrast to our model, the network described in the Matveev et al. [13] paper does not contain short-term synaptic depression, and burst termination is instead accomplished through the dynamics of a slow T-type calcium current.
The Poincaré map derived here replicates the results from numerical simulations of the full two-cell ODE system: Given the strength of maximum conductance between the two neurons, fixed points of our map predict the type and period of n : n patterns, the switch between burst solutions of different periods, as well as the occurrence of co-existent solutions.In addition to proving the existence and stability of fixed points, our map shows that fixed points are created via a fold bifurcation of maps.Finally, we use our map to derive algebraic conditions that allow us to predict parameter values of the maximum conductance at which n : n solutions bifurcate to (n + 1) : (n + 1) solutions, and vice versa.Because our map is fully explicit, it lays the framework for studying the effects of other model parameters on network dynamics without the need to run expensive numerical integrations of the ODEs.
This paper is organized as follows.First, we introduce the network of two neurons, and describe the properties of single cell and synapse dynamics.We use numerical simulations of the network to provide an intuition for the range of possible burst dynamics the system can produce.Next, we state and justify the simplifying assumptions that are necessary for the map construction.Finally, we analytically derive the first return map of the depression variable as well as the conditions that are required for stable n : n solutions.We end this work with a discussion.

MATERIALS AND METHODS
We consider a pair of identical Morris-Lecar neurons [12], with parameters from Bose and Booth [6].The Morris-Lecar model is a set of two first-order differential equations that describe the membrane dynamics of a spiking neuron.The depolarisation is modeled by an instantaneous calcium current, and the hyperpolarisation by a slow potassium current and a leak current.The membrane potential v i and potassium activation w i of neuron i (i, j = 1, 2) is described by: Here v s is the inhibitory reversal potential, and ḡ and s j are the maximal synaptic conductance and the synaptic gating, respectively, constituting the total inhibitory conductance ḡs j from neuron j to neuron i. Function f (v i , w i ) describes the membrane currents of a single cell: The currents include a constant current I, and three ionic currents: an instantaneous calcium current, a potassium current, and a leak current, with respective reversal potentials v Ca , v K , and v L , as well as maximum conductances g Ca , g K , and g L .The function h(v i , w i ) models the kinetics of the potassium gating variable w i , and is given by The steady-state activation functions m ∞ and w ∞ as well as the default model parameters are described in the Supplementary Material 1.
The dynamics of the synaptic interactions between the neurons are governed by a synaptic gating variable s i and a depression variable Variable d i describes a firing rate dependent depletion mechanism that governs the amount of depression acting on the synapse.The model is agnostic with respect to the exact mechanism of this depletion, be it pre-or post-synaptic.When the voltage of cell i is above firing threshold (v i > v θ ), variable d i decays with time constant τ b , and recovers with time constant τ a when voltage is below firing threshold (v i < v θ ).Since the synaptic inhibition occurs on a much faster timescale than synaptic depression, we assume that s i is instantaneously reset to d i whenever v i increases above v θ , where it remains throughout v i > v θ .Whenever v i < v θ , the synaptic variable decays exponentially with time constant τ κ .The equations for the depression model are identical to the Bose and Booth [14] model.These equations are a mathematically tractable simplification of the established phenomenological depression model previously described by Tsodyks and Markram [15].When the total inhibitory conductance ḡs j is constant, the membrane dynamics are determined by the cubic v-nullcline v ∞ (v i ) and the sigmoid w-nullcline w ∞ (v i ), satisfying vi = 0 and ẇi = 0, respectively.In case of no inhibition (ḡ = 0), the two curves intersect near the local minimum of v ∞ to the left of v θ (commonly referred to as "left knee" of v ∞ ), creating an unstable fixed point p f with a surrounding stable limit cycle of period T = T a + T s (Figure 1A).Here T a is the amount of time the membrane potential spends above firing threshold (v i > v θ ), while T s is the time it spends below firing threshold (v i < v θ ).Trajectories along that limit cycle have the familiar shape of the action potential (Figure 1B).Applying a constant nonzero inhibition, e.g., by letting s j = 1 and ḡ > 0, moves the cubic v ∞ with the ensuing unstable fixed point down w ∞ in the (v i , w i ) -plane.When ḡ is large enough, the fixed point moves past the left knee and becomes stable via a subcritical Andoronov-Hopf bifurcation, attracting all previously periodic trajectories.
In the following section we will refer to the value of the total conductance ḡs j at the bifurcation point as g bif .
The two-cell network model is numerically integrated using an adaptive step-size integrator for stiff differential equations implemented with XPPAUT [16] and controlled through the Python packages SciPy [17] and PyXPP [18].The following mathematical analysis is performed on the equations of a single cell.Unless required for clarity, we will therefore omit the subscripts i, j from here on.

Anti-phase Burst Solutions
Short-term synaptic depression of inhibition in a half-center oscillator acts as a burst termination mechanism [7] and is known to produce n : n anti-phase burst solutions of varying period.Such n : n solutions consist of cells firing bursts of n spikes in alternation.Figure 2D shows the timecourse of a typical 4 : 4 burst.While one cell is firing a burst it provides an inhibitory conductance to the other cell, preventing it from firing.
Therefore, at any given moment one cell is spiking while the other is suppressed and does not spike.We will refer to the currently firing cell as "active" and we will call the suppressed cell "silent."Additionally, we will distinguish between two phases of a n : n solution: We will refer to the time interval when a cell is firing as the "active phase, " and we will call the remaining duration of a cycle, when a cell is not firing, the "silent phase." With each action potential of the active cell, short-term depression leads to a decrease of d, and consequently of s.If d depresses faster at spike time than it can recover in the inter-spike-intervals (ISIs), the total synaptic conductance ḡs will eventually become sufficiently small to allow for the silent cell to be released [19,20] and start firing, thus inhibiting the previously active cell.While a cell is silent its depression variable can recover.Once the silent cell becomes active again its synaptic inhibition will be sufficient to terminate the burst of the previously active cell and commence a new cycle.As previously demonstrated by Bose and Booth [6], in a two-cell reciprocally inhibitory network with synaptic depression the coupling strength ḡ determines the type of n : n solution.Increasing ḡ produces higher n : n burst solutions with more spikes per burst and a longer cycle period.Figure 2 shows numerically stable n : n solutions for varying values of ḡ.For small values of ḡ the network produces anti-phase spiking 1 : 1 solutions (Figure 2A).As ḡ is increased the network generates solutions of increasing n, that is 2 : 2 (Figure 2B), 3 : 3 (Figure 2C), and 4 : 4 (Figure 2D).When ḡ is sufficiently large (Figure 2E), one of the cells continuously spikes at its uncoupled period T while the other cell remains fully suppressed.Depending on the initial conditions either of the two cells can become the suppressed cell, which is why the suppressed solution is numerically bistable.
Branches of numerically stable n : n solutions and their associated limit cycle period for varying values of ḡ are depicted in Figure 3 (see Supplementary Material 2 for algorithm description).Not only do higher n : n solutions branches require stronger coupling ḡ, but also within n : n branches the period increases with ḡ.In line with Bose and Booth [6] we find small overlaps between solution branches indicating numerical bistability, for example such as between the 2 : 2 and 3 : 3 solution branches.Branches of higher n : n burst solutions occur on increasingly smaller intervals of ḡ, for instance is the ḡ interval of the 5 : 5 branch shorter than that of the 4 : 4 branch and so on.The interval between the 5 : 5 branch and the suppressed solution (region between dotted lines in Figure 3) not only contains even higher numerically stable n : n solutions, such as 11 : 11 bursts, but also other non-symmetric n : m solutions as well as irregular, non-periodic solutions.However, the analysis in the following sections will only be concerned with the numerically stable and symmetric n : n solutions.

Mathematical Analysis of Two-Cell Network
The goal of the following mathematical analysis is to reduce the complexity of the eight-dimensional system to a more tractable problem.As we will explain, we do this by approximating the full dynamics by a reduced system that describes the evolution of the depression variable d of either of the two cells.We will construct the solution of d in a piecewise manner from one spike to the next, first during the active phase, and then during the silent phase.This construction will require two assumptions about the membrane and synaptic dynamics.The first assumption states that during a burst the active cell fires at its uncoupled period T, which simplifies the construction of the solution of d.The second assumption states that once the inhibitory conductance acting on the silent cell drops below a critical threshold, the cell is immediately released and fires.The second assumption is necessary to predict the release time of the silent cell, which allows us to model the recovery of d during the silent phase.In other words, the second assumption requires that the release of the silent cell from inhibition depends only on the timecourse of the inhibition, and not on the membrane dynamics of the silent cell.The approximate validity of both assumptions can be observed in coupled relaxation-oscillator types of neurons such as the Morris-Lecar model we use, and will be numerically verified below.Both assumptions were first used in Bose and Booth [6] to derive algebraic conditions that guarantee the periodicity of the depression variable for different n : n solutions.However, here we will use these assumptions to construct a Poincaré map of d, which will provide a geometric intuition for the dynamics of the full two-cell network and its dependence on model parameters.
Our first assumption about the model states that the active cell fires at its uncoupled period T, that is, during the active phase of a burst we have ISI = T. Solution profiles in Figure 2 suggest that the ISIs are indeed approximately constant.Numerically computing ISIs for all stable n : n solutions in Figure 3 reveals that ISIs differ by at most 1 ms from the intrinsic firing period T ≈ 376 ms.Assuming ISI = T seems reasonable given that inhibition acting on the silent cell decays exponentially on a much shorter timescale τ κ than the duration of the ISI.Therefore, once the silent cell is released its trajectory quickly approaches the spiking limit cycle.Naturally the above assumption requires a sufficiently small τ κ , and fails when τ κ is large.In the Supplementary Material 3 we numerically explore how different values of τ κ affect the ISIs of the active cell.Finally, assuming ISI = T allows us to ignore the non-linear membrane dynamics during the active phase, and to construct the evolution of the synaptic variables iteratively from spike to spike.
Our second assumption states that the silent cell is released and spikes as soon as the total inhibitory conductance ḡs acting on it drops below some threshold value.We call this critical threshold value the "release conductance, " and define it as the value of ḡs at the time when the voltage of the silent cell first crosses the firing threshold v θ , that is when that cell is released and fires its first spike.Recall that when a cell is silent its v-and w-nullclines intersect at a stable fixed point and ḡs > g bif .A sufficient condition for the silent cell to be released is therefore ḡs < g bif .However, depending on the topology of the stable manifold, the (v, w)-trajectory of the silent cell can escape the stable fixed point and allow the cell to produce a spike for ḡs > g bif .In this case the value of the release conductance depends on the type of n : n solution and the coupling strength ḡ.For any stable n : n solution in Figure 3   conductance converges to some constant conductance value g ⋆ ≈ 0.0068 mS/cm 2 .Here g ⋆ is the value of ḡs at the end of a cycle of a suppressed solution, just before the active cell spikes.Using g ⋆ as a constant approximation for the release conductance will allow us to formulate a scalar condition that predicts the release time of the silent cell.Moreover, using g ⋆ is convenient because its exact value can be derived explicitly, as will be shown in the following section.
Assuming a constant release conductance for all n : n solutions will naturally introduce some error in the prediction of the release time of the silent cell.We can compute that error for any associated solution in Figure 4A by calculating the time interval between the first spike of the silent cell and the time when ḡs first crosses g ⋆ .We will call this time interval the "release delay." Figure 4B shows the numerically computed graph of such release delays.For n > 1 the absolute delays are smaller than 2 ms.Therefore, using as a constant release condition for all n : n solutions allows us to accurately predict the timing of the release of the silent cell.And to simplify the terminology, from now on we will refer to Equation ( 7) simply as the "release condition."In summary: We assume that the release condition is sufficient to predict when the silent cell is released.Due to the symmetry of n : n solutions the release occurs at exactly half the period of the full cycle.The release time therefore uniquely determines the type of n : n solution.Furthermore, computation of the release time does not depend on the membrane nor the synaptic dynamics of the silent cell.Instead, the solution of the synaptic variable s of the active cell is sufficient to predict when ḡs = g ⋆ is satisfied.Finally, the value of s at each spike time is determined by the evolution of the depression variable d of the active cell.Constructing a solution of d during the active phase of either cell will therefore uniquely determine the solution of the full eightdimensional network.However, finding the solution d requires us to know the initial value d(0) at the start of a cycle at t = 0.In the next section we will construct a scalar return map that tracks these initial values d(0) from cycle to cycle of stable n : n solutions.

Construction of the Scalar Poincaré Map
In this section we construct the scalar Poincaré map n : d ⋆ → d ⋆ .Here the discrete variable d ⋆ tracks the values of the continuous depression variable d at the beginning of each n : n burst.The map n therefore describes the evolution of d, of either of the two cells, from the beginning of one cycle to the beginning of the next cycle.To simplify the map construction we will assume that an active cell fires exactly n times before it becomes silent.We will construct n by evolving d first during the active phase and then during the silent phase of the n : n limit cycle.The terms "active" and "silent" phases will be defined in terms of the state of the depression variable.During the active phase the depression variable of the active cell both decays and recovers, while during the silent phase it only recovers.First, let us give explicit definitions of the active and silent phases of a burst.A schematic illustration of both phases is given in Figure 5.
Suppose that at t = 0 cell 1 becomes active with some initial d(0).Cell 1 then fires n spikes at the uncoupled period T = T a + T s .Let s(t) and d(t) be the corresponding solutions of the synaptic and depression variables of cell 1.After n spikes the total conductance ḡs(t) acting on the silent cell 2 has decayed sufficiently to satisfy the release condition (Equation 7).That is at some time t = (n − 1)T + T a + t, where t < T s will be determined below, we have ḡs(t) = g ⋆ [6].Cell 2 is then released and prevents cell 1 from further spiking.Once released, cell 2 also fires n spikes until cell 1 becomes active once again.Let P n denote the full cycle period of a n : n solution: We can now define the active and silent phases of cell 1 explicitly.
The active phase of a burst is the interval that lasts from the first spike time up until the beginning of the silent phase of the last spike, that is for time 0 < t < (n − 1)T + T a .During the active phase of cell 1, the silent cell 2 is inhibited sufficiently strong to prevent it from firing, hence ḡs > g ⋆ .The silent phase of cell 1 is the remaining duration of the cycle when the cell is not firing, that is for (n − 1)T + T a < t < P n .The silent phase lasts for (n − 1)T + T a + 2 t.
Note that only the silent phase depends on t, which will play a central role in the construction of n .From Equation (8) t can be computed as We can use Equation ( 9) and the numerically computed bifurcation diagram of the period for stable n : n solutions in Figure 3 to obtain the graph of t as a function of ḡ (Figure 6).Each continuous branch of t is monotonically increasing and corresponds to a n : n burst: Stronger coupling ḡ increases the total synaptic conductance ḡs that acts on the silent cell, thus delaying its release.It is easy to see that for any n-branch we have t < T s : Once t crosses T s , the active cell can "squeeze in" an additional spike and the solutions bifurcate into a (n+1) : (n+1) burst.
Distinguishing between the active and silent phases of a n : n cycle allows us to describe the dynamics of the depression variable d explicitly for each phase.As can be seen from Figure 5C, during the active phase d depresses when v > v θ and recovers when v < v θ .In contrast, during the silent phase d only recovers and does not depress.Given the initial d ⋆ = d(0) at the beginning of the cycle and the number of spikes in the active phase n, we can now construct the burst map n .The map is a composition of two maps.Map models the evolution of d in the active phase.F n takes an initial value d ⋆ and calculates t.Map models the recovery of d in the silent phase.Given some t map Q n computes d ⋆ at the start of the next cycle.
Our aim in the following analysis is to elucidate the properties of n and to understand the structure of its parameter space by exploring how the stable and unstable fixed points of n are created.To that effect it will be useful to include not only positive, but also negative values of d ⋆ to the domain of n .But it is important to add that values d ⋆ < 0 are biologically impossible as the depression variable models a finite pool of neurotransmitters, and therefore must be positive.Because n maps first from d ⋆ to t, and then back to d ⋆ , we will also consider negative values of t, interpreting them as n : n solutions with partially overlapping bursts.As will become evident, t < 0 is only a formal violation of the biological realism of the map n , as numerically stable n : n solutions of the full system of ODEs only exist for t > 0. We start the construction of n by first considering the active phase and building the map F n .At each spike time t k where d(t k ) = d k , variable d decays first for the duration of T a , as described by the solution to Equation (5).At t = t k + T a we have The depression variable then recovers for T s until t k+1 , where for 0 < t < T s : By substituting t = T s we can build a linear map that models the depression of d from spike time t k to the subsequent spike time t k+1 during the active phase: where to keep the notation simple we let Given constant T a and T s , the derived parameter λ determines how much the synapses depresses when v > v θ , while ρ determines how much it recovers when v < v θ .Since 0 < λ, ρ < 1, the map in Equation ( 15) is increasing and contracting, with a fixed point at where 0 < d s < 1.The value d s is the maximum depression value that can be observed in the suppressed solution where the active cell fires at its uncoupled period T (see Figure 2E).Using the release condition in Equation ( 7) allows us to derive the value of the minimum coupling strength that will produce the full suppressed solution, denoted as ḡs .Solving Equation (7) for s(t) Each continuous branch is associated with a stable n : n burst solution.
Increasing ḡ increases t until the solutions bifurcate at t ≈ T s .
with t = T s and setting the initial value s(0) = d s λ then gives us the aforementioned approximation of the release conductance g ⋆ : By substituting the definition of d s in Equation ( 18) and rearranging, we can also write ḡs as a function of λ and ρ: Note that the above dependence of ḡs on λ is linear and monotonically decreasing.Increasing λ reduces the strength of the depression of the active cell.This in turn allows the active cell to fully suppress the silent cell at smaller values of ḡ.Solving Equation (15) gives us the linear map δ n , that for some initial d ⋆ computes the depression at the nth spike time, that is d(t n ): Since λ < 1, function δ n is a linearly increasing function of d ⋆ with a fixed point at d s for all n.Having identified d after n spikes, we can now use the release condition ḡs = g ⋆ (Equation 7) to find t.At the last (nth) spike of the active phase at time t n = (n − 1)T the synapse variable s is set to the respective value of d(t n ) = δ n (d ⋆ ), and mirrors the value of d for the duration of T a .At the end of the active phase at time t n + T a variable d has decayed to δ n (d ⋆ )λ, therefore Finally s decays exponentially for t < T s .Solving (Equation 6) with initial condition s(0) = δ n (d ⋆ )λ yields: Substituting s( t) into s of the release condition (Equation 7) gives then Our assumption of the release condition guarantees that the silent cell 2 spikes and becomes active when ḡs−g ⋆ crosses zero.Solving (Equation 24) for t allows us to compute t as a function of d ⋆ , which defines the map F n : Figure 7A shows F n for various n, which is a strict monotonically increasing function of d ⋆ as well as ḡ.Larger values of d ⋆ and ḡ, respectively, cause stronger inhibition of the silent cell, and therefore prolong its release time and the associated t.Map F n is defined on d ⋆ > d a , where d a is a vertical asymptote found by solving δ n (d ⋆ ) = 0 in Equation ( 21) for d ⋆ , which yields We now turn to the construction of map Q n , which describes the recovery of the depression variable during the silent phase.As we have identified earlier, the recovery of d in the silent phase of a n : n solution starts at time t n + T a and lasts for the duration of (n − 1)T + T a + 2 t.Substituting that duration into the solution of d (Equation 5) with the initial condition d(0) = δ n (d ⋆ )λ yields the map Q n : We can find δ n (d ⋆ ), i.e., the value of d at the nth spike time, by rearranging the release condition in Equation ( 24): Map Q n is shown in Figure 7B for various values n.Note that Q n is monotonically increasing as larger values t imply a longer recovery time, and hence Q n grows without bound.All curves Q n intersect at some t = τ κ ln ḡ/g ⋆ where As we will show in the next section, all fixed points of the full map n occur for d ⋆ < 1.We will therefore restrict the domain of Q n to (−∞, τ κ ln ḡ/(g ⋆ ) ) and the codomain to (−∞, 1).Additionally, while values t > T will be helpful in exploring the geometry of n , recall from Figure 6 that in the flow system n : n solutions bifurcate into (n + 1) : (n + 1) solutions exactly when t = T s , and we will address this concern in the last part of our map analysis.
Having found F n and Q n , we can now construct the full map where we substituted τ = 2τ k /τ a .Recall that δ n (d ⋆ ) and g ⋆ are obtained from Equations ( 21) and (19), respectively.Since d is the slowest variable of the system and τ a ≫ τ κ , we will also assume τ < 1. Figure 8A depicts n for various n.Intersections of n with the diagonal are fixed points of the map. Figure 8B shows 2 with n = 2. Varying the synaptic strength ḡ moves the curves n up and down the (d ⋆ , n )-plane.For ḡ < 0.0015 mS/cm 2 map 2 has no fixed points.As ḡ is increased to ḡ ≈ 0.0015 mS/cm 2 , curve 2 coalesces with the diagonal tangentially.When ḡ > 0.0015 mS/cm 2 , a pair of fixed points emerge, one stable and one unstable fixed point, indicating the occurrence of a fold bifurcation of maps.
n is monotonically increasing with respect to ḡ and also d ⋆ : The monotonicity of n w.r.t.ḡ is evident from Equation ( 30), while the monotonicity w.r.t.d ⋆ follows from the monotonicity of both Q n and F n .In the following sections we will heavily rely on this monotonicity property of n .Just as F n , curves n spawn at the asymptote d a (Equation 26), and because fixed points of n lie in (d a , 1).

Existence and Stability of Fixed Points
We introduce the fixed point notation . The existence of fixed points d ⋆ f for ḡ sufficiently large can be shown from the strict monotonicity of n with respect to ḡ and d ⋆ (Equations 32, 31), as well as the fact that the slope of n is monotonically decreasing, In the limit d ⋆ → d a the value of n decreases without bound for any ḡ > 0. In the limit ḡ → 0, n also decreases without bound, but as ḡ → ∞ values of n approach 1.It follows from Equation (31) and the intermediate value theorem that for some ḡ large enough n intersects the diagonal.Moreover, because n and its slope are monotonic with respect to d ⋆ , there exists some critical fixed point (d ⋆ b , ḡb ) where n aligns with the diagonal tangentially with (36)

Fold Bifurcations of Maps
Fixed points of n satisfy the fixed point equation As we have already shown, for ḡ > ḡb (n) solutions to Equation (37) exist in pairs of stable and unstable fixed points.Solving (Equation 37) explicitly for d ⋆ is not trivial, but solving for ḡ is straightforward and given by ḡ = G n (d ⋆ ), where is defined for d ⋆ < 1 and δ n (d ⋆ ) > 0. Plotting d ⋆ against ḡ gives the fixed point curves, which are shown in Figure 9A.Note the typical quadratic shape of a fold bifurcation of maps.It is also evident that the fold bifurcations occur for increasingly smaller ḡ as n is increased.Moreover, the graph suggests that for n > 1 unstable fixed points have negative values of d ⋆ .Equation (38) also allows us to find the critical fixed point connected with the fold bifurcation, namely Function G n is strictly monotonic on the respective intervals of d ⋆ f that correspond to the stable and unstable fixed points, that is which allows us to express the stable and unstable fixed points as the inverse of G n on their respective intervals of d ⋆ f .Because we are primarily interested in the stable fixed points Function φ n (ḡ) is also monotonic, and is therefore straightforward to compute numerically.We use the Python package Pynverse [21] for that purpose.
Having found the stable fixed points d ⋆ f as a function of ḡ, we can now compute the associated cycle period.Recall that the period is given by Equation ( 8), which can be written as a function of ḡ: where map F n (Equation 25) calculates t given a stable fixed point d ⋆ f = φ n (ḡ). Figure 9B shows the period P n (ḡ) computed from Equation (44) versus the cycle period of stable n : n solutions, computed from numerically integrating the full system of ODEs.The overlap between blue and orange curves suggests that stable fixed points of n accurately predict the cycle period of stable solutions of the flow system.
It is evident from Figure 9A that φ n is strictly increasing with ḡ.This property follows directly from the quadratic normal form of the fold bifurcation, but can also be shown using implicit differentiation and the fixed point equation n [φ n (ḡ), ḡ] = 0 in Equation (37).For The inequality follows from ∂ n /∂ ḡ > 0 and the fact that 45) allows us to explain why the period P n increases with ḡ, as seen in Figure 9B.Differentiating P n gives: where the partial derivatives of F n (d ⋆ f , ḡ) are: Fold bifurcation diagrams of stable (continuous curves) and unstable (dotted curves) fixed points of n for varying n. (B) Cycle periods computed from stable fixed points of n (blue), and the corresponding periods from stable n : n solutions acquired via numerical integration of the system of ODEs (orange).
Equation ( 46) and ( 47) have an intuitive biological interpretation: Increasing the coupling strength between the neurons leads to overall stronger inhibition of the silent cell, which delays its release and leads to a longer cycle period.The latter allows more time for the synapse to depress in the active phase and recover in the silent phase, resulting in overall larger values of d ⋆ f , that is weaker depression at the burst onset.
While fixed points of our Poincaré map predict the cycle period of the flow system excellently, its construction relies on the strong assumption that the active phase contains exactly n spikes.As is evident from Figure 9B this assumption is clearly violated in the flow system, as stable n : n bursts exists only on certain parameter intervals of ḡ.The multi-stability of fixed points of maps n in Figure 9B does therefore not imply a similar multi-stability of the flow system.In the last sub-section we will analyze the mechanisms that guide how the stable n : n are created and destroyed, and use our previous analysis to derive the corresponding parameter intervals of ḡ where such solutions exist.

Stable Solution Branch Borders
Let ḡL (n) and ḡR (n) denote the left and right parameter borders on ḡ where stable n : n solutions exist.That is, as ḡ is increased stable n : n solutions are created at ḡL (n) and destroyed at ḡR (n).When ḡ is reduced beyond ḡL (n), n : n solutions bifurcate into (n − 1) : (n − 1) solutions, while when ḡ is increased beyond ḡR (n), n : n solutions bifurcate into (n + 1) : (n + 1) solutions.Let us briefly recap our observations regarding ḡL (n) and ḡR (n) from the numerical bifurcation diagram in Figure 9B.For n > 1 there are the following relations: ḡL (n) < ḡR (n + 1) and ḡR (n) < ḡL (n + 1), ( 50) Equations (49, 50) are self-explanatory.Equation (51) formally describes occurrence of co-existence between stable n : n and (n+ 1) : (n + 1) solutions.Equation ( 52) implies that the parameter interval on ḡ of n : n solutions decreases with n, in other words, bursts with more spikes occur on increasingly smaller intervals of the coupling strength.All of the above relations are reminiscent of the bifurcation scenario of type period increment with coexistent attractors, first described for piecewise-linear scalar maps with a single discontinuity by Avrutin and colleagues [e.g., see [22][23][24].While our maps n are fully continuous, the above observation suggests that a different piecewise-linear scalar map that captures such period increment dynamics of the full system exist.We will explore what such a map might look like in the discussion.
Let us now find algebraic equations that will allow us to calculate the critical parameters ḡL (n) and ḡR (n) associated with the left and right n : n branch borders.Recall that the period P n derived from the fixed points of n is an increasing function of ḡ (Equation 46).That is, as the coupling strength increases, it takes longer for the total synaptic conductance to fall below the value of the release conductance, which delays the release of the silent cell, and t becomes larger.When t > T s , the active cell can produce another spike and the solution bifurcates into a (n + 1) : (n + 1) solution.Note, however, that at ḡL (n) the bifurcation into a (n − 1) : (n − 1) does not occur at t = 0.Here the mechanism is different: A sufficient reduction of ḡ causes the total synaptic conductance to drop below the release conductance in the previous ISI, which allows the silent cell to be released one spike earlier.
Using the above reasoning we can now formulate the conditions for both bifurcations at ḡL (n) and ḡR (n).As in the previous sections, we will only restrict ourselves to the analysis of the stable fixed points given implicitly by d ⋆ f = φ n (ḡ) (Equation 43).At the right bifurcation border ḡR (n) we have t = T s , and after substituting our F n map (Equation 25) this translates into which lets us define a function whose root is the desired right bifurcation border ḡR (n).In case of the left bifurcation border at ḡL (n), the release condition is satisfied just before the active cell has produced its nth spike, where total synaptic conductance is given by which can be rewritten as a function whose root is ḡL (n).Both R n and L n are increasing with respect to ḡ, which makes finding their roots numerically straightforward.
Figure 10 shows the period P n (ḡ) as predicted by the fixed points of n (Equation 44) plotted on their respective intervals ḡ ∈ [ḡ L (n), ḡR (n)] (blue), as well as the cycle period acquired from numerical integration of the full system of ODEs (orange).
Here g L (n) and ḡR (n) were computed from Equations ( 56) and (54), respectively.Note that the width of n : n branches decreases with n, which confirms the inequality in Equation (52).That is, bursts with more spikes occur on increasingly smaller intervals of ḡ, which can be interpreted as a lost of robustness with respect to the coupling strength of long-cyclic solutions.We also note the occurrence of bistability between pairs of n : n and (n + 1) : (n + 1) branches, which also confirms our initial observation in Equation (51).As previously observed in Figure 9B our maps prediction of the cycle period is accurate.However, the mismatch in the left and right branch borders is significant.This mismatch might be due to the millisecond release delay error (Figure 4B) induced by our assumption of a constant release conductance for all n : n solutions (see Equation 7).Another explanation for the border mismatch could be that our assumptions on the time scales of (v, w) vs s-and d-dynamics do not hold near the stability borders, and that they can only be captured by more complex approximations.Nevertheless, our map allows approximate extrapolation of the cycle period and the respective bifurcation borders where numerical integration of the ODEs would require a very small time step.

DISCUSSION
Synaptic depression of inhibition is believed to play an important role in the generation of rhythmic activity involved in many motor rhythms such as in leech swimming [25] and leech heart beat [26], and in the lobster pyloric system [1,2].In inhibitory half-center CPGs, such as believed to be found in the struggling network of Xenopus tapdoles, synaptic depression can act as a burst termination mechanism, enabling the alternation of bursting between the two sides of the CPG [11].Modeling can shed light on the underlying mathematical principles that enable the generation of such anti-phase bursts, and help identify the components that control this rhythm allowing it to switch between different patterns.
To study the mechanisms of burst generation in halfcenter CPGs we have analyzed a neuronal model network that consists of a pair of inhibitory neurons that undergo a frequency dependent synaptic depression.When the strength of synaptic inhibition between the neurons is varied, such a simple network can display a range of different n : n burst patterns.Using the timescale disparity between neuronal and synaptic dynamics, we have reduced the network model of eight ODEs to a scalar first return map n of the slow depression variable d.This map n is a composition of two maps, F n and Q n , that model the evolution of the depression during the active and silent phases of n : n solutions respectively.Both F n and Q n maps are constructed by using the dynamics of a single uncoupled neuron.Fixed points of n are created in pairs through a fold bifurcation of maps, where the stable fixed point correspond to stable n : n burst solutions of the full two-cell system of ODEs.The results from our onedimensional map match excellently with numerical simulation of the full network.Our results are also in line with Brown's [7] rhythmogenesis hypothesis, namely that synaptic depression of inhibition is a mechanism by which anti-phase bursting may arise.
We have studied n : n solutions assuming that the synaptic coupling ḡ between the two cells is symmetrical.However, Bose and Booth [6] have shown that asymmetrical coupling (ḡ 1 , ḡ2 ) can result in network solutions of type m : n, where one cell fires m spikes, while the other n spikes.It is conceivable that our map construction can be extended to also capture such m : n solutions.Remember, in the case of symmetrical coupling with n : n solutions, the timecourse of the depression variables d 1 and d 2 were in anti-phase, and it was therefore sufficient to track only one of the two variables.To capture the full network dynamics in case of asymmetrical coupling one would also have to account for burst patterns of type m : n, where the solutions of the depression variables d 1 and d 2 are not simply time-shifted versions of each other.To do that, one could track the state of both variables by constructing a two-dimensional Poincaré map (d 1 , d 2 ).While geometrical interpretation of two-dimensional maps remains challenging, there exist a number of recent studies which have employed novel geometrical analysis methods to understand the dynamics of two-dimensional maps of small neuronal networks [27][28][29].Generally speaking, our map construction approach is applicable to any small network, even with more than two neurons.As long as the network dynamics occur on separable timescales the main challenges to the map construction lie in identifying the slowest variables, and finding an appropriate, simplified description of their respective timecourses.In theory, the reduction approach can be also applied to neuronal systems with more than two timescales [e.g., see 30].
In tadpoles, struggling is believed to be initiated by an increase in the firing frequency of reciprocally inhibitory commisural interneurons, which has been hypothesized to lead to stronger synaptic depression of inhibition and result in the iconic antiphase bursting [11].It would therefore be interesting to study how varying the cell intrinsic firing period T could affect the network rhythm.While we have laid out the framework to perform such an investigation, due to the choice of neural model we have avoided varying T. Recall that T is a derived parameter in the Morris and Lecar [12] model, and can therefore not be varied in isolation of other model parameters.This makes verifying any analytical results from our map analysis via numerical integration of the ODEs difficult.A more abstract model such as the quadratic integrate-and-fire model [31] allows varying T independently of other model parameters, and could be more fitting for such an investigation.
Our simulations of the network showed that n : n solutions lose robustness as their period is increased.That is, solutions with a larger cycle period occur on increasingly smaller intervals of the coupling strength.We were able to replicate this finding by numerically finding the respective left and right borders of stable n : n branches of fixed points of n , and showing that the distance between these borders shrinks with n.We have also noted the resemblance of our bifurcation diagram to one where such n : n branches are created via the bifurcation scenario of type period-increment with co-existent attractors, first described for scalar linear maps with a discontinuity [24,32].It is worthwhile noting that the bifurcations of piecewise linear maps studied by Avrutin et al. [32] result from a "reinjection" mechanism [33].
Here the orbit of a map performs multiple iterations on one side of the discontinuity, before jumping to the other side and being reinjected back into the initial side of the discontinuity.The stark difference of such a map to our map is that reinjection allows a single scalar map to produce periodic solutions of varying periods.In contrast, we rely on n different maps n to describe the burst dynamics without explicitly capturing the period increment dynamics.It is therefore conceivable that despite the complexity and non-linearity of the dynamics of our two-cell network, a single piecewise-linear map might be already sufficient to capture the mechanisms that shape the parameter space of the full system.In their discussion, Bose and Booth [6] briefly outline ideas about how such a linear map could be constructed.
In addition to stable n : n solutions, the numerical continuation by Bose and Booth [6] also revealed branches of unstable n : n solutions.While we have identified fold bifurcations of our burst map, we have not found corresponding bifurcations of the flow ODE system, and have generally ignored the significance of unstable map fixed points.However, the quadratic nature of the period bifurcation curve is reminiscent of a saddle-node on an invariant circle (SNIC) bifurcation, where the oscillation period lengthens and finally becomes infinite as a limit cycle coalesces with a saddle point.SNIC bifurcations have been studied in great detail [e.g., 34], and a next step would be to provide a rigorous explanation of not only the map dynamics, but also of the flow dynamics of the ODE system.
We have shown that when the strength of the maximum synaptic conductance is varied, synaptic depression of inhibition can enable our two-cell network to produce burst solutions of different periods.This result is in line with the idea that one role of synaptic depression in the nervous system may be to allow a finite size neuronal network to participate in different tasks by producing a large number of rhythms [6,11,35].To change from one rhythm to another would only require a reconfiguration of the network through changes in synaptic coupling strength.Thus short-term synaptic depression of inhibition may provide means for a network to adapt to environmental challenges without changing its topology, that is without the introduction or removal of neurons.

FIGURE 1 |
FIGURE 1 | Periodic solution of ML model neuron.(A) Projection of limit cycle onto (v, w)-phase plane with v-nullcline (blue, v ∞ ) and w-nullcline (orange, w ∞ ).Unstable fixed point p f is indicated by an orange dot; firing threshold v θ is denoted by a dashed line.(B) Corresponding voltage trace v(t) of an action potential.
we can compute an associated release conductance numerically by recording the value of ḡs at the time of the first spike of the silent cell.Such values of the release conductance are shown in Figure 4A, and the graph suggests that as n increases, the value of the release

FIGURE 3 |
FIGURE 3 | Numerically computed bifurcation diagram of the cycle period of stable n : n solutions for increasing coupling strength ḡ.Regions of bistability are indicated by light blue vertical stripes.Dashed lines show the interval between the 5 : 5 and the suppressed solution, where higher period n : n solutions occur on increasingly smaller intervals of ḡ.

FIGURE 4 |
FIGURE 4 | Numerically computed values of the release conductance (A) and release delay (B) for various n : n solutions and values ḡ.The dashed line indicates the analytical approximation of the release conductance by g ⋆ .

FIGURE 5 |
FIGURE 5 | Schematic diagram of the active and silent phases for a 3 : 3 solution.(A) Membrane potentials of cell 1 (v 1 blue) and cell 2 (v 2 orange).Gray patches depict t intervals.(B) Total synaptic conductance of cell 1 ( ḡs 1 ) as it crosses the release conductance g ⋆ .(C) Solution d 1 (t) of depression variable of cell 1, during active (blue) and silent phases (orange).

FIGURE 6 |
FIGURE 6 | Numerically computed bifurcation diagram of t for varying ḡ.Each continuous branch is associated with a stable n : n burst solution.Increasing ḡ increases t until the solutions bifurcate at t ≈ T s .

FIGURE 7 |FIGURE 8 |
FIGURE 7 | Maps F n (A) and Q n (B) for ḡ = 0.5 mS/cm 2 and n = 1, 2, 3, 4. Curves F n intersect at d ⋆ = d s which is indicated by a dashed vertical line.Curves Q n intersect at t = τ κ ln ḡ/g ⋆ .

FIGURE 10 |
FIGURE 10 | Bifurcation diagram of the period of stable n : n solutions computed analytically from fixed points of n , plotted on the respective intervals of ḡ ∈ [ ḡL (n), ḡR (n)] (blue), and computed from numerical integrations of the ODEs (orange).