A biophysical model explains the spontaneous bursting behavior in the developing retina

During early development, waves of activity propagate across the retina and play a key role in the proper wiring of the early visual system. During a particular phase of the retina development (stage II) these waves are triggered by a transient network of neurons, called Starburst Amacrine Cells (SACs), showing a bursting activity which disappears upon further maturation. The underlying mechanisms of the spontaneous bursting and the transient excitability of immature SACs are not completely clear yet. While several models have attempted to reproduce retinal waves, none of them is able to mimic the rhythmic autonomous bursting of individual SACs and reveal how these cells change their intrinsic properties during development. Here, we introduce a mathematical model, grounded on biophysics, which enables us to reproduce the bursting activity of SACs and to propose a plausible, generic and robust, mechanism that generates it. The core parameters controlling repetitive firing are fast depolarizing V-gated calcium channels and hyperpolarizing V-gated potassium channels. The quiescent phase of bursting is controlled by a slow after hyperpolarization (sAHP), mediated by calcium-dependent potassium channels. Based on a bifurcation analysis we show how biophysical parameters, regulating calcium and potassium activity, control the spontaneously occurring fast oscillatory activity followed by long refractory periods in individual SACs. We make a testable experimental prediction on the role of voltage-dependent potassium channels on the excitability properties of SACs and on the evolution of this excitability along development. We also propose an explanation on how SACs can exhibit a large variability in their bursting periods, as observed experimentally within a SACs network as well as across different species, yet based on a simple, unique, mechanism. As we discuss, these observations at the cellular level have a deep impact on the retinal waves description.

