A scalar Poincaré map for struggling in Xenopus tadpoles

Synaptic plasticity is widely found in many areas of the central nervous system. In particular, it is believed that synaptic depression can act as a mechanism to allow simple networks to generate a range of different firing patterns. The simplicity of the locomotor circuit in hatchling Xenopus tadpoles provides an excellent place to understand such basic neuronal mechanisms. Depending on the nature of the external stimulus, tadpoles can generate two types of behaviours: swimming when touched and slower, stronger struggling movements when held. Struggling is associated with rhythmic bends of the body and is accompanied by anti-phase bursts in neurons on each side of the spinal cord. Bursting in struggling is thought to be governed by a short-term synaptic depression of inhibition. To better understand burst generation in struggling, we study a minimal network of two neurons coupled through depressing inhibitory synapses. Depending on the strength of the synaptic conductance between the two neurons, such a 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. Using a fast/slow analysis, we reduce the multidimensional network equations to a 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 prove that the map has an infinite number of stable fixed points for a finite coupling strength interval, suggesting that the full two-cell network also can produce n – n bursts for arbitrarily large n. Our findings further support the hypothesis that synaptic depression can enrich the variety of activity patterns a neuronal network generates.

burst to the next, and whose fixed points correspond to stable burst solu-71 tions of the two-cell network. Our map construction is motivated by the burst 72 length map of a T-type calcium current, utilised in Matveev et al. (2007), 73 which approximates the anti-phase bursting dynamics of a network of two 74 coupled Morris-Lecar ( Morris and Lecar, 1981) neurons. This network does 75 not contain short-term synaptic depression, and burst termination is instead 76 accomplished through the dynamics of a slow T-type calcium current which 77 makes the single cells bursty, even in absence of network interactions. 78 This paper is organised as follows. First, we introduce the network of two 79 neurons, and describe the properties of the single cell and synapse dynamics. 80 We use numerical simulations of the network to provide an intuition for the 81 range of possible burst dynamics the system can produce. Next, we state  potential v i and potassium activation w i of neuron i (i, j = 1, 2) is described (1) (2) Here v s is the inhibitory reversal potential, andḡ and s j are the maxi-100 mal synaptic conductance and the synaptic gating, respectively, constituting 101 the total inhibitory conductanceḡs j from neuron j to neuron i. Function f (v i , w i ) describes the membrane currents of a cell that include a constant 103 current I, and three ionic currents: an instantaneous calcium current, a 104 potassium current, and a leak current, with respective reversal potentials 105 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 107 w i , and is given by The asymptotic functions m ∞ and w ∞ for the calcium and potassium con-109 ductances, respectively, are given by All model parameters are given in Table 1. 111 Our mathematical analysis is performed on the equations of a single cell, 112 which is why we will omit the subscripts i, j from here on. When the cells 113 are uncoupled (ḡ = 0), the membrane dynamics are determined by the cubic and the sigmoid w-nullcline w ∞ (v), satisfyingv = 0 andẇ = 0, respectively.

116
The two curves intersect along the middle branch of v ∞ , creating an unstable 117 fixed point p f with a surrounding stable limit cycle of period T ( Figure   118 1A). Trajectories along that limit cycle have the familiar shape of the action 119 potential ( Figure 1B). The trajectory of an action potential can be dissected  and all periodic trajectories relax into an equilibrium (see Figure 2). There  Thus, when a constant inhibitory conductance is applied, g acts as a release 137 conductance: a cell is prevented from tonic spiking whenḡs > g and it will 138 fire spikes whenḡs < g (Bose and Booth, 2011). The dynamics of the synaptic interactions between the neurons is gov-140 erned by a synaptic gating variable s, and a depression variable d.
Variable d describes a firing rate dependent depletion mechanism that governs 142 the amount of depression acting on the synapse; our model is agnostic with 143 respect to the exact mechanism of this depletion, be it pre-or post-synaptic.

144
At spike time when the voltage v crosses a synaptic threshold v θ , the synaptic 145 gating variable s is reset to the current value of the depression variable d, 146 and d is scaled down by the depression strength parameter λ (0 < λ < 1).

147
This results in the following discontinuous reset rule on spike time: that is while v > v θ for a decay time constant τ β . To simplify the mathematical 155 analysis we have reduced these dynamics to a single-reset event at spike time 156 according to (11). It is easy to see that this simplification is natural given become evident from our results that this simplification does not qualitatively change the dynamics of the network compared to the original equations with bursts of n spikes in alternation (see Figure 3A). While one cell is firing a 166 burst it provides an inhibitory conductance to the other cell thus preventing 167 it from firing. Therefore, at any given moment one cell is spiking while the 168 other is inhibited. Consistent with Bose and Booth (2011) we will refer to 169 the currently firing cell as "free" and we will call the inhibited cell "quiet".

