Pairing cellular and synaptic dynamics into building blocks of rhythmic neural circuits. A tutorial

In this study we focus on two subnetworks common in the circuitry of swim central pattern generators (CPGs) in the sea slugs, Melibe leonina and Dendronotus iris and show that they are independently capable of stably producing emergent network bursting. This observation raises the question of whether the coordination of redundant bursting mechanisms plays a role in the generation of rhythm and its regulation in the given swim CPGs. To address this question, we investigate two pairwise rhythm-generating networks and examine the properties of their fundamental components: cellular and synaptic, which are crucial for proper network assembly and its stable function. We perform a slow-fast decomposition analysis of cellular dynamics and highlight its significant bifurcations occurring in isolated and coupled neurons. A novel model for slow synapses with high filtering efficiency and temporal delay is also introduced and examined. Our findings demonstrate the existence of two modes of oscillation in bicellular rhythm-generating networks with network hysteresis: i) a half-center oscillator and ii) an excitatory-inhibitory pair. These 2-cell networks offer potential as common building blocks combined in modular organization of larger neural circuits preserving robust network hysteresis.

(5) the dynamics of its inactivation gating variable h is given by , where (6) and ( 7) b h (V ) = 1 1 + e (55 V s )/10 , and (10) The fast potassium I K -current is given by the equation with the reversal potential E K = 75mV and the maximal conductance set as g K = 0.3nS.The dynamics of inactivation gating variable are described by and ( 14) , where a n (V ) = 0.01 55 V s e (55 V s )/10 1 (16) The rapidly depolarizing h-current is given by with E h = 70mV, and g h = 0.0006nS; the dynamics of its y-activation probability is described by whereas the leak current is given by with E L = 40mV, and g L = 0.003nS.The sub-group of slow currents in the model includes the TTXresistant sodium and calcium I T -current given by with E I = 30mV and g T = 0.01nS, where the dynamics of its slow activation variable are described by with the time constant t x set as 100 or 235ms in this study.The slowest outward Ca 2+ activated K + current given by with E K = 75mV, and g KCa = 0.03, where the dynamics of intracellular calcium concentration is governed by with the Nernst reversal potential E Ca = 140mV, and small constants r = 0.0003ms 1 and K c = 0.0085mV 1 .

APPENDIX II: IMPLEMENTATION DETAILS FOR FAST-SLOW DECOMPOSITION
In this appendix, we describe the computational approaches used to determine the relative positions of the slow nullclines, the SNIC bifurcation curve, and the average quantities needed to locate periodic orbits, as well as to create the heat maps used in the figures of the main text.

Approach I: Equations of slow nullclines
To locate the positions of the slow nullclines x ¨= 0 and [Ca] ¨= 0 in the ([Ca] x)-phase plane in the figures above, we first parametrize the membrane voltage V n within an equilibrium range [ 70; +20]mV with some step size.Next, solve the equilibrium state equation V ¨= 0 for the I KCa -current: with the corresponding static functions h ô (V n ), n ô (V n ), and m ô (V n ) for equilibria in the currents above, and next using Eq. ( 24) find the parameterized array To determine the position of the nullcline [Ca] ¨= 0, first solve the equation V ¨= 0 for the I T -current: and then from Eq. ( 25) find and use the ordered pairs

Approach II: Averaging approach to identify periodic orbits
To appraise and quantify the fast spiking behavior of the model via its fast subsystem, we need to determine the location of a corresponding stable orbit on the M po -manifold and its 2D-projection, as depicted in Figs.4A and B,respectively. Following Refs.[62,81,82], one numerically solves the following averaged system adopted from Eqs. ( 10)-( 11) for a pair (Öxã ò , Ö[Ca]ã ò ): where V po (t), x po (t), and [Ca po (t)] are the coordinates of the periodic orbit with period T , situated on the critical 2D M po -manifold at the given (D Ca , D x )-parameter values.Here Öxã is the average quantity from Eq. ( 12).The solution (Öxã ò , Ö[Ca]ã ò ) of the average system above corresponds to "gravity center"of the sought periodic orbit in the phase pace of the swim interneuron model.
In practice, there is a (computationally) cheap trick for finding a loose approximation of the location and the gravity center of such a tonic spiking orbit in the phase space of the swim interneuron model.By setting t x = 1000s or larger values one makes the dynamic x-variable slow so that it cannot echo fast voltage oscillations, but traces out the sought average Öxã nullcline in the ([Ca], x)-phase plane.Its intersection with the calcium nullcline [Ca] ¨= 0 pin-points a location of the sought periodic orbit on the tonic-spiking manifold M po .