Retinal waves, observed in many vertebrate species -chicks 1 , ferrets 2 , mice 3 , turtles 4 , macaques 5 etc. are spontaneous bursts of activity propagating in the developing retina and playing a fundamental role in shaping the visual system and retinal circuitry. They emerge due to the conjunction of intrinsic single-cell properties (excitability and long refractoriness) and network interactions 6 . In the developing retina, waves evolve during three consecutive stages, mainly characterized by different types of synaptic transmission; gap junctions (stage I), acetylcholine (stage II) and glutamate (stage III) 7 and are mediated by transient networks of specific cell types. An important type of such cells are the Starburst Amacrine Cells (SACs), which are involved in direction selectivity in the mature retina 8,9 thanks to their special asymmetric synaptic distribution of excitatory and inhibitory synapses. However, during early development, these cells have a different role; they are responsible for eliciting stage II retinal waves via their spontaneous bursting behavior induced by autonomous intrinsic cellular mechanisms 6 . This activity disappears completely upon maturation 6 . During stage II, SACs constitute a transient network of autonomous bursters which, through cholinergic coupling, achieve opportunistic local synchrony. Under conditions yet to be identified, this local area of synchronized bursters can give rise to propagating waves. Although stage II retinal waves are quite thoroughly studied experimentally 1-3 as well as theoretically by several modeling approaches [10][11][12][13] (see Methods for a short review of these models and a comparison to ours) there are still important questions not yet addressed by existing studies. So far, previous models do not allow to understand how biophysical parameters, like conductances, impact the activity of SACs during this stage. This question is important as some of these parameters evolve during development, or can be controlled pharmacologically, with drastic effects on individual SACs bursting 6,14 . Finally, and although stage II retinal waves are ubiquitous, experiments reveal a vast variability in the bursting period of SACs across species. Are SACs fundamentally different from one species to another, or is there a common mechanism of bursting under which the period is quite sensitive to physiological parameters variations? In the spirit of previous works [15][16][17] using dynamical systems and bifurcation theory to characterize neural systems, we develop a thorough modeling of the SACs bursting which accurately reproduces the experimental observations made in the seminal papers by Zheng et al. 6,14 , several of them for the first time. In addition to new theoretical predictions, we also propose a simple scenario for the SACs bursting mechanism and its evolution during development, that could be tested experimentally. Bursting is thoroughly studied in the neuronal modeling literature 18 , as well as the possible bifurcations scenarios generating it. Through our analysis, we identify which parameters are essential to control the transient bursting of immature SACs occurring naturally along development or pharmacologically. As parameters variation is restricted to small ranges by biophysical constraints, our hypothesis is that, during development, SACs are close to a bifurcation point, consequently making them sensitive to small physiological parameters variations.
More precisely, we propose to describe and explain the experimental observations of Zheng et al. 6,14 on individual immature SACs dynamics, in the following theoretical framework: SACs burst because, during stage II, they are close to a (saddle-node) bifurcation, driving the cell from a rest state to fast oscillations. The bifurcation is triggered by small variations of the membrane potential, due to noise or other cellular effect. After this bifurcation, SACs stay in the regime of fast repetitive firing, triggering the activation of slow calcium-gated potassium channels generating a slow hyperpolarization current (sAHP). This current eventually drives the cell to a homoclinic bifurcation. As a consequence, cells return back to the quiescent state allowing a long after hyperpolarization phase. This bifurcation analysis reveals a mechanism for SACs rhythmic bursting, interpretable with respect to biophysics and the experiments of Zheng et al. 6,14 on individual immature SACs dynamics. In addition, we identify the fast voltage-gated potassium channels as a possible candidate for the fast hyperpolarizing channels involved in bursting (to our best knowledge not yet been experimentally determined). Based on 19,20 we propose they could be of Kv 3 type. This leads us to a conjecture on the role of the Kv 3 potassium conductance variation which may act on the transient intrinsic properties of SACs excitability during development. Furthermore, the closeness to a bifurcation has strong consequences on the interburst intervals of SACs. It is naturally associated with a high variability which can explain the wide range of interburst periods observed within a network 21 , as well as across species 10 . As we show, a small external current (~5 pA) strongly impacts the interburst intervals of SACs. This implies, that, in a network of cholinergically coupled SACs, the bursting period has a wide intrinsic variability, resulting from simple bifurcation arguments, without requiring exogenous mechanisms. In the discussion section, we briefly argue how a set of such inhomogeneous bursters may achieve synchrony and wave generations without any additional mechanism than the existing cholinergic interaction. The general discussion on how cholinergic interaction between nearby autonomous bursters generates wave-like, population-level behavior propagating in a similar way as is observed experimentally will be addressed in a forthcoming paper (see [22][23][24].
The paper is organized as follows. The section 'Results' introduces the model definition, the bifurcation analysis, and the main conclusions resulting from this analysis, including experimental predictions. The section 'Discussion' is mainly dedicated to the biological interpretation of our theoretical results, as well as their impact on the description of retinal waves. Finally, in Methods, we provide the technical aspects of our work including the biophysical derivation of the model's equations, as well as mathematical computations used in the paper. We also include in this section a short review on the comparison of existing models.

Results
During retinal development, Starburst Amacrine Cells (SACs), which are cells forming a lattice right above the ganglion cells layer, are found to exhibit spontaneous intrinsic rhythmic bursting activity, which disappears upon maturation 6 . In the following, we derive a biophysical set of equations, describing the bursting activity of immature SACs, following experimental observations regarding the ionic currents involved in their bursting 6 . We calibrate the equations parameters in order to reproduce the bursting characteristics found in 6,14 . The analysis of our model leads to our central hypothesis, that the rest state of SACs lies near a bifurcation point, which means that small variations of parameters could lead to abrupt changes in the dynamics of the cell. This hypothesis, which is novel for immature SACs, jointly explains A) the biophysical mechanism of bursting in the developing SACs, B) accounts for the vast variability of bursting characteristics of SACs across species and C) predicts a possible mechanism on how SACs lose their excitability upon maturation.
Let us now start by deriving our model along with the biological justifications for our equations.
A biophysical model for bursting immature sACs. In order to address the mechanism of bursting and how the physiological parameters of SACs change during development we need to model their underlying cellular mechanisms. In the case of immature SACs, the two key biophysical mechanisms associated with the emergence of spontaneous bursting during early development are 6 : • fast repetitive bursts of spikes mainly controlled by fast voltage-gated channels; • prolonged After HyperPolarizations (AHPs) modulating fast oscillations, controlled by Ca 2+ -gated K + channels.
In order to account for (i) and (ii), we model immature SACs activity with a conductance based model of the Morris-Lecar type 25 , an electronic-circuit-equivalent equation for the local membrane potential, with additional non-linear currents. Our model involves 5 variables, whose evolution is controlled by a set of non-linear Scientific RepoRts | (2019) 9:1859 | https://doi.org/10.1038/s41598-018-38299-4 differential equations with 3 timescales: (1) fast variables (~10 ms) V(t), the local membrane potential, and, N(t), the gating variable for fast voltage-gated K + channels; (2) medium variable (~2 s) C(t), the intracellular Ca 2+ concentration; (3) slow variables (~10 s) R(t) and S(t), the gating variables for slow Ca 2+ -gated K + channels. All parameters values and the auxiliary functions involved are found in Methods. The membrane voltage V(t) obeys: where C m is the membrane capacitance, I L = −g L (V − V L ) is the leak current, with g L , leak conductance and V L , leak reversal potential. The fast repetitive firing (see (i)) during the active phase of bursting generally results from the competition between a depolarizing and a hyperpolarizing current. Experiments in 6 have specified voltage-gated Ca 2+ channels as the depolarizing component of bursting in SACs. However, the ionic channels related to the fast hyperpolarizing component of SACs bursting have not yet been identified experimentally. In this work, we propose that fast voltage-gated K + channels play this role (see also Discussion for further justification). Therefore, the terms I C and I K , respectively corresponding to Ca 2+ and K + currents, are generating the fast oscillations in our model.   (19) in Methods) when the constant current I tot is varied in a wide range. Each deep blue curve corresponds to the projection of a trajectory of the fast subsystem onto the V axis, for a fixed value of I tot . (B) Interpretation of all the different dynamical regimes of the bifurcation diagram. (C) Same as in A, except that the various trajectories are now drawn in the (N, V) plane and not only projected onto the V axis. Blue. Evolution of the limit cycles when I tot varies. Red. Homoclinic orbit. Bursting stops by an homoclinisation of the limit cycle. Note that 6 have shown that voltage-gated Na + channels do not participate in the bursting mechanism of immature SACs (bursting activity of SACs was not altered upon tetrodoxin -TTX-application), thus, dynamics of Na + channels will not be considered in our modeling (see also Discussion).
The voltage-gated Ca 2+ current is described by a Morris-Lecar model 25 where: is the voltage dependent conductance of the Ca 2+ channel (see Eq. (21) in Methods). The fast voltage-gated K + channel is modeled as: where the evolution of the fast voltage-gated K + channel gating variable N(t) is given by: N Λ(V) and ∞ N V ( ) are given by Eqs (21) and (22) in Methods. The form of equations (3) and (4) is based on 25 . However, the timescale of the variable N is much faster (~ms) in our model that in the original Morris-Lecar model (~s), calibrated in order to capture the frequency of the fast repetitive firing of SACs, which is around 20 Hz 6 (see also Bifurcation Analysis).
The long refractoriness in-between consecutive bursts is controlled by a slow After Hyper-Polarization (sAHP) K + current, I sAHP . It was observed by 6 that I sAHP is mediated by Ca 2+ -gated K + channels, and that it resembles the sAHP observed by Abel et al. 26 , generated by specific channels called SK. Following these tracks we propose a modeling of "SK"-like channel (as named in 6 ) based on 12 for the structure of the equations and 27 for the calcium . Red lines correspond to stable fixed points, gray and pink lines to unstable fixed points, green lines are the extremal values of voltage fast oscillations. Upon this bifurcation diagram, obtained from the reduced model (Morris-Lecar with I sAHP = I tot = constant), we superimpose the trajectory of the full system in the plane (I sAHP , V) (deep blue). The blue arrow indicates the direction of the flow. In the fast oscillations regime V is varying periodically, with a fast period, inducing a fast variation of the term V − V K in the sAHP current, explaining the diagonal motion. Despite these fast oscillations one sees that the bifurcation is driven by the slow motion on the branches obtained from the bifurcation analysis assuming a constant current. The conductance g sAHP R 4 is varying slowly, with the time scale of R, explaining the slow left-wise shift of the trajectory until the homoclinic bifurcation point is reached. The cyan vertical line corresponds to the value of I ext = 0. (b) Noise driven bursting. The dynamics has now a stable fixed point (purple circle). An additional current, I ext = −4 pA (vertical cyan line) is now present. Bursting cannot occur without an exogenous excitation, here noise. The amount of noise (controlled by σ) necessary to induce bursting depends on the distance between the fixed point (controlled by I ext ) and the bifurcation point I SN 1 . dynamics. The mechanism of the opening of Ca 2+ -gated ionic channels channels is analyzed in detail in the "Methods" section. In order to simplify the cascade of chemical reactions taking place while opening the sAHP channels, we approximate the channel dynamics by reducing the process into two discrete steps: a) Four ions of Ca 2+ bind to a second messenger protein called calmodulin, forming a saturated calmodulin complex, CaM; b) CaM binds to each of the four intracellular subunits of the channel to open it (see Fig. 1). This process is mapped to our model through three variables: 1) the variable C which models the intracellular calcium concentration and mainly controls the gating variables of the sAHP channels, 2) the variable S which models the fraction of the saturated calmodulin and 3) the variable R which models the fraction of bounded terminals. This gating mechanism is sketched in Fig. 1. The sAHP current takes the form: where g sAHP is the maximum sAHP conductance. Indeed, 4 bound terminals are needed to open a Ca 2+ -gated K + channel, thus the corresponding conductance is g sAHP R 4 , involving a fourth order nonlinearity.
We model the gating mechanism of the Ca 2+ gated K + channels as follows. The gating variable R(t) obeys: R R the fraction of saturated calmodulin concentration S(t) obeys: S S 4 and the intracellular calcium concentration C(t) obeys: The derivation of the equations (6), (7) and (8) is fully justified in Methods. Finally, ξ t is a white noise whose amplitude is constant with time and controlled by σ and I ext is an external current. For subsequent analysis we introduce the current I tot = I sAHP + I ext . (1) Fast spiking occurs due to the competition between voltage gated Ca 2+ (excitatory current) and K + channels (inhibitory current). (2) Calcium load increases during the rapid firing phase while the voltage is in a high voltage state. This leads to a slow increase of sAHP. When sAHP is large enough there is a sharp decay of the voltage (bifurcation) and the hyperpolarisation phase of the cell starts. As voltage is low, the calcium load starts to slowly decrease. (3) During the slow offloading calcium stage we observe a refractory phase. The decrease of calcium concentration induces a slow decay of the sAHP. When calcium is small enough rapid firing starts again via a new bifurcation. Dynamical changes with respect to parameters variations. In our model, the variable V, N evolve with a fast time scale of the order of a few milliseconds under the influence of the current I sAHP whose conductance is slow, being driven by slow variables C, S, R (time scale -several seconds). As a consequence, the main material of this paper relies on a thorough analysis of the bifurcation structure in the fast Morris-Lecar dynamics. There already exist remarkable numerical bifurcation analyses of the Morris-Lecar model (see for example the very detailed work of 28 ) which has however many parameters: changing their range of variations can dramatically impact the dynamics. Our study on SACs bursting considers quite different range of parameters than the studies we know about (more related to the original paper investigating the barnacle giant muscle 25 ). For example 28 , investigates the effect of a positive external current whereas we consider a negative external current, with quite different effects on    We distinguished two classes of parameters. The first one consists of modeling parameters constraining the Morris-Lecar dynamics, namely V 1 , V 2 , V 3 , V 4 (see Methods for the corresponding equations). We verified the robustness of the scenario described in the paper when varying these parameters. Especially, V 3 is the half activation potential of the fast potassium channels, playing a central role in SACs excitability as shown later in the text. The second class is the set of parameters constrained by biophysics such as conductances, reversal potentials, capacitance,... Most of them have been fixed based on the biophysical literature and are not considered to vary. Here, we choose to study the variation of 4 of these parameters: I tot , because it controls bursting, and g K , g C and V 3 because, as we show later, their variations allows us to reproduce experimental facts, that have not been reproduced before by any model. The most striking consequence is a conjecture on the role of the potassium conductance during development.