170
Additionally, we will distinguish between two phases of a n − n solution: We 171 will refer to the burst duration of a cell as the "free phase", and we will 172 call the remaining duration of a cycle, when a cell is not spiking, the "quiet 173 phase".

174
With each action potential of the free cell, short-term depression leads to 175 a step-wise decrease of d, and consequently of s. If d depresses faster at spike 176 time than it can recover in the inter-spike-intervals (ISIs), the total synaptic 177 conductanceḡs will eventually become sufficiently small to allow the quiet 178 cell to be released and start firing, thus inhibiting the previously free cell.

179
While a cell is quiet its associated depression variable can recover. Once the 180 quiet cell is released its inhibition will be sufficient to terminate the burst of 181 the previously free cell and commence a new cycle. Figure       The goal of the following mathematical analysis will be to reduce the 205 complexity of the eight-dimensional system to some easily tractable quantity.

206
As we will see later this quantity is the value of the depression variable d of 207 one of the two cells. We will then construct the evolution of d for the free and terms of d will allow us to explore the parameter space and understand how  The following mathematical analysis is greatly simplified by assuming 219 that the ISIs during the free phase are of constant duration T . That is to 220 say, once a cell becomes free, it fires at its uncoupled period T . This as-221 sumption seems reasonable given that (1) the inhibitory conductance decays 222 exponentially to zero during the quiet phase of a cycle, and (2) the fact that 223 many biophysical neuron models, such as the Morris-Lecar model, will fire 224 at approximately their uncoupled period T once released from inhibition.

225
It is evident from the burst solutions presented in Figure 4 that the re-226 spective ISIs in the free phase are nearly constant. We further verify this 227 observation by computing the associated ISIs of stable solutions found in 228 Figure 5A. Figure 5B illustrates the bifurcation diagram of ISIs for varying 229 ḡ. Naturally, burst solutions with n > 1 contain two types of ISIs: Long ISIs 230 during the quiet phase, and short ISIs in the free phase. Note that in the 231 latter ISIs are nearly constant, ISI ≈ T , affirming the above assumption.

232
The assumption ISI = T allows us to ignore the membrane dynamics of the 233 free cell, and construct the evolution of its synaptic variables iteratively from 234 spike to spike. By tracking the exact values of d and s of the free cell we 235 will determine when the quiet cell is released and becomes free, terminating 236 the burst of the previously free cell. In the next section we will define such a 237 "release condition" by exploiting the discrepancy in timescales between the 238 fast voltage dynamics, and the slow potassium and synaptic dynamics. timescales (Fenichel, 1979). Consider equations (1) and (2) for an uncoupled 246 Morris-Lecar neuron withḡ = 0. By rescaling f =f / we obtain the generic 247 equations of a fast/slow system: By analysing these equations in the limit → 0, we will find an explicit 249 description of the action potential that is of a lower dimension than the 250 original system. Consider two timescales: a slow timescale t that describes 251 the evolution of v and w in the silent and active phases, and a fast timescale 252 τ that is used to describe the jumps up and down. We initially inspect the 253 fast timescale τ : By first substituting t = τ , and then setting = 0 we 254 acquire the singularly perturbed system of equations Here the prime symbol " " represents the derivative with respect to τ . The while for the jump down w = w R . In the singular limit we can therefore 261 simplify the jump dynamics to the scalar equation: Now, consider the slow timescale t by setting = 0: The first equation describes the v-nullcline and the second equation models and v R (w) describe the v-nullcline on the left and right branches, that is for 269 v < v L and v > v R respectively. Then the slow timescale dynamics are given In sum, geometric singular perturbation theory allows us to reduce the 272 two-dimensional uncoupled Morris-Lecar model to two scalar differential equa-273 tions for the fast and slow phases of an action potential respectively, at least 274 in the limit → 0. In the above analysis we assumed uncoupled cells, that is 275ḡ = 0. In the following section we will extend the analysis to include synaptic 276 interaction in the network whenḡ > 0. Our assumption for the free cell states that it fires at its uncoupled period 279 (ISI = T ). We can therefore disregard weak inhibition (ḡs < g ) in our can be extended to include these piecewise-continuous synaptic dynamics.