Bi-parametric sweep of neuronal dynamics
The second approach employs contour lines identified in bi-parametric sweeps of the fast subsystem, where the x and [Ca] variables become frozen to be treated as control parameters in the fast subsystem of the swim interneuron model.First, we create a grid (of [300] ✓ [300]-resolution) of ordered ([Ca], x)pairs over the slow phase plane.For each pair, a phase trajectory of the model is computed and the initial transient is discarded.The trajectory is classified as either a spiker if the voltage has three or more zero crossings, and quiescent otherwise.Each case is treated separately but in both cases, four quantities are characterized: spike frequency, and average values ÖV ã, Öxã and Ö[Ca]ã.When the system is quiescent, the frequency is 0, the average voltage is the final voltage, which is simply plugged into the x and [Ca] state equations to determine Öxã and Ö[Ca]ã-values.
In the case of tonic-spiking activity, we use the final three zero crossings of the entire trajectory to determine the established spike frequency and the average voltage.The x-and [Ca]-coordinates are calculated for every point in the sampling of the limit cycle and subsequently averaged.This order of operations is required due to nonlinearities in the model dynamics.The nullclines are plotted as the zero contours of the resulting Öxã and Ö[Ca]ã arrays, while the SNIC-curve itself is evaluated as the zero contours of the frequency array.The average voltage is only used to generate the corresponding heat-map color-painting in the phase plane with a red hue for its oscillatory part above the SNIC-curve, and below it with a blue hue representing quiescent phases.

APPENDIX III: SEQUENTIAL PHASE PLANE ANALYSIS OF NETWORK PERIODIC ORBITS
In order to validate the mechanism underlying rhythmogenesis in the two network models we visualize the movement of the nullclines over the course of a network bursting cycle.We partitioned a network bursting trajectory into subsets where all S(t) were monotonic in time and visualized the transformations of the geometry as the network oscillations unfolded.During our investigation of the excitatory-inhibitory module, we noticed that it was also possible to obtain network oscillations with frozen x and [Ca]-variables, suggesting that the mechanism of rhythmogenesis was not a product of the slow phase dynamics alone.Further investigation into the excitatory-inhibitory module is required.
We showed above that in the case of the HCO network, the onset of network oscillations emerging from voltage transients is coordinated by the cyclicly moving x-nullclines x ¨= 0 in the slow ([Ca], x)-phase plane for both neurons.Here, we include a supplementary Fig. 1 to further detail and validate our analysis.Its panels A and D for neurons 1 and 2, respectively, demonstrate representative snapshots of time-varying Figure 1.Illustration of the dynamic evolution of the slow nullclines, equilibrium x ¨= 0, average Öxã, and calcium [Ca] ¨= 0, in the slow [Ca, x]-phase plane overlaid with the thick semicircular bursting orbits for both neurons in the HCO network over its half period.In all panels, the color spectrum, varying from the blue to the red color, is meant to communicate time progression, whereas black and white lines are associated with neurons 1 and 2, resp.The thick lines in Panel B represent the smoothed the voltage traces V (t), whereas the thin one represent the original ones.The time progressions of the synaptic variables S post (t) for both neurons are presented in Panel C. Panels A and D demonstrate small deviations of the calcium nullcline [Ca] ¨= 0 over a burst cycle, while the perturbed equilibrium nullclines x ¨= 0 in the post-synaptic neurons are shifted by the inhibition far leftwards in the phase plane compared with unperturbed ones on the right in the pre-synaptic neurons.
rearrangements of the nullclines, equilibrium x ¨= 0, average Öxã and calcium [Ca] ¨= 0, in a chronological sequence over a half of each burst cycle.Since the HCO network is symmetric, the neurons will next switch their roles to complete the full bursting cycle.Close examination yields an interpretation by breaking a half-cycle sequence into two stages: Stage I.As the active pre-synaptic neuron 2 (its initial/final positions given by blue/red dots in Panel D) drifts rightwards and approaches the SNIC-curve (not shown) along the average Öxã nullcline, its spike frequency decreases so that the synaptic probability S 2 (t) (black line in Panel C) drops to zero.The concurrent decay of inhibition projected from neuron 2 onto neuron 1 (its initial/final positions given by blue/red dots in Panel A) lets the corresponding x 1 -variable rise, which in turn makes the average Öxã nullcline for neuron 2 shift to the right.Stage II.After the x 1 -variable reaches the average Öxã nullcline (color-coded red line in Panel A), neuron 1 starts spiking and hence inhibiting neuron 2. The inhibition causes the average Öxã nullcline for neuron 2 to shift leftwards so that the x 2 variable drop down close to zero and that the corresponding phase trajectory of neuron 2 is about to start traversing along the forced (red line) equilibrium nullcline x ¨= 0 in the phase plane in Fig. 1D, and so forth.
To smooth a V (t) trace of the active spiking neurons, it is first partitioned to ensure that each spike belongs to a unique subset by monitoring a threshold set around -35mV.A maximum spike time of 1s determines whether each subset contains a spike or a quiescent episode of the burst.Spike-containing subsets are replaced with their voltage averages, while quiescent voltage level are left intact.The partitioning and averaging for each S(t) follow their corresponding presynaptic voltages.Finally, a composition of two moving averages with a 0.5s window is applied as a filter for further smothering.