Bifurcation analysis.
Bursting is an alternation between a rest state and repetitive firing, often modulated by slow voltage or Ca 2+ dependent processes 18 . In order to address which are the biophysical mechanisms for a) triggering, and, b) sustaining bursting in immature SACs, we will perform bifurcations analysis. This is a powerful mathematical tool providing a compact, geometric representation, of dynamical systems behavior when Figure 9. Bifurcation diagram as a function of I tot , g K . At the bottom center, the bifurcation map is plotted. Points are sampling points obtained by a direct inspection of the phase portrait whereas continuous lines correspond to fit. "Hc" means "homoclinic bifurcation" and "SN" stands for "saddle-node". The capital letters (A-E) inside correspond to the typical phase portraits surrounding the bifurcation map. In these phase portraits, continuous lines correspond to nullclines (N V , in black, is the V-nullcline whereas N N , in gray, is the N-nullcline). 'Si' stands for 'Sink' , 'Sd' for 'Saddle' , 'SF' for 'Stable Focus' and 'UF' for 'Unstable Focus' . Sink and stable focus correspond to stable rest state: a small perturbation about this state decay exponentially fast. On the opposite, saddle and unstable focus are unstable. For saddles we show the stable and unstable directions (black lines) as well as the stable (cyan) and unstable (red) manifolds W s , W u . 'SPo' means 'Stable Periodic Orbit' . It corresponds to fast oscillations, plotted in dark blue. In region (A), there is a low voltage stable state. In region (B), a stable state with low voltage coexists with a limit cycle (fast oscillations) separated by an unstable state. When the cell is in the low voltage state, a large enough perturbation (e.g. noise or other cells action) leads it to fast oscillations. In region (C), two stable states, one with low voltage and one with high voltage, coexist separated by an unstable point. In region (D), the cell only exhibits a fast oscillations regime, reached whatever the initial condition. Finally, in region (E), there is a stable state with high voltage. When varying parameters, continuous change occurs in the interior of a region (e.g. the period of oscillations continuously varies), whereas crossing the bifurcation lines leads to abrupt changes. parameters are modified. While, in general, a continuous variation of a parameter induces a continuous variation of the dynamical system's solutions, the variation of parameters around bifurcation points induces, in contrast, strong and abrupt changes in the dynamics. We first perform a bifurcation analysis of the fast V, N dynamics in the presence of the current I tot , constant in time, and used as a bifurcation parameter. The bifurcation diagram is represented in Fig. 2A (drawn using MATCONT 29 ), where we have explored a wide range of variation of I tot (−70 to +310 pA). For the analysis, we consider an example with the potassium g K = 10 nS, and calcium conductance g C = 12 nS (see Robustness with respect to parameters variations). We reveal a rich dynamical profile upon the variation of an external current where for I tot < −3.7 pA, stable rest state co-exists with unstable fixed points, for −3.7 < I tot < 250 pA, cell oscillates constantly due to a limit cycle and for I tot > 250 pA only a stable rest state at a high voltage exists (see Fig. 2B). We refer the reader to e.g. 30 for the dynamical systems terminology used in the next paragraphs.
In more detail, when − .
rest mV, coexisting with two unstable fixed points, a saddle and an unstable focus (see Fig. 9A for the phase portrait of this case).
tot H c , there is an homoclinic bifurcation giving rise to a stable limit cycle. In the range pA, the stable rest state coexists then with a stable limit cycle, corresponding to fast oscillations, and an unstable fixed point.
pA the stable rest state coalesces with the middle unstable branch and both disappear with a saddle-node bifurcation (SN 1 ). For > I I tot S N 1 the dynamics has only one attractor, the limit cycle, corresponding to fast oscillations. This cycle eventually disappears by a Hopf bifurcation at I tot = 250 pA (Fig. 2). However, this value of external current is quite beyond the range of plausible values during the bursting activity of SACs.
We observe that the rest state of SACs, corresponding to I tot = 0 pA, lies near a saddle-node bifurcation point within a narrow region of a few pA, which is our main emerging hypothesis. This leads to important consequences for SACs described in the following.
Bifurcations analysis reveals two possible drives for triggering bursting. The bifurcations analysis for the fast subsystem shown in Fig. 2, reveals all possible dynamical regimes upon the variation of a constant external current within a wide interval of our choice. Thanks to the timescales separation in our model (fast, medium, slow), we make the approximation that the bifurcation diagram in Fig. 2A, does not change significantly when we consider the full model case. However, in the full model I tot is not constant, but time-dependent, driven by sAHP. By superimposing the solution of the full model and the bifurcation diagram of the fast subsystem, we reveal how the trajectory driven by the slow variations of I sAHP matches the bifurcation diagram, in the plane (I tot , V). See Fig. 3 showing a zoom of the bifurcation diagram shown in Fig. 2A, to our regime of interest (stable branch: red, neutral saddles: gray, unstable branch: pink, green traces showing the maximum and minimum of V along the limit cycle), with the superimposed trajectories of the full system (blue curves). Analyzing how the slow variations of I sAHP generate bursting, we actually exhibit two possible drives of triggering bursting: (i) Dynamical drive, induced by the mere dynamics of neurons, without external influence. In our model, the joint fast dynamics of Ca 2+ and K + channels generates fast oscillations while the slow AHP, mediated by Ca 2+ gated K + channels, modulates slow oscillations 6 . The conjunction of these two mechanisms generates bursting. Note that in this case, I ext = 0 in Eq. (1), meaning there is no external current (Fig. 3 left). In this scenario, cells burst periodically, with a frequency controlled by the characteristic times τ R , τ S of variables R and S respectively. From a dynamical systems perspective, there are two important bifurcations associated with bursting: a) bifurcation of the rest state that leads repetitive firing and b) bifurcation of a spiking attractor that mV whereas the heat map shows calcium peaks lasting more than 1 s in this region. This is due to the particular structure of the phase portrait of the Morris-Lecar model in this region and the effect of noise, as shown in (c). This phase portrait has two saddles (Sd) down (Sd d ) and up (Sd u ), and two sinks (Si), down (Si d ) and up (Si u ). The unstable manifold (W u , in red) of Sd d connects to Si u , on the right, and to Si d , on the left; The unstable manifold of Sd u connects to Si u , on the right, and to Si d , on the left. Thanks to noise, dynamics can leave Si d and reach Si u where it stays until sAHP rising drive it back to Si u . This is illustrated by the blue trajectory (complete system with noise) which follows the unstable manifolds of the Morris Lecar model. The inset shows the trajectory of the voltage for the complete system in the presence of noise. (ii) (d) While the bifurcation diagram predicts bursting for V 3 ~ −28 mV, g K ~ 12 nS, the heat map does not show bursting in this region (black). This is again due to noise. When the cell is on low rest state (Si) noise drives it to the limit cycle, but after a few oscillations, noise drive the cell back to the rest state Sd (see inset for the voltage trajectory). This effect depends on the distance between the two fixed points Si, Sd, fixed by the model parameters, and on the noise intensity σ (here 4 pA ms 1 2  leads to a rest state 18 . There are two possible types of bifurcations linked to such type of behavior (1) Homoclinic saddle-node bifurcation and (2) Hopf bifurcation. We have therefore shown that (1) holds in our case. The first bifurcation associated with the breaking of the rest state is the saddle-node bifurcation. The "spiking attractor", here the limit cycle, is created by an homoclinisation. According to the classification of bursting made by Izhikevich in 18 , this corresponds to a "square-wave" point-cycle planar burster. (ii) Noise-drive, induced by an external excitation, typically noise. In order to move to this regime, it suffices to add an external constant current < = − .
to Eq. (1). Then, the deterministic dynamics converges towards a stable fixed point, the cell ends up in a steady rest state and cannot switch to rapid firing. In this case, the bifurcation diagram (Fig. 3 right) has now two fixed points: V rest , which is stable (purple point), and V u which is unstable (intersection of the gray branch and the vertical blue line). However, in the presence of noise, sufficiently close to the bifurcation point (depending on the noise amplitude), the random fluctuations around the rest state leads eventually the cell into the rapid firing regime (Fig. 3 right). Note that when I ext is sufficiently large in absolute value the total current I tot is smaller than I SN 1 whatever the (negative) value of I sAHP .
Taken together, since noise is always present in neurons, distinguishing those two scenarios for triggering bursting is not easy experimentally but still possible (see Discussion). Maybe both scenarios are possible in the same network of SACs, since it would correspond to a slight variability in the intrinsic parameters which control the rest state of each SAC, placing it on the right or left of the bifurcation within a narrow region (~10 pA according to our model). We are actually able to capture both the dynamically and the noise driven bursting regimes in a single, unifying model. In the presence of noise the dynamically driven bursting regime still exists but the period has now random fluctuations. In fact, the noise smooths the transition between those 2 regimes. Also, what our analysis reveals is that, in the noise driven regime, bursting will occur only if the noise has a sufficiently large amplitude compared to the distance between the rest state and the bifurcation point. On the other hand, in the dynamically driven regime, bursting will occur only if the rest state is sufficiently close to the bifurcation in order for the amplitude of the negative sAHP current to be sufficient to compensate the depolarizing current (I C ) and enter the hyperpolarization phase. Otherwise, if the rest state is further on right from the bifurcation, SACs would exhibit only a repetitive firing regime without hyperpolarization. In other words, in both regimes, the cell has to be close enough to the bifurcation point to observe bursting. Bifurcations analysis reveals a biophysical mechanism for SACs sustained bursting. Our analysis reveals also a possible biophysical mechanism to explain how bursting activity is sustained periodically in immature SACs. Once fast repetitive firing is elicited in SACs, either through i) dynamical, or ii) noise drive, the biophysical cyclic mechanism which produces bursting is the same. From the bifurcation analysis we see that once the cell crosses the saddle-node SN 1 , it is firing rapidly (limit cycle). Being in a high voltage state, Ca 2+ loads, leading to a slow rise of sAHP conductance. This corresponds, in the bifurcation diagram, to a motion towards negative current values, with a slow variation due to the sAHP conductance driven by the slow variable R, and fast oscillations due the fast variable V − V K appearing in the sAHP current I sAHP = −g sAHP R 4 (V − V K ). As revealed by a bifurcation analysis where the sAHP conductance is the bifurcation parameter (not shown), and as shown in Fig. 3, these fast oscillations have no effect on the bifurcation, what matters is the slow variation. The increase of sAHP conductance continues until the cell reaches the homoclinic bifurcation where rapid firing stop. I tot = I sAHP is very negative now, leading to the hyperpolarization phase. Then, because voltage is low, [Ca 2+ ] decreases and I sAHP drops down. Eventually, SAC crosses the SN 1 bifurcation point and starts firing again either due to dynamics in one case or noise. This is illustrated in Fig. 3. The deep blue traces are the trajectory of a burst in the plane (I tot = I sAHP , V), for the full system. One sees in particular the fast oscillations coming from the fast variable V − V K . Nevertheless, for both proposed scenarios the biophysical mechanism for sustaining bursting is the same. Although the conclusion drawn on the model is entirely based on a mathematical analysis via bifurcation theory, we have sketched the main mechanism of SACs bursting in Fig. 4, for readers not familiar with bifurcation theory. Note that Fig. 4 is only qualitative. In contrast, bifurcation analysis gives us more insights on bursting mechanism, where actually two regimes are possible, and quantitative results, especially on the parameters determining the transitions inducing bursting and the range of these parameters. To our best knowledge this is the first model able to provide a mechanism to account for the repetitive firing of SACs, which is the condition under which calcium controls the after-hyperpolarization phase. Figure 13. Bifurcation diagram as a function of g C , g K in the noise driven bursting regime. Same representation as Fig. 9. In region (A), there is a unique stable rest states with low voltage. In region (B) two stable rest state with low voltage coexist separated by an unstable point; hence dynamics is bistable. In region (C), a stable rest state with low voltage coexists with a limit cycle (fast oscillations) separated by an unstable state. When the cell is in the low voltage rest state, a large enough perturbation leads it to fast oscillations. In region (D) the cell only exhibits a fast oscillations regime, reached whatever the initial condition. In region (E) two stable rest states, one with low voltage and one with high voltage, coexist separated by an unstable point. Region (F) has only one high voltage stable rest state. Finally, in region (G), there is a low voltage rest stable state. explaining the wide range of interburst intervals (IBI) across species by a unique mechanism.
Spontaneous bursting activity in immature SACs has been consistently observed across various species i.e. mice, rabbits, chicks, turtles, macaques etc., although, the measured mean interburst interval (τ IBI ) strongly varies [1][2][3][4] . This experimental observation raises the question whether there exist several bursting mechanisms changing from species, explaining these variations in interburst. Bifurcations theory provides another explanation though. As we showed, in our model, bursting involves a saddle-node bifurcation. It is also known from 31 , that interburst is known to be sensitive to parameters variations near this type of bifurcation. In Fig. 5, we illustrate this effect by numerically computing the heat map of the mean interburst intervals of SACs, while varying the parameters g K and g C , for both bursting scenarios (noise and dynamic driven). We show that the IBI increases when approaching the Saddle-Node line. Note that we define a burst numerically by a thresholding on Ca 2+ : a burst occurs when this concentration exceeds 150 nM during a time larger than 1 s. There is a "critical" region (yellow zone) where interburst increases (up to several minutes). This is nicely explained by a bifurcation analysis (see section Robustness with respect to parameters variations). The irregular shape of this region -alternation between yellow and black regions -is due to our numerical procedure. We generate 20 trajectories of duration 2000 s and count the number of bursts in each trajectory. In our simulation, when going close to the critical region we observe that this sampling is not sufficient; we have rare bursts giving poor statistics, and, in some cases (black regions in-between yellow ones), we observe no burst in the sample.
To illustrate the sensitivity of τ IBI to a parameter variation, we now vary the external current I ext (Fig. 6). This indeed induces strong variations in τ IBI , compatible with the variations observed across species, which increases monotonously and sharply as I ext decreases, following a hyperbola whose form has been derived analytically (see Methods and Fig. 6, blue trace). As shown in the paper 31 , hyperbola is one of the 2 possibilities, corresponding to the case where the noise amplitude is small compared to the excitability. We also observe a sharp transition from a bursting to a non-bursting regime, where τ IBI = 0, which corresponds to the loss of SACs excitability.
Based on the shape of our theoretical curve, we show that τ IBI exhibits a strong asymptotic behavior around a very narrow regime of This result suggests that, depending on slight variations of physiological conditions, a SAC could exhibit very variable bursting periods, providing an explanation how different species exhibit variability in the interburst interval of immature SACs. In real SACs, the inhibitory current I ext could correspond to a change in the intrinsic inhibitory ionic currents of individual SAC or to external inhibitory inputs from other layers of the retina. It could also mimic changes in synaptic inhibition (e.g. GABA A maturation, changes in [Cl] in ). These factors could drive the disappearance of bursting activity and subsequently the loss of SACs excitability. On the other hand, an increase of I ext leads to a gradual decrease of τ IBI towards zero, meaning that in this regime of parameters, SACs tend to burst repeatedly, without a refractory period. This scenario is not observed under normal physiological conditions, but could be tested experimentally by varying I ext as an externally applied current. Finally, in the experimental paper 32 , it is found that the bursting period increases upon maturation (experiment in P1-P2 rabbit SAC). In light of the present model-driven analysis, this could be linked to a change of intrinsic properties of SACs during development.
To compare with experimental results, we show in the same figure τ IBI for different species as found in the literature (Fig. 6). We have not been able to find τ IBI for turtles and chicks. We found the mean inter-wave intervals τ IWI instead. In order to extrapolate to τ IBI in these cases, we used a common constant scaling factor of 3, based on the ratio between τ IBI and τ IWI for rabbits and mice found in 6 and 21 respectively. Besides this narrow regime a decrease of I ext results in a sharp disappearance of bursting activity, which means that bursting activity stops after a limit value.
To conclude, our study suggests that bursting activity in immature SACs could share a common mechanism across species. The same mechanism could also induce interburst variability within a network SACs, thus with cells bursting with a different period. Such variability has been implemented in 21 by varying the decay time of the sAHP from a Gaussian distribution. In our case, the variability arises without adding an extra mechanism and is induced by the hyperbolic behavior observed near a saddle-node bifurcation upon a tiny variation of a single parameter, I ext , interpreted as an effective depolarizing/hyperpolarizing current. Taken together, the vast variability of SACs bursting features accross different species, could be explained simply by our hypothesis that the rest state of SACs lie within a narrow regime around a saddle-node bifurcation point.
Characterizing the effect of noise on the probability distribution of bursts. In 6 , recordings show that the bursting periods of SACs are not regular and have a certain probability distribution. We obtain this effect by adding the noise ξ to the dynamics of the voltage V (see Eq. (1)). Note that noise is always present in real neurons due e.g. to fluctuations in ionic currents from the random opening of ionic channels. What we want to illustrate here is the role played by the intensity of this noise, controlled by the parameter σ. For small σ, dynamics has fluctuations around the deterministic trajectory with little effects during the bursting phase. In contrast, low additive noise during the slow (after-hyperpolarization) dynamics, is enough to accelerate the start of a burst. For higher values of σ the bursting period of cells decreases dramatically. As shown in Fig. 7 the presence of noise has also a drastic impact on the shape of the interbursts intervals distribution. Experimentally, maybe due to the lack of sufficiently large samples in 6 , the exact shape of the distribution of the interburst intervals is not sharply defined. Therefore, the comparison with our theoretical results is difficult. However, the experimental distribution obtained by these authors is definitely not exponential. As we observe an exponential distribution for large enough σ values, Fig. 7, this remark provides us with an upper bound on the level of noise which should be no greater than 8 pA ms 1/2 . Note that the left-wise peak in the bimodal distribution (σ ~ 8 pA ms 1/2 ), corresponds to short bursts. Similar bursts have been observed in the spontaneous activity of retinal ganglion cells during waves, in turtle and in mouse, see e.g. control in Figs 3 and 7 in 4 (E. Sernagor private communication). the role of the fast potassium conductance in bursting activity. We now address in further detail the potential role of the fast voltage-gated potassium channels used in our model to produce bursting (see also the work of E. Marder and collaborators, in a different context 33 ). Zheng et al. 6 propose that the ionic channels mainly involved in the bursting activity of SACs during early development are voltage-gated Ca 2+ channels. In this work the hyperpolarizing current involved is not characterized and these authors don't mention potassium channels before the end where they perform experiments with TEA (see below for more details). We proposed above fast voltage-gated K + channels as a source of fast inhibition necessary for the active phase of bursting of SACs. We now justify this claim based on several experiments made by Zheng et al., interpreted in the context of our bifurcation analysis.
We first show (Fig. 8) how our model accurately reproduces a key experiment of 6 . Here, the authors artificially control the triggering of fast oscillations by applying a short current pulse (I ext = 150 pA for 60 ms) to individual immature SACs. Also, upon the pharmacological application of Cd 2+ , which blocks all Ca 2+ related channels (voltage-gated Ca 2+ and sAHP), they show that no oscillatory activity is triggered upon stimulation, but only a raise in the plateau of the level of the voltage. The corresponding figure of their experiment in the paper 6 has been reproduced in Fig. 8A (with the kind authorization of the authors). As shown in Fig. 8B, we are able to reproduce reliably this result using the fast K + channels described in Eq. (3). Particularly, we are able to simulate the emergence of fast oscillations during a short current pulse, with an AHP phase after the end of the pulse (see gray curve in Fig. 8B). We emulate as well the disappearance of the oscillations observed by these authors, upon blocking all Ca 2+ related channels (voltage-gated Ca 2+ and sAHP), by setting, in our model, the corresponding conductances to zero (see black curve in Fig. 8B). This is interpreted as follows from the bifurcation diagram Fig. 9, in the plane I tot , g K . The simulation has been done in the noise driven bursting regime where the cell, initially in its rest state, is in region A. A pulse of current I ext = 150 pA drives the cell in region D where it fast spikes. Removing the current pulse drives back the cell in region A. Here we observe a small sAHP hyperpolarizing the cell. In order to observe the same behavior in the dynamically driven bursting regime, one needs to clamp the voltage so as to maintain the cell in the rest state.
These observations support therefore our proposition that the hyperpolarizing component of the fast oscillations observed in immature SACs are driven by fast voltage-gated K + channels.