288
Our singular equations for the fast system predict that v will be at its The fast time-scale reduction then guarantees that on the slow timescale 301 the trajectory of the quiet cell will jump up instantaneously, becoming free 302 and firing an action potential. As will become evident in the next section, 303ẇ ṡ indeed holds for relaxation oscillator type models and fast inhibitory 304 synapses, such as our configuration of the Morris-Lecar model.

305
In summary: The release condition is sufficient to predict when the quiet cell is released. Due to the symmetry In the next section we will construct a scalar return map that tracks these 318 initial values d(0) from cycle to cycle of stable n − n solutions. the cycle when cell 2 is bursting, that is when (n−1)T < t < 2(∆t+(n−1)T ).

347
Note that only the quiet phase depends on ∆t.

348
The quantity ∆t will play a central role in the construction of Π n . From We can use (24) and the numerically computed bifurcation diagram of P 352 in Figure 5A to obtain the graph of ∆t as a function ofḡ (Figure 7). The Distinguishing between the free and quiet phases of a cycle allows us to 359 describe the dynamics of the depression variable d explicitly for each phase.

360
As can be seen from Figure 6C, during the free phase d depresses at spike  the reset d(t k ) ← λd(t − k ) (equation (11)). For the duration of the subsequent 371 ISI (t k < t < T + t k ), the recovery of d is described by the solution to (9): By substituting t = T we can build a linear map that models the depression 373 of d from spike time t k to the subsequent spike time t k+1 during the free 374 phase: where to keep the notation simple we let ρ = e −T /τ d in (26). Since 0 < λ, ρ < 376 1, the map (26) is increasing and contracting, with a fixed point at where 0 < d sup < 1. As we will see later, d sup is associated with fixed points 378 of the burst map Π n corresponding to the suppressed solution (see Figure   379 4E). Solving (26) gives us the linear map δ n : d 1 → d n , that for some initial 380 d 1 = d(0) at t = 0 computes d n at the nth spike time t n : where 382 α(n) = (λρ) n−1 , Note that 0 < α, β < 1. Figure 8 shows δ n (d 1 ) for varying n. Since λ < 1, 383 function δ n is a linearly increasing function of d 1 , decreasing with n. Natu-384 rally, δ n has a fixed point at d sup for all n. Having identified d after n spikes, 385 we can now use the release conditionḡs = g (equation (22)) to find ∆t. At 386 the last spike of the free phase at time t n = (n − 1)T , synapse variable s is 387 reset according to equation (10): Here d(t − n ) = δ n (d(0)), and t − n is the time just before the nth spike. After where ∆t < T . Now, our fast-slow analysis guarantees that the quiet cell 392 2 will spike and become free whenḡs − g crosses zero. We can therefore 393 formulate the release condition as: Finally, solving for ∆t and substituting δ n (d k ) from the equation (28)    To elucidate the properties of F n we will consider negative values of d.

398
However, it is important to note that values d < 0 are biologically impossible 399 as d models a finite pool of neurotransmitter, and therefore must be positive.  Here d n is the value of the depression variable at the last spike t n of the free 426 phase. Given ∆t, we solve the release condition (22) for d n : Map Q n is shown in Figure 10 for various values n. Map Q n (∆t) is monoton-428 ically increasing as larger values ∆t imply a longer recovery time, and hence 429 Q n grows without bound. All curves Q n intersect at some ∆t = τ s ln ḡ g λ 430 where 431 Q n τ s ln ḡ g λ = 1.
As we will show in the next section, all fixed points of the full map Π n occur 432 for d k < 1. We will therefore restrict the domain of Q n to (−∞, ln (ḡ g λ )τ s ) 433 and the codomain to (−∞, 1). Having found F n and Q n , we can now construct the full map Π n (d k ) = 436 Q n [F n (d k )]:

Composition Map
where we substituted τ = 2τ s /τ d . Since d is the slowest variable and τ d > τ s , 438 we assume τ < 1. Figure 11A depicts Π n for various n. Just as F n , map 439 Π n is defined on (d a (n), 1). Intersections of Π n with the diagonal at some d Then the derivative with respect toḡ is calculated as We can see that the sign of DḡΦ n is determined by the sign of 1 − λδ n (d ).

451
Solving 1 − λδ n (d ) = 0 for d by applying the definition of δ n in equation Simple calculations reveal that d φ (n) is strictly increasing with n: Since d φ (n = 1) = λ −1 > 1 for λ ∈ (0, 1), we can conclude that The graph of the inverse ofḡ(d ) is shown in Figure 12 (orange line).