exploring the role of the potassium conductance in the loss of sACs excitability upon maturation.
SACs in the retina of vertebrates lose their ability to spontaneously burst once they reach a certain stage of development different for each species. This transient excitability is a key process in the developing retina and the shaping of the visual system. It is not yet clear which physiological properties of SACs change upon maturation causing the abrupt change from autonomous bursting to rest state. So far, to our best knowledge, no experiment has studied in detail the biophysical properties of immature SACs at the level of each ionic channels involved during this window of development. However, there are some indirect experimental implications which help us extrapolate a possible scenario on how SACs change their properties upon maturation. In 20 Ozaita et al. show that mature SACs not only stop bursting but also cannot be depolarized beyond −20 mV. A specific type of inhibitory voltage gated K + channels, Kv 3 , is responsible for this property, providing an electric shunt to SACs somas. In the contrary, in the developing retina, SACs are autonomous bursters, depolarised beyond mature SACs (~−5 mV) 6 . Therefore, upon maturation, the characteristics of voltage gated K + channels evolve, leading to a drastic change in SACs activity. In our model, these observations can be reproduced by the variation of two parameters; (i) the conductance g K and (ii) the half-activation potential V 3 of the voltage gated K + channels.
Here a point of caution is necessary. The modeling chosen here, in the same spirit as Morris-Lecar or Hodgkin-Huxley model, characterizes the behavior of a small piece of a cellular membrane, large enough to contain many ionic channels so that the description in terms of activity variables (N, R, S) is correct. As a consequence, parameters such as g K , V 3 are phenomenological, depending on membrane characteristics (e.g. the density of ionic channels of a given type), and on channels characteristics. When we vary a parameter such as V 3 , it means that we mathematically change the threshold of the term mimicking fast-K channels in the model. The biophysical interpretation of this variation is more difficult. The easiest interpretation is that it corresponds to a change of threshold response of the corresponding channels, but other explanations could be also possible. Assume for example that there are two types of fast K + channels (e.g. kV3.1 and kV3.2) on the membrane, with different densities and reactivity. Then, a mere change of their density during development would impact the value of the phenomenological parameters g K , V 3 of our model.
With our model, we are able to reproduce an experiment by Zheng et al. 6 . The results of their experiment are shown in Fig. 10, left column (with the kind authorization of the authors), whereas the model results are presented in Fig. 10, right column. Our results are based on a bifurcation analysis in the plane V 3 − g K , which reveals a wide region of parameters where bursting takes place. In Fig. 11 we show, on the top left, the bifurcation diagram in the plane V 3 , g K and on the top right, a heat map, similar to Fig. 5 Fig. 10), this spontaneous activity disappears at P6 (second row, left). We can easily reproduce this observation (Fig. 10, first and second row, right) by moving from region C (fast oscillations) in the bifurcation diagram of Fig. 11 to another region. For example, a simple transition to region A is obtained by decreasing the half-activation potential V 3 of the fast K + channels, leading to the saturation of SACs depolarization below −20 mV observed in mature SACs. A concomitant variation of g K is also possible leading to stop bursting as well. To illustrate this scenario, we set g K = 10 nS and V 3 = −34 mV as an indicative example (see Fig. 10). • How can bursting be restored in further mature SACs? Zheng et al. 6 show that bursting activity can be re-initiated pharmacologically in P8 (stage III) and P22 retinas with chemical agents upon i) blocking all synaptic connections (gap junctions, cholinergic, gabaergic and glutamatergic synapses), especially inhibition by applying an antagonist cocktail and ii) decreasing the conductance of the fast inhibitory K + channels, by applying tetraethylammonium (TEA). This suggests that upon maturation, bursting is in fact suppressed by alterations in the intrinsic properties of individual SACs. Particularly, in the same preparation (P8 isolated SACs), the authors apply a voltage-gated K + channel blocker tetraethylammonium (TEA), which results in restoring bursting activity in isolated mature SACs (Fig. 10, third row, left). A potential type of potassium channels involved is the TEA -sensitive K + of the Kv 3 family. Therefore, this experiment emphasizes once again the potential role of fast potassium channels and especially the Kv 3 family. We mimic the blocking of these channels in our model by decreasing the conductance g K . When g K is small enough, we reproduce the bursting restoration (Fig. 10, third row, right). This example corresponds, in Fig. 11, to a motion from (g K , V 3 ) = (10 nS, −34 mV) (region B) to (g K , V 3 ) = (8 nS, −34 mV) (region C). Note that the chosen value V 3 = −34 mV allows us to fix the maximal depolarization to −20 mV, the value observed by Ozaita et al. 20 for further mature SACs. Remark also that the lower branch of the heat map Fig. 11 (region between G and B in the bifurcation map) corresponds to interburst intervals of the order of 8-16 min, so according to SACs recording not to a biophysically plausible bursting regime. • How can bursting be restored in adult SACs? Upon further maturation, Zheng et al. show that the restoration of bursting in P22 SACs (late stage III-before eye-opening) depends on the interplay between the change of the intrinsic properties of K + channels strong inhibition (TEA application) and blocking the strong inhibitory (gabaergic) input induced by other amacrine cells at this phase (Fig. 10, fourth and fifth row, left). Neither of these conditions suffice alone for the bursting restoration. To model the gabaergic inhibitory input to coupled SACs at P22, we add an inhibitory current I ext = −10 pA whereas the effect of TEA is modeled by decreasing the g K conductance to 4.5 nS. This conductance is not set to 0 because TEA blocks only one subtype of Kv 3 channels while there exist several of them (see below) in the developing retina. The whole operation has the effect of suppressing bursting (red trace top) and corresponds moving from region B to region A in Fig. 9.
Bursting re-initiation is obtained by setting back the current I ext to zero.
Along these lines, we have been able to reproduce the effects of the pharmacological manipulation performed by Zheng et al. at the three separate ages; P4, P8 and P22 in Fig. 10.
Taken together, during development, SACs exhibit a transient bursting activity which is lost upon maturation. In the context of our bifurcation analysis, this process is interpreted as follows: immature SACs rest state lies close to a bifurcation point in order to facilitate the bursting activity which is essential for the generation of retinal waves and the shaping of the visual system. Upon maturation and once this functionality is no further useful for the retinal circuitry, we predict that SACs rest state is naturally driven towards the stable branch (further left in Fig. 2A), away from the bifurcation point, where bursting generation is no longer possible (we remind that noise could drive the cell to cross the bifurcation but, the further away, the less biophysically plausible this scenario is regarding the noise amplitude). In our model, this essentially corresponds to adding a hyperpolarizing (negative) current as an input. For identifying the source of such current, we propose a conjecture on the possible role of the Kv 3 potassium channels during development, based on our analysis, which could reveal a mechanism of how SACs lose their excitability.

A conjecture on the role of the Kv 3 channels in the loss of sACs excitability upon maturation.
We predict that bursting occurs in immature SACs, during a certain age window while fast potassium channels are under-expressed. Upon development, we make the hypothesis that potassium channels increase their expression and spontaneous bursting stops when the hyperpolarizing K + dominate the dynamics over the depolarizing Ca 2+ channels. Also, we suggest that upon maturation, the half activation potential V 3 decreases in order to provide a saturation in mature SACs activity as shown in 20 . This theoretical prediction could be experimentally tested, by following the expression of the fast K + channels in SACs during development and showing a potential increase of their activity along maturation (see discussion section). Moreover, we propose that the fast voltage gated K + channels involved is the Kv 3 family. More specifically, two subtypes of fast voltage gated K + channels have been identified 19,20 , both belonging to the Kv 3 family: i) I K delayed rectifier currents, sensitive to TEA (tetraethylammonium) emitted by subunits Kv 3 .1 and Kv 3 .2 and ii) I A A-type currents, sensitive to 4-amino-pyridine (4AP), emitted by subunits Kv 3 .3 and Kv 3 .4 19 which are not sensitive to TEA. The apparition and evolution of such channels in SACs during development has not been studied yet. However, in the experiment of 6 , the application of TEA in immature SACs (rabbit ~ P8) was crucial to the re-initiation of bursting activity. This finding implies that, during development TEA sensitive K + channels are already expressed in SACs. On this basis, we propose that the specific type of K + channels responsible for the spontaneous bursting activity during development and eventually its loss upon maturation, is the I K rectifier currents of the To conclude, we propose that, during early development (<P6), Kv 3 channels are under-expressed, allowing a competition between inhibition (K + ) and excitation (Ca 2+ ), leading SACs to burst (see Fig. 10, blue trace). Upon maturation, the expression of these channels could evolve increasingly, leading to stronger inhibition, dominating fully the competition of inhibitory/excitatory channels, by suppressing oscillations completely. Therefore, we suggest that the evolution of the expression of fast Kv 3 channels could be part of the transient process that leads to a complete loss of excitability of mature SAC. Our theoretical results indicate that the level of the expression of Kv 3 channels could increase gradually upon maturation, which essentially means that the conductance of these channels would increase as well as other physiological properties of the channels such as kinetics parameters like the half-activation potential (fixing the characteristic activation sigmoid of the channel). This type of experiment could elucidate the exact role of Kv 3 channels in the intrinsic properties of the excitability of SACs.
Robustness with respect to parameters variations. We now demonstrate the robustness of this analysis to variations of our main parameters. We already showed the effect of varying I tot and g K in Fig. 9. As bursting is a competition between calcium and potassium dynamics, we have also checked the structure of bifurcations when varying the calcium conductance g C and the potassium conductance g K around the values we have fixed in our model. We have used different software (XPPAUT 34 , MATCONT 29 , PyCont 35 ) and we were not able to get all bifurcations with them (especially homoclinic bifurcations). So we obtained the bifurcation diagram by a direct inspection of the phase portraits with respect to parameters variations. The points shown in the bifurcation maps are sampling points whereas continuous lines correspond to fit. We note that Hopf and Saddle Node bifurcations were obtained identically using the bifurcation analysis software quoted above.
The resulting bifurcation diagrams, in the two regime of bursting, are drawn in Fig. 12 (dynamically driven bursting) and 13 (noise driven bursting). In the noise driven bursting regime there is a region (G in the Fig. 13), delimited by two homoclinic bifurcations lines, where the cell is bistable without fast oscillations. Note that this region does not exist in the dynamically driven bursting regime, hence Figs 12 and 13 have a different x axis range. (The goal is to enlarge region G, which does not appear in Fig. 6).