460
Note the typical quadratic shape of a fold bifurcation. It is also evident

Fold Bifurcations of Π n 471
The derivative of Π n (d k ) is given by It is straightforward to show that Π n is strictly monotonous with DΠ n > 0 473 and D 2 Π n < 0. In the limit d k → d a the slope DΠ n is unbounded, and since In the limit n → ∞ we have Since Π n (d sup ) is monotonically increasing with n, it follows that for any 495ḡ <ḡ sup there exists some N > 0 such that Π n (d sup ) > d sup for any n > N .

496
The monotonicity of Π n and DΠ n with respect to d k then implies that such 497 maps Π n have a pair of fixed points on (d a , 1). In other words, for any finite 498 interval ofḡ on (0,ḡ sup ) there exists an infinite number of maps Π n that have 499 a stable fixed point and are associated with a burst limit cycle of the full 500 network. In other One example of a long-cyclic 11 − 11 burst can be seen 501 in Figure 13. Note that while the map predicts stable n − n solutions for 502 anyḡ <ḡ sup sufficiently large, these solutions occur on increasingly small 503 parameter intervals ofḡ, which makes identifying them in the ODE system 504 numerically challenging.

505
Deriving the exactḡ for each n associated with the fold bifurcations is 506 cumbersome. However, in the limit λ → 1 the equations (48) and (42) be-507 come more manageable. To provide a better intuition for the map we will 508 therefore demonstrate the fold bifurcations in the theoretical case of total ab-509 sence of depression, that is when λ = 1. This scenario can be associated with 510 a full suppressed solution. In this case Π n | λ=1 simply models the exponential 511 recovery of d according to equation (9) for 0 < t < 2(∆t + (n − 1)T ) = P : Equation (51) has a fixed point at d = 1 for all n > 0, which allows us to 513 find ∆t = F n (1). We obtain: Parameter valuesḡ that are associated with the fold bifurcations satisfy 515 equation DΠ n (d = 1)| λ=1 = 1, solving forḡ gives whereḡ f (n) denotes the fold bifurcation value ofḡ as a function of n. Figure   517 14 depicts the values ofḡ f (n) for increasing n in the limit λ → 1 (orange), 518 as well whenḡ f (n) is numerically computed for λ = 0.5 (blue) using a con- holds, and cell 1 can produce one more spike. In this case, L will be at least 552 of length T and we set: We can find the value of the depression variable at the next spike of cell 1 at 554 time t = T by evolving d according to equation (25): Now we are in a position to repeat the preceding procedure and test once 556 more whether the quiet cell will remain quiet for a duration T . Assume that 557 after n spikes d will be sufficiently depressed to release the quiet cell before 558 the subsequent spike of the free cell. The exact time of that release will be 559 then given by the previously accumulated burst duration L and the inter- This iterative procedure allows us to define F recursively. Given some initial Note that as we update the input of the recursive call to F with we also implicitly track the evolution of d. Figure 15 shows the partially Since n is unknown, we can use the remainder to find ∆t: Then map Q is defined by: Here d n is the value of the depression variable at the last spike of the free 579 phase. We can find d n as in equation (40), that is Figure 16 shows the partially continuous map Q : L → d k+1 and the con- of Π n regarding monotony, slope, and dependence onḡ also hold for Π on 595 the respective continuous intervals. Most importantly, it follows that forḡ 596 sufficiently large, Π has at least two fixed points, one of which is stable.

597
Stable fixed points d of Π are associated with stable n − n solutions of the 598 flow system, and can be interpreted as the maximum values of d at the start 599 of the cycle (see Figure 18B). 600 We have previously shown that for any finite range of coupling strengths 601 withḡ < g sup , there exists an infinite number of maps Π n that have a stable 602 fixed point. This also means that Π has stable fixed points linked to stable 603 n − n solutions for any n sufficiently large. In other words, given that our 604 reduction assumptions hold, a two-cell reciprocally inhibitory network with 605 depression can produce an infinite number of anti-phase burst patterns (of 606 increasing period) on the coupling strength interval (0,ḡ sup ). 607 We also observe that there exist some values ofḡ, such asḡ = 0.76 608 in Figure 17A, where Π has two stable fixed points. Figure 18A Figure 18B  shown the existence of bistability in a very similar network for n < 3. 616 We compute these fixed points numerically to derive the full bifurcation 617 diagram. Figure 19 illustrates two bifurcation diagrams: One diagram that 618 has been used in Figure 5 (blue lines) and computed by simulating the orig-

629
To study the mechanisms of burst generation in tadpole struggling we