Discussion
In this paper we have proposed a model for spontaneous bursting of SACs, one of the key ingredients for the initiation of stage II retinal waves. To our best knowledge, this is the first model able to characterize SACs dynamics in great detail and to show that biophysical properties of SACs are sufficient to explain intrinsic bursting and refractoriness in these cells. The model can reproduce spontaneously occurring fast oscillatory activity followed by long refractory periods in individual SACs. It is explained by bifurcations arising purely from intrinsic excitable membrane properties. SACs can spontaneously switch to strong oscillatory activity followed by refractoriness, reflecting what has been reported experimentally in immature SACs during the period of retinal waves. The bursting frequency is determined by a Ca-activated K conductance responsible for the sAHP. This led us to propose several conjectures and possible experiments, directly inspired from the model analysis.
We would like now to develop other aspects, not considered in the main text.
Biophysical simplifications. Although our model is able to explain fairly well and reproduce several aspects in the physiology of bursting SACs, as well as to propose new experimental predictions it is yet simplified with respect to the real cellular processes in SACs, lacking potentially certain biological mechanisms with important consequences. Especially, there are several important ion channels (e.g. I Na ) 36 and transporters in immature SACs that are not considered in our mathematical description. These currents make changes in the intracellular ion concentrations (ex Cl − concentration 37 ) and make changes in the reversal potentials (e.g. V L ) of the ion channel currents. A change in V L results in a change of the effective current g L V L in eq. (1) and therefore can be absorbed in I ext . Here, we justify further our approximations: • The role of sodium channels. In 6 , it has been shown that voltage-gated Na + channels do not contribute to the bursting activity of immature SACs, since upon the application of tetrodoxin (TTX), their activity remained unaltered. Note that the paper 36 (10 years anterior to 6 ) mention the role of Na + but it seems only responsible for the initial peak of voltage arising when the burst starts. This point is confirmed in 6 , see Fig. 4 of this paper. We have therefore decided to remove Na + channels from our modeling. Eliminating Na + channels allows us to consider a more simplified model than Hodgkin-Huxley, yet a canonical model in computational neuroscience, the Morris-Lecar model. • On calcium dynamics. In 6,12 , the authors showed that the long refractory period of SACs is controlled by Ca 2+ -gated K + channels, probably SK-like, as those observed in pyramidal cells in hippocampal neurons studied in 26 . Based on these experimental observations, we have proposed a model of sAHP dynamics. Although our modeling is widely inspired from 12 , we get deeper in the biophysical description of Calcium dynamics, based on the seminal work by M. Graubner and collaborators 27,38 (see Methods). Additionally, the parameters were tuned to match the SK channels experimental characterization made by Abel et al. 26 . Yet, we have simplified the dynamics of NCX exchangers and we neglected the ionic pumps currents. These simplifications are detailed in Methods.
The Interburst variation and the interpretation of I ext . Spontaneous bursting activity in immature SAC has been observed consistently across species i.e. mice, rabbits, chicks etc., but with different characteristic bursting periods. In Fig. 6, we show the relationship of the average bursting period τ IBI with respect to the value of the tunable current I ext whose effect could be translated on how an effective inhibition affects the characteristics of the bursting activity and even stops it completely. The source of such an inhibition could be linked either with intrinsic cell transient properties or inhibitory inputs from other layers of the retina upon maturation. Although it is unlikely that changes in intrinsic properties of SACs would occur over the short period of stage II, such changes could be observed across species. As the interburst intervals (measured in retinal ganglion cells) increase and then decrease (in mouse) during Stage II waves (see 3 ) the effective inhibition current I ext is more likely arising from changes in synaptic inhibition (e.g. GABA A maturation or changes in the inner chloride concentration).

How identical cells can display variability in behaviour.
One of the key result of this study is that SACs display bursting because they are close to a bifurcation point. As a consequence, the characteristics of bursting (IBI, amplitude) vary strongly upon tiny variations of physiological parameters such as … g g I V , , , K C ext 3 and noise can induce huge fluctuations in their bursting period. In a network of SACs those parameters may fluctuate. Especially, the presence of acetylcholine inducing a cholinergic current impacts the bursting period as was experimentally observed by Zheng et al. in 6,14 and confirmed by our model 22,23 . Therefore, a network of identical coupled SACs is expected to display a huge variability: SACs are not identical bursters. This variability has been invoked and modeled in a paper by Ford and Feller 21 where they show how it impacts the retinal waves dynamics and provide more realistic results than an homogeneous model. In their paper, inhomogeneity is imposed ad hoc (Gaussian distribution of IBIs) and does not evolve in time. In contrast we propose here that (i) inhomogeneity purely results from the biophysical mechanisms leading SACs to burst; (ii) this heterogeneity evolves in space and time, depending on the collective dynamics (e.g. acetylcholine concentration varies in space and time and induces a variability in a given cell bursting period). In particular, one could have two populations of cells: one population acting as dynamically driven bursters with a well defined bursting period, and one population acting as noise driven bursters with a large variation of the period. Those populations would act in a different way on the collective dynamics of waves (see next section).

How inhomogeneous bursters can produce retinal waves?
To our best knowledge, this model is the first able to characterize SACs dynamics in great detail and to show that biophysical properties of SACs are sufficient to explain intrinsic bursting and refractoriness in these cells. Yet, this work was concerned with individual neuron dynamics a crucial activity for the emergence of a propagating wave activity. A next step therefore consists of studying a network of coupled SACs. The crucial point is to understand how a network of coupled inhomogeneous bursters can achieve the third ingredient necessary to initiate retinal wave: synchrony. In the retina, at stage II, SACs are coupled via cholinergic synapses, whose strength evolve during development 6 . Modeling the cholinergic coupling, we observed indeed that SACs are able to achieve local synchrony and generate waves for sufficiently large value of the cholinergic conductance, compatible with experiments 6 . This is further developed in a forthcoming paper but we give here the main ideas. Local synchrony results from the cholinergic coupling via the classical mechanism of non linear synchronization between inhomogeneous oscillators, called "Arnold' tongues" 39 . For a sufficiently large cholinergic conductance, compatible with experiments, this synchrony leads to a wave propagation. This propagation can be described, in our model, via a transport equation allowing to compute the wave characteristics (speed, size, duration) as well as the effect of sAHP on dynamics. Because of its slow time scale, quite slower than the characteristic time of the wave propagation, sAHP somewhat imprints the medium where the wave has to propagates leading to transient spatial structures observed in Ford et al. 40 . The full analysis of these aspects will be done in a forthcoming paper (see [22][23][24]. Pharmacological control. Although this model has many parameters, bifurcation analysis allows us to highlight several of them, g C , g K , V 3 and I ext , controlling important aspects of dynamics, in direct links with experiments. Varying those parameters have a deep impact on SACs dynamics as shown e.g. in Fig. 10. This is possible via a pharmacological control (e.g. TEA to vary K conductance, Cd to vary Ca conductance, synaptic input changing I ext …). It would be interesting to confirm experimentally our bifurcation diagram with such experiments.
In addition, this study opens the possibility that mature SACs are still potential bursters. Then, by a suitable pharmacological treatment, guided by the bifurcation analysis, they could start to burst in mature retinas. Would this restore a wave activity as recently observed in mature mice retinas 41 . This deserves however further investigations. Especially, SACs are not connected to each other in the adult retina (cholinergic synapses disappear at the end of stage II 6 ). In adults SACs make synapses with bipolar cells axons, with retinal ganglion cells, and with some other amacrines types (inhibitory) providing several potential coupling mechanism to propagate waves.
proposing new experiments towards the validation of our theoretical predictions. How to validate experimentally a bifurcation diagram for SACs dynamics? It would be very elucidating to conduct experiments in order to calibrate accurately SACs dynamical repertoire with respect to our bifurcation analysis shown in Fig. 2A. It would require to apply several consequent constant current pulses to an isolated immature SAC, within a wide range, and create an "experimental bifurcation diagram", showing the observed behaviors of the cell upon the variation of a constant current. Such an experiment, would shed light also on understanding the two possible drives for triggering SACs bursting which is a result of our modeling. Also, it would help validate our hypothesis about SACs being close to a bifurcation during development.
How to validate experimentally the possible role of K + channels on the loss of SACs excitability? One of the predictions that emerged from our model is that the Kv 3 channels, that are expressed in the SACs at adult stage 20 , may not be expressed at stage II. To test this hypothesis, immunolabeling experiments are needed in the mouse retina, using an antibody labeling the Kv channels, both in adult and immature retina in order to compare and deduct whether the localization of Kv channels is weaker or inexistent in the immature retina (P3-P10 for mice). Along with immunochemistry experiments, actual measurements of the K + currents at different stages of development (we thank one of the reviewers for this comment) would elucidate their possible transient expression during development and their role in the lose of SACs excitability. and we obtain a kinetic equation: where C is the intracellular calcium concentration. Setting τ = Note that equation (7) has a similar form to the one proposed by Hennig et al. in 12 . However, in their model our term C 4 is replaced by Binding of calmodulin to SK terminals. This corresponds to a reaction: where F is the density of free terminals, B the density of bounded terminals, P FB (P BF ) the transition rate from free to bound (bound to free).
Calling R the fraction of bounded terminals, τ = we end up with a kinetic equation: Finally, since R is the probability that a terminal is open and since 4 terminals must be open to let the SK channel open, the sAHP conductance is g sAHP R 4 .
Note that equation (6), is similar to Eq. (3) in 12 , with a remarkable difference: in our model there is no direct dependence on calcium concentration whereas the term αC in 12 corresponds to a direct binding of Ca 2+ to a terminal. Note that, taking the quite large value of the parameter α (~2400) considered by these authors, their equation is essentially equivalent to τ α α α with a steady state = α α + R C C 1 very close to 1 whenever αC is quite larger than 1. In this case, S plays essentially no role.
Calcium concentration. Both variables R and S are driven by intracellular Ca 2+ concentration dynamics, given by: 1 is the surface to volume ratio that accounts for the localization of the channels at the surface of the membrane, n Ca = 2 is the calcium valence and F = 96500 C mol − The first term in Eq. (13), r I S C corresponds to an increase of internal Ca 2+ concentration upon calcium influx (current I C (V)) generated by spikes or, in experiments, by voltage clamp. As V Ca = 50 mV this current is positive The second term is −rJ X (C), where is the efflux current density through sodium-calcium exchanger (NCX). It corresponds to an outward current through NCX exchangers, contributing to restoring the initial Ca 2+ concentration. Here, ρ X is the surface density of the NCX membrane proteins. We take ρ X = 100 μm −2 = 10 12 dm −2 (from 27 ). I X = 4.8 × 10 −19 Cms −1 = 4.8 × 10 − 16 Cs −1 is the maximum ionic flux through a single NCX channel. This corresponds to 3+ charges (1.5 Ca 2+ ) per ms; H X = 1.8 μM = 1.8 × 10 −6 M is the half activation concentration 38 . Also, in (13), J p (C) is the current density of Ca 2+ pumps. We shall neglect this term from now on. Finally, L is the leakage surface current density representing the residual conductivity of the plasma membrane. We have not been able to find its value in the literature.
To summarize Eq. (13) becomes: We also assume that NCX current term is approximated by a linear term. This is valid if one assumes that calcium concentration  C H X which is the case in our simulations. In order to match (15) with the form (12) ((8) in the text), we set: Ca Ca τ = . C rL (18) Ca 0 How does the average interburst depend on the parameter I ext . In this section we analyze how the interburst interval τ IBI depends on the current I ext . In particular, we justify Eq. (9). As explained in the text there are two distinct regimes where the interburst interval τ IBI is constrained by different factors. In the noise driven bursting regime ( < I I ext S N 1 ), the stochastic dynamics around the rest state is described, with a good approximation, by a Ornstein-Uhlenbeck process (upon neglecting the non linearities coming from g C and g K , which are small in this range of membrane potential). Therefore, fluctuations of V are Gaussian with a mean V rest and a variance proportional to σ, the noise intensity. Thus, the probability to cross the bifurcation threshold and enter into the bursting regime can be easily computed. This is a sigmoid with a slope proportional to σ 1 . As a consequence, bursting can take place when < I I ext S N 1 , because of noise, but the probability to burst decreases rapidly when I ext decreases, after a threshold value depending on σ. Formally, the bifurcation value I SN 1 is shifted to I c , explaining this term in Eq. (9). For I ext < I c the probability to induce bursting by noise is so small that we consider it vanishes. Here τ IBI is not defined anymore but we set it to zero in Eq. (9) by convention.
In the dynamically driven bursting regime ( > I I ext S N 1 ), τ IBI is constrained by two times: (i) the time t 1 to get in a neighborhood of the bifurcation point V SN 1 in Fig. 3 where bursting starts (essentially determined by the sAHP) and (ii) the time t 2 to leave this neighborhood. When I ext is quite larger than I SN 1 , τ IBI is largely dominated by t 1 .
Here, τ IBI decreases slowly as I ext increases. Indeed, the higher I ext the more neurons are prone to bursting. When I ext becomes too large, I sAHP is not enough to stop the oscillatory phase and τ IBI = 0.
Close to I SN 1 , we can use the normal form of the saddle-node bifurcation for the variable V: where V s is the rest state at the bifurcation point and α depends on the other parameters of the model. In this approximation, there are two branches of fixed points, existing only when (1 )  Table 1. Having integrated over all the surface of the membrane we omit the surface units.
Calibrating parameters from experiments. All parameters values are calibrated with respect to biophysics, found in the literature or fitted from experimental curves in 14,26 and 6 . Morris-lecar tuning parameters V 1 , V 2 , V 3 , V 4 were calibrated (see Fig. 8), so as to reproduce the experiment of 6 (Fig. 4a), where the authors investigate the ionic mechanisms of the fast oscillations. Note that the bursting regime is robust to (small) variations of these parameters (results not shown). We tuned the sAHP parameters taking into account the analogy with SK channels studied in 26 (fit not shown). Also, we note that the intensity of sAHP observed by Abel et al. in pyramidal neurons (of order 150 pA) is quite bigger than in stage II SAC. In our model, this means a lower sAHP conductance g sAHP (g sAHP = 2 nS).
Auxiliary functions. The dimensionless auxiliary functions involved in the dynamical equations appearing in the model definition are:   3 4 Parameters. The parameters used in the model are displayed in Table 2.

Physical quantity Units
Comparison with existing models. In this section we shortly revisit models of SACs activity in the stage II and compare them to our model, (see [11][12][13]40,42,43 , for a review see 10 ). We would first like to remark that these existing models are devoted to describe wave activity and do not focus on thoroughly describing individual SACs dynamics. Especially, none of the models we know describe the biophysical mechanisms of SACs bursting activity and the role played by biophysical parameters. Instead their focus was more on having a relatively simple description of the cell activity with a minimal set of tunable parameters (a notable example is Butts et al. model 43 which has two free parameters governing the waves properties). Finally, our model is the first one, able to characterize SACs dynamics in great detail and to show that biophysical properties of SACs are sufficient to explain intrinsic bursting and refractoriness in these cells.
The closest model to ours has been proposed by M. Hennig and collaborators (refered to as "Hennig model" 12 , see also the extension by Ford and Feller and the recent paper of Xu et al. 44 ). Actually, our model has been widely inspired by this work with several notable differences. As exposed in the Methods section our biophysical analysis of sAHP dynamics leads to equations and parameters values departing from Hennig model. Additionally, Hennig model does not consider a fast potassium dynamics and there is no fast oscillation. The mechanism that mimics SACs bursting is a switch from low membrane potential level to high one. This switch is determined by an exogenous shot noise i.e. a voltage dependent rate modulated Poisson process with a slow decay. This activity is maintained long enough so that sAHP can be activated, enabling the cell to return to rest. In our model shot noise is not necessary to trigger activity. Instead, a cell can spontaneously switch to the bursting state, where it stays until the sAHP produced by its activity leads it back to the rest state. By "spontaneous" we mean literally "happening or done in a sudden way, without any planning or without being forced/without premeditation". In our model this sudden switch is a bifurcation, induced by the mere cells dynamics. The presence of a fast (Brownian) noise facilitates this transition, but, there is no need for a shot noise. The cell stays in the bursting state by its mere dynamics, even when it is isolated.

Parameter
Physical value C m 22 pF g L 2 nS g C nS [3,20] g K nS [1,20] Table 2. Range of values for the parameters used in the paper.