Bruno Cessac

Dora Matzakou-Karvouniari

The non linear dynamics of retinal waves

Keywords: Retina development, neuroscience modelling, nonlinear dynamics, bifurcations

We investigate the dynamics of stage II retinal waves via a dynamical system, grounded on biophysics, and analysed with bifurcation theory. We show how the nonlinear cells coupling and bifurcation structure explain how waves start, propagate, interact and stop. Especially, we analyse how the existence of a small region in parameter space, where dynamics returns in a recurrent way, gives rise to a very rich dynamics. In this context, we propose a non linear transport equation characterizing the waves propagation and interaction.

Introduction

The visual system is an important part of the central nervous system. Although its functionality seems effortless, it carries out complex tasks including the reception of light and the formation of visual representations, the identification and categorization of visual objects, computing distances to and between 5 objects and guiding body movements in relation to the environment [76,[START_REF] Chalupa | The visual neurosciences[END_REF][START_REF] Daw | How Vision Works: The Physiological Mechanisms Behind What We See[END_REF]. How the visual system acquires such abilities during development, pre-and postnatal, is a fascinating question.

Vision starts at the retina, a light-sensitive tissue at the back of the eye that covers about 65 percent of its interior surface [START_REF] Besharse | The Retina and its Disorders[END_REF]. The overwhelming capac-10 ity of the retina to convert complex visual scenes into spike trains -that send information to the visual cortex -is largely due to its layered structure and to retical analysis of the model's structural stability (stability of a behaviour with respect to parameters variations). This is actually a salient question as physiological parameters such as the reversal potential for GABA receptors or the efficacy of cholinergic synapses [START_REF] Zheng | A developmental switch in the excitability and function of the starburst network in the mammalian retina[END_REF] are known to change during development, thereby inducing potential changes in wave dynamics.

Especially, in [START_REF] Zheng | A developmental switch in the excitability and function of the starburst network in the mammalian retina[END_REF], Zheng et al. have experimentally shown that the spontaneous stage II retinal waves are mediated by a transient network of SACs, connected through excitatory cholinergic connections [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF], which are formed only during a developmental window up to their complete disappearance. Especially, the intensity of the acetylcholine coupling is monotonously decreasing with time (see Fig. 3B in [START_REF] Zheng | A developmental switch in the excitability and function of the starburst network in the mammalian retina[END_REF] and table 1 below). How do stage II retinal waves dynamics depend on this parameter? How do their characteristics (e.g. size, duration, speed) evolve when acetylcholine coupling decreases? What could be the impact on visual system development? In the present paper, we investigate the question of the dynamic origin and structure of stage II retinal waves, their distribution, and the stability of their behaviour, during the stage II developmental window, upon a decreasing acetylcholine coupling compatible with experimental findings. This paper is the continuation of a study starting with [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF] and fully exposed in the PhD thesis [START_REF] Karvouniari | Retinal waves: theory, numerics, experiments[END_REF] aiming to: [START_REF] Abel | Relationships between intracellularneocortical pyramidal neurons calcium and afterhyperpolarizations in neocortical pyramidal neurons[END_REF] propose a detailed biophysical modeling of the cellular mechanisms of the spontaneous activity in immature SACs, addressing questions directly accessible and linked to pharmacological manipulations on retinal waves; (2) develop a mathematical analysis of this model, based on dynamical systems and bifurcation theory so as to propose generic mechanisms for wave initiation, propagation and stop. In addition, such analysis allows to question the structural stability of the model; (3) validate the working hypothesis that SACs produce a wide spectrum of waves duration, size, interwave period, across species, thanks to a single fact: they are close to specific bifurcations generating what Izhikevich calls Type 1 excitability [START_REF] Izhikevich | Neural excitability, spiking, and bursting[END_REF][START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability And Bursting[END_REF]. Although point [START_REF] Alili | Representations of the first hitting time density of an ornstein-uhlenbeck process[END_REF] and (3) might look contradictory we will show here that this is not the case.

The paper [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF] was devoted to the study of single, isolated cells dynamics.

The detailed modelling of individual SACs dynamics as autonomous, rhythmic bursters and the mathematical analysis of our dynamical system using bifurcation theory helped us identify the key parameters which control bursting in immature SACs. Especially, we exhibited that few biophysical parameters regulating calcium and potassium activity control bursting and we proposed a testable experimental prediction on the role of voltage-dependent potassium channels on the transitory excitability properties of SACs along development.

We also proposed an explanation on how SACs can exhibit a large variability in their bursting periods across different species , as observed experimentally, yet based on a simple, unique, mechanism.

In the present paper, extending the single neuron dynamics, we model in detail the mutual cholinergic synaptic connections between SACs, ending up exploring the mechanisms of SACs synchronization and waves. In section 2, we give a brief account of retinal waves properties, for non expert readers, before introducing a multi-dimensional model of SACs featuring their intrinsic dynamics (bursting and hyperpolarization) and their non linear coupling via acetylcholine. We argue that the dynamics of SACs waves is essentially controlled by two parameters slowly evolving in time: one, G A , controlling the excitatory cholinergic cells coupling, and the other, G S , controlling cell's hyperpolarisation and refractoriness. Thanks to a bifurcation diagram in the space {G S , G A }, we analyse, in section 3, how waves start, propagate, interact and stop. Mostly, the variety observed in waves dynamics comes from the fact that they start in a tiny region in the space {G S , G A } delimited by bifurcation lines, where waves initiation is quite sensitive to perturbations such as noise, which is a feature of type 1 excitable systems [START_REF] Izhikevich | Neural excitability, spiking, and bursting[END_REF][START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability And Bursting[END_REF]. Although this region is tight, the slow dynamics returns to it in a recurrent way, regenerating the potentiality to trigger new waves sensitive to perturbations with a wide variation in waves characteristics. In addition, this scenario holds on an interval of acetylcholine coupling, explaining the apparent contradiction between a variability induced by closeness to bifurcations and stability to variations in developmental parameters. On these tracks we numerically investigate, in section 4, the effect of varying acetylcholine coupling during development, in a range of parameters extracted from the experimental literature. Simulations are done in a 1-dimensional lattice, with nearest neighbours interactions. To extend our analysis to larger dimensions we derive, in section 5, transport equations for G S , G A , now considered as propagating fields shaping waves dynamics. From this, we are able to compute the wave speed as a function of acetylcholine coupling, as well as to show the existence of a critical value of this coupling, below which no wave can propagate. As we argue, these transport equations bare interesting analogies with the Kardar-Parisi-Zhang (KPZ) equations of interface growth [START_REF] Kardar | Dynamic scaling of growing interfaces[END_REF] on one hand, and Self-Organized Criticality [START_REF] Bak | Self-organized criticality: An explanation of 1/f noise[END_REF] on the other hand, opening up perspectives for future research.

Stage II retinal waves model

The biophysics of retinal waves

Retinal waves are bursts of activity occurring spontaneously in the developing retina of vertebrate species, contributing to the shaping of the visual system organization [START_REF] Wong | Transient period of correlated bursting activity during development of the mammalian retina[END_REF][START_REF] Firth | Retinal waves: mechanisms and function in visual system development[END_REF][START_REF] Sernagor | Chapter 49 -retinal waves: Underlying cellular mechanisms and theoretical considerations[END_REF][START_REF] Ford | Assembly and disassembly of a retinal cholinergic network[END_REF]. They are characterized by localized groups of neurons becoming simultaneously active, initiated at random points and propagating at speeds ranging from 100 µm/s (mouse, [START_REF] Singer | Potentiation of L-Type Calcium Channels Reveals Nonsynaptic Mechanisms that Correlate Spontaneous Activity in the Developing Mammalian Retina[END_REF], [START_REF] Maccione | Following the ontogeny of retinal waves: pan-retinal recordings of population dynamics in the neonatal mouse[END_REF]) up to 400 µm/s (chick, [START_REF] Sernagor | Differential effects of acetylcholine and glutamate blockade on the spa-tiotemporal dynamics of retinal waves[END_REF]), with changing boundaries, dependent on local refractoriness [START_REF] Ford | Assembly and disassembly of a retinal cholinergic network[END_REF][START_REF] Ford | Cellular mechanisms underlying spatiotemporal features of cholinergic retinal waves[END_REF]. This activity, slowly spreading across the retina, is an inherent property of the retinal network [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF]. More precisely, the generation of waves requires three conditions [START_REF] Gjorgjieva | Modeling developmental patterns of spontaneous activity[END_REF][START_REF] Ford | Assembly and disassembly of a retinal cholinergic network[END_REF]: (C1) A source of depolarization for wave initiation ("How do waves start?").

Given that there is no external input (e.g. from visual stimulation in the early retina), there must be some intrinsic mechanism by which neurons become active1 .

(C2) A network of excitatory interactions for propagation ("How do waves propagate?"). Once some neurons become spontaneously active, how do they excite neighboring neurons? (C3) A source of inhibition that limits the spatial extent of waves and dictates the minimum interval between them ("How do waves stop ? ").

Wave activity begins in the early development, long before the retina is responsive to light. It emerges due to several biophysical mechanisms, which change during development, dividing retinal waves maturation into 3 stages (I, II, III) [START_REF] Sernagor | Chapter 49 -retinal waves: Underlying cellular mechanisms and theoretical considerations[END_REF]. Each stage, mostly studied in mammals, is characterized by a certain type of network interaction (condition C2): gap junctions for stage I; cholinergic transmission for stage II; and glutamatergic transmission for stage III. In this work, we focus on stage II.

During this period, the principal mechanism for transmission is due to the neurotransmitter acetylcholine (Ach) with nicotinic receptors [START_REF] Feller | Dynamic processes shape spatiotemporal properties of retinal waves[END_REF][START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF][START_REF] Sernagor | Chapter 49 -retinal waves: Underlying cellular mechanisms and theoretical considerations[END_REF]. The first functional cholinergic connections are formed around birth, firstly at the level of Starburst Amacrine Cells (SACs). The emergence of waves depends on cellular mechanisms studied by [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF] (for rabbits), where it is found that Starburst Amacrine Cells emit spontaneous intrinsic calcium bursts (condition C1). When bursting a SAC emits acetylcholine thereby increasing the level of excitability of its neighbours. When the bursts occurring in a neighbourhood of a quiet SACs are synchronized, they eventually lead it to burst, inducing a wave propagation. Stage II waves are therefore spatiotemporal phenomena resulting from the spatial coupling of local bursters (SACs) via acetylcholine (condition C2). In addition, a strong after hyperpolarization current (sAHP) induces a long refractory period for the recently active neurons, preventing the propagation of a new wave for a period of order of a minute, in the region where a wave has recently spread [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF]. More generally, sAHP plays a prominent role

has not yet been identified, act as pacemaker triggering waves [START_REF] De Montigny | Evidence for novel transient clusters of cholinergic ganglion cells in the neonatal mouse retina[END_REF]. We shall not consider this aspect in our work although it could be included in the proposed formalism.

in shaping waves periodicity and boundaries (condition C3), as we discuss in this paper.

The model

Among the 3 stages of retinal waves, stage II is the one which received the most attention from modellers [START_REF] Feller | Dynamic processes shape spatiotemporal properties of retinal waves[END_REF][START_REF] Butts | Retinal waves are governed by collective network properties[END_REF][START_REF] Godfrey | Retinal wave behavior through activitydependent refractory periods[END_REF][START_REF] Godfrey | Theoretical models of spontaneous activity generation and propagation in the developing retina[END_REF][START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF][START_REF] Gjorgjieva | Modeling developmental patterns of spontaneous activity[END_REF][START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF] (see [START_REF] Kähne | Gap junctions set the speed and nucleation rate of stage i retinal waves[END_REF] for a recent model on stage I, and [START_REF] Choi | Intrinsic bursting of aii amacrine cells underlies oscillations in the rd1 mouse retina[END_REF] for a model of stage III). The model we propose is inspired by the work of [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF] and [START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF], with strong modifications, in the model itself and in the methods used for study. It was introduced in [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF] at the level of single cells (no network) where it was thoroughly justified on biophysical grounds, leading to experimental predictions. The approach to study its dynamics, based on a bifurcation analysis of an individual SAC dynamics was original too. Here, we extend the study to a network of SACs, with an analysis still grounded as well on bifurcation theory. In this section, we define the equations ruling the model's dynamics.

Note that SACs display a regular tiling of the retina with a disk-shaped dendritic tree [START_REF] Tauchi | The shape and arrangement of the cholinergic neurons in the rabbit retina[END_REF]. Their average distance is a 50 µm. We thus define a network of Starburst Amacrine Cells distributed on a regular lattice in Z d , d = 1, 2, with lattice spacing a, where sites/neurons are labelled with an index i = 1 . . . N . The state variables characterizing neuron i are: V i (t), the membrane potential, N i (t), the gating variable for fast K + channels, R i (t) and S i (t), the gating variables for slow Ca 2+ -gated K + channels, C i (t), the intracellular Ca 2+ concentration, A i (t), the extracellular acetylcholine concentration emitted by neuron i. We note B i the set of SACs which are in synaptic contact with SAC i. A SAC is not connected to itself, (i ∈ B i ). The structure of neighbours is nearest neighbours. This is mainly because it allows to write explicit transport equations in term of a Laplacian (see section 5). More complex structures have been considered in [START_REF] Karvouniari | Retinal waves: theory, numerics, experiments[END_REF] or [START_REF] Karvouniari | Mathematical and experimental studies on retinal waves[END_REF] for this model. We come back to this point in the discussion section.

The parameters values are carefully calibrated with respect to biophysics, found in the literature or fitted from experimental curves in [START_REF] Abel | Relationships between intracellularneocortical pyramidal neurons calcium and afterhyperpolarizations in neocortical pyramidal neurons[END_REF], [START_REF] Zheng | A developmental switch in the excitability and function of the starburst network in the mammalian retina[END_REF] and [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF] (see [START_REF] Karvouniari | Retinal waves: theory, numerics, experiments[END_REF] for detail on the fitting procedures). Parameters values and units are given in appendix 6.3.

Dynamics has several time scales, from fast (a few milliseconds) to slow (a few seconds) to very slow (about one minute). These multiple scales are at the source of the dynamic complexity and the core of the present analysis, so we present the equations ruling the evolution of SACs taking into account this time scale separation. First, the membrane potential of SAC i obeys:

C m dVi dt = I Li + I Ci + I Ki + I Si + I Ai ; (1)
where C m is the membrane capacitance,

I Li = -g L ( V i -V L ) is a leak current,
with leak conductance g L and leak reversal potential V L . The current:

I Ci = -g C M ∞ (V i )(V i -V C ), (2)
is a calcium current of the Morris-Lecar form [START_REF] Morris | Voltage oscillations in the barnacle giant muscle fiber[END_REF], with maximal calcium con-

ductance g C , calcium reversal potential V C and M ∞ (V ) = 1 2 [1 + tanh( V -V (1) V (2) )]. The fast potassium current I Ki takes the form I Ki (V i , N i ) = -g K N i (V i -V K ),
where the evolution of the voltage-gated K + channel gating variable N i (t) is given by:

τ N dN i dt = Λ(V i ) ( N ∞ (V i ) -N i ) , (3) with
Λ(V i ) = cosh( Vi-V (3) 2V (4) ), and N ∞ (V i ) = 1 2 [1 + tanh( V -V (3) V (4) )].
Here τ N , the characteristic time of the activation variable N , is of order 5 ms. I Ki is therefore quite fast, in contrast to the original Morris-Lecar paper, calibrated in order to capture the frequency of the fast repetitive firing of SACs, (around 20 Hz [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF]).

V (1) , V (2) , V (3) , V (4) are tuning parameters whose value is given in the appendix 6.3.

The slow After Hyperpolarization (sAHP) current

I Si = -G Si ( V i -V K ) is
a calcium-gated potassium current, whose effective conductance is:

G Si = g S R 4 i , (4)
g S being the maximal sAHP conductance, evolves very slowly (time scale of order one minute). This conductance is controlled by a cascade of mechanisms involving the entrance of calcium in the cell, and, thus, the increase of the intra-cellular calcium concentration, C i , upon membrane potential increase. Increasing C i increases the fraction, S i , of saturated calmodulin that binds to slow calcium-gated potassium receptors. This increases the fraction, R i , of bounded terminals in the corresponding potassium channels, thereby increasing G Si . There is a power 4 because 4 bound terminals are needed to open a Ca 2+gated K + channel. This large exponent, 4, is actually quite important for the waves dynamics, as developed below.

The corresponding equations are [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF]:

         τ C dCi dt = -α C H X C i + C (0) + δ C I Ci (V i ); τ S dSi dt = α S C 4 i (1 -S i ) -S i ; τ R dRi dt = α R S i (1 -R i ) -R i ; (5)
where I Ci is given by [START_REF] Alili | Representations of the first hitting time density of an ornstein-uhlenbeck process[END_REF]. The value of the tuning parameters α C , C (0) , H X , δ C , α S , α R used in this paper are given in the appendix 6.3. The time scale τ C is of order 2 s whereas τ S , τ R are of order one minute in real SACs. Actually, in the simulation displayed in section 4.1 we will consider shorter time scales, of order 10 s, to reduce the computational time.

All the terms in (1) considered up to now correspond to local, uncoupled, dynamics of SACs. In addition, SACS are coupled via the Ach current:

I Ai = -G Ai (V i -V A ), (6)
with:

G Ai = g A j∈Bi A 2 j γ A + A 2 j , (7)
where g A is the maximal Ach conductance and γ A = K 2 d , K d is the halfactivation constant of Ach production (∼ 1 -4 nM). The effective Ach conductance depends on the Ach concentration produced by SACs j connected to i. The exponent 2 comes from the physiology of nicotinic receptors: two Ach molecules have to bind to the receptor to open the corresponding ion channel (permeable to sodium and potassium, resulting in a reversal potential V A close to 0 mV). The Ach concentration A j depends on the pre-synaptic cell's voltage via the differential equation:

dA j dt = -µ A A j + β A T A (V j ), (8)
where µ A is a degradation coefficient (of order 2 s -1 ), β A is the maximal production rate of A (in nM s -1 ), and the production term is:

T A (V j ) = 1 1 + e -κ A (Vj -V A ) . (9)
The parameters µ A , κ A , β A , V A have been fitted from experiments [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF][START_REF] Karvouniari | Retinal waves: theory, numerics, experiments[END_REF].

They are uniform on the SACs population (independent of i). The set of equations ( 1)-( 9) constitutes a non linear dynamical system of dimension 7 × N on a spatial lattice of N SACs. All parameters value are reported in appendix 6.3.

Note that the sAHP conductance, G Si , depends on cell i intrinsic activity, whereas the Ach conductance, G Ai , depends on the activity of cells connected to cell i. Also note the notations: g S , g A are physiological parameters (maximal conductance), whereas G S , G A are dynamical variables depending non linearly on the system evolution. The time scales of evolution are rather different: G A has a time scale of order seconds (due to the evolution of A, eq. ( 8)), G S has a time scale of order a minute (due to the evolution of R). Finally, later in the paper (section 4.1) we are considering the evolution of g A corresponding, physiologically, to a time scale of order a day (see Tab. 1).

Bifurcation analysis

The key observation affording the model analysis is that the evolution of the effective conductance G Si (sAHP) and G Ai (Ach) of a SAC, i, is slow compared to the fast dynamics of V i , N i . The idea is thus to consider G Si , G Ai as slow parameters tuning the evolution of the fast dynamics (1), [START_REF] Arinyo-I Prats | Traveling pulses in type-i excitable media[END_REF]. In addition, we will see that the pattern of slow-time evolution of G Si , G Ai during a wave is pretty similar for all cells. So, the central idea is to consider first equations ( 1)- [START_REF] Arinyo-I Prats | Traveling pulses in type-i excitable media[END_REF] with tunable parameters, the sAHP conductance, G S , and the Ach conductance, G A , independent of the cell. The core of this analysis is then displayed in the bifurcation diagram of Fig. 1. We explain the diagram in this section. Later, we will let G S , G A slowly evolve in time for each cell, so as to describe the dynamics of waves initiation, propagation and stop.

Let us now comment about Fig. [START_REF] Arinyo-I Prats | Traveling pulses in type-i excitable media[END_REF][START_REF] Moreno-Speigelberg | Bifurcation structure of traveling pulses in type-i excitable media[END_REF] dealing with traveling pulses in Type 1 excitable media with some important differences though: (1) the presence of noise in our model; [START_REF] Alili | Representations of the first hitting time density of an ornstein-uhlenbeck process[END_REF] the fact that our parameters G S , G A are slowly evolving according to the fast dynamics;

(3) the presence of long lasting refractory regions shaping waves propagation.

Although it could actually be interesting to see what happens in our model when entering this region, note that this would require values of g A about 10

times higher than what we have inferred from experimental data (see Table 1).

Excluding this region, Fig. When varying G S , G A the transition between regions correspond to the bifurcation lines displayed in Fig. 1 top. "SN" refers to Saddle-Node, "Hc" to homoclinic, "Hopf" to Hopf bifurcation, and "SNIC" to Saddle-Node on an Invariant Cycle (see [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields[END_REF] or [START_REF] Ermentrout | Foundations of Mathematical Neuroscience[END_REF] for a classification of bifurcations). One goes from region A to region B by the saddle-node bifurcation SN 1 where the Sd and UF of region A coalesce and disappear. For the transition from A to C there are two situations: (1), crossing region D. Here, the stable and unstable manifold of Sd give rise to a limit cycle by homoclinisation (H c ), then, Si and Sd coalesce by the saddle-node bifurcation SN 2 , leaving a unique attractor, the stable periodic orbit of region C. (2) Direct transition from A to C by a SNIC bifurcation (green line in Fig. 1). Movies illustrating these bifurcations are available on the web page https://team.inria.fr/biovision/bifurcations-map-ofthe-retinal-waves-model/.

The equation characterizing the bifurcation line of the transition between region D and C via the SN 2 bifurcation will be used later. It is given by [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF]:

-G A ( V Si -V A ) -G S ( V Si -V K ) = I SN , (10)
where V Si < 0 is the rest voltage on the sink in region D. V Si is actually a "SN" stands for "Saddle-Node bifurcation", "Hc" for "Homoclinic bifurcation", "Hopf" for "Hopf bifurcation", "SNIC" for "Saddle-Node on an Invariant Cycle","SNH" for Saddle-Node-Homoclinic point (also called saddle-node separatrix loop [START_REF] Schecter | The saddle-node separatrix-loop bifurcation[END_REF], SNSL), "BT" stands for Bogdanov-Takens. At the SNH codim-2 bifurcation point, there is a Hc bifurcation line connecting to the BT point with a SNIC emerging from the opposite side of the SNH in which the SN2 and Hc are tangent. At the BT point, there is, on one side the SN1 line and, on the other side the Hc, the SN1 and a Hopf line, tangent at the BT point (see https://team.inria.fr/biovision/bifurcations-map-of-the-retinal-wavesmodel/ for a zoom on this region, explanations and movies). Inset: Zoom on region D, in log scale.

The Hc and SN2 line meet tangentially at the SNH point and a SNIC line emerges on the other side. 2. Phase portraits of regions A, B, D, C with nullclines (N V , in black, is the V nullcline, N N , in grey, is the N nullcline), stable manifold of hyperbolic fixed points (Ws in cyan), unstable manifold of the saddle-fixed points (Wu in red), stable periodic orbit (in blue). "Si" means "Sink", "Sd" means "Saddle", "UF" means "Unstable Focus", "SPo" means "Stable Periodic orbit". function of G S , G A given by eq. ( 11) below. The current I SN ∼ 0.3 pA [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF]. A similar equation holds for the homoclinic transition where I SN is replaced by I Hc , with a value quite close to I SN .

In region D there is bi-stability between a rest state and an excited oscillatory state (in contrast to region A where there is bistability between two stationary states, one with low and the other with high voltage). Noise allows to switch from rest to oscillatory state, as developed below. This corresponds to what Izhikevich call Type 1 excitability [START_REF] Izhikevich | Neural excitability, spiking, and bursting[END_REF][START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability And Bursting[END_REF]. Now, in the bifurcation diagram, which has been represented in large range of parameters values G S , G A for completeness, region D looks quite small and may appear as probably irrelevant, but we show that it is not. Recall first that G S corresponds, in the model, to 4 where g S is the maximal sAHP conductance, whereas R ∈ [0, 1] is a probability. For our parameters value, in the absence of waves, cells are at rest with R ∼ 0.3 (see Fig. 2), so that G S ∼ 0.0162 nS for g S = 2 nS, and G S ∼ 0.081 nS for g S = 10 nS. In addition, at rest too, (V ∼ -70 mV), T A (V ) ∼ 0 (eq. ( 8)), so that the Ach production is essentially zero, therefore G A ∼ 0. This implies that, at rest, cells are in region D. More generally, we argue below that there is long period, at the end of sAHP time course, where the cell lies in region D (see also Fig. 6). As we will see, this has a big impact on wave initiation and propagation.

g S R

The mechanism of bursting

Bursting, the succession of fast oscillations and rest state, is a well known phenomenon, widely studied in the literature of mathematical neuroscience [START_REF] Izhikevich | Dynamical Systems in Neuroscience: The Geometry of Excitability And Bursting[END_REF][START_REF] Ermentrout | Foundations of Mathematical Neuroscience[END_REF]. Here we describe its mechanism in our model (for more detail, see [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF]). Consider a cell, at rest (non oscillating), thus either in region A, B of the bifurcation map or in the rest state (Si) of region D. This cell can be lead to an excited state with fast oscillations (stable limit cycle in region C or D) from different scenarios discussed in the next section. Here, we are just interested in what happens when the cell starts to produce these fast oscillations.

In this situation, the average voltage is high (∼ -30 mV, see in Fig. 1, the phase portraits of region C, D). The rising in voltage increases the intracellular calcium concentration and triggers the generation of an increasing sAHP current, according to eq. ( 5). In the bifurcation diagram this corresponds to an increase in G S , which eventually leads the cell to region A or B where no stable limit cycle exists. Thus, fast oscillations stop and the cell is feeling a strong, slow, after hyperpolarization current that slowly decreases its voltage. As a consequence, its intracellular calcium production decreases, thus sAHP, leading the cell slowly back to its rest value. This process takes around one minute during which the cell is hyper-polarized. This is how bursting is produced (see Fig. 3 in [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF]).

Note that this corresponds to a pseudo-cycle, shown in Fig. 2, in the subspace of variables C, R, S, controlling the sAHP production. The following remark is important. The values taken by R during the cycle, R ∈ [0.3, 0.7], correspond, for g S = 10 nS, to G S ∈ [0.081, 2]. As already mentioned, R is quite slow compared to Ach. So, during the re-polarization phase, G A has a very low value close to 0 (see e.g. Fig. 6). Thus, when getting back to rest with a low G A value the SAC successively crosses region B, A then D. When a cell is in region B or A the only possibility to excite it is to shift it to region C, requiring that G A is large enough, namely, that neighbours cells are bursting and g A , the maximal Ach conductance is very large. Otherwise, it is not possible to excite it. When a cell is in region D it can be excited, either by increasing G A and shifting it to region C, or, by jumping from the rest state (Si) to the period orbit, e.g. by a small amount of noise as developed in section 3.1.1. All in all, this implies that, on its way back to rest after hyperpolarization, a cell undergoes two sub-phases. An absolute refractory (AR) phase where it is not possible to excite it and a relative refractory phase (RR) where it can be excited again. In particular, in region D, cells are in the RR phase. Actually, the cell spends a relatively long time in the relative refractory phase. Although it can be excited anywhere in this phase, the level of excitation required to make it enter in the fast oscillations strongly depends on its location in region D (see Fig. 4), i.e. on the actual value of G S and G A . "FP" means "Fixed Point", "Repol" means "Repolarization". Here, we have plotted one piece of trajectory corresponding to a single bursting period. In general, successive burst correspond to non superimposing trajectories, due to the effect of coupled cells.

From single cell dynamics to waves

From the previous analysis we have learned that a SAC can either be: (1) at rest (region D); (2) excited with fast oscillations (region C or D); (3) hyperpolarized, either in the absolute refractory period where it is not possible to excite it (region A or B), or, in the relative refractory phase (region D, stable fixed point). When excited, each SAC undergoes the excitation-hyperpolarization cycle (Fig. 2) before eventually returning to rest.

We now analyse how this individual dynamics, coupled via Ach, induces retinal waves and controls waves initiation, propagation and stop.

The mechanism of wave initiation

Let us consider a situation where all cells are at rest so that they cannot burst without an external influence. By assumption, no cell is excited and no external influence can come from other cells. Actually, as the dynamics (1)-( 9) is deterministic, a SAC at rest would therefore stay at rest indefinitely, without a possibility to trigger a wave. However, there is always noise in real neurons due e.g. to fluctuations in ionic currents from the random opening of ionic channels.

We model this by a white noise ξ(t), whose amplitude is constant with time and controlled by a parameter η. In our formulation η ξ(t) has the physical dimension of a current, hence the white noise ξ(t) has a physical dimension of ms -1 2 , η is in pA ms 1 2 . The influence of this noise is studied in detail in the paper [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF] for single and isolated neurons. Here, we outline the main consequences on retinal waves. Another scenario of wave generation, not requiring noise, is discussed in the discussion section 6.1.

Noise Induced Bursting

In region D, neurons are bistable: a neuron in the rest state (Si) can go to the excited oscillatory state (sPO) provided it receives a kick (here, noise) allowing its state, in the fast-dynamics phase space V -N , to jump above the Sd barrier and to be trapped in the attraction basin of the stable periodic orbit (Fig. 1).

Crossing the saddle barrier leads to the stable periodic orbit where the neuron stays a certain time, producing fast oscillations, returning eventually back to the Si attraction basin. So, in region D we observe periods of fast oscillations punctuated by return to rest. Note that, in region A, it is possible to observe isolated spikes, with fast return to rest, because noise can make the trajectory going beyond the saddle barrier. As there is no attractor beyond this barrier, the trajectory eventually returns back to the attraction basin of Si, following the unstable manifold of Sd, which is attracted to Si. This takes a bit of time though. In the trajectory, this effect is manifested by one or a few short spikes in the voltage. We want to differentiate these two situations focusing on bursting cells. Indeed, to anticipate on waves propagation, we need SACs to be in the excited state long enough to produce Ach emission so as to eventually excite the neighbours. As the characteristic time of Ach growth is of order 1 -2s (a few µ -1 A = 0.53 s), the cell must be excited for a duration of this order. This does not happen in region A. In region D it does happen, although it depends on the noise intensity.

A good, numerical and automatic indicator to ensure that a SAC is excited long enough, is to fix a threshold θ on calcium concentration (here θ = 4 C (0) , where C (0) is the baseline calcium concentration). We consider that a burst arises at time t when C(t) > 4 C (0) . As we checked, this condition is a good criterion to ensure that a sAHP current is generated. Note that the characteristic time for calcium rising, τ C = 2 s (eq. ( 2)) is of the same order as the characteristic time for Ach production, so this criterion is also valid to ensure that SAC is excited long enough to produce Ach.

Bursting probability map

Here we want to compute the probability that calcium dynamics crosses the threshold θ ensuring the generation of Ach and sAHP for the full dynamics. As in the previous section, we consider G S , G A as constant.

The rest state corresponds to dV dt = 0 and dN dt = 0 (fast dynamics) in the absence of noise, and for fixed values of G S , G A . This corresponds to a voltage:

V Si = g L V L + g C M * V C + g K N * V K + G S V K + G A V A g L + g C + g K + G S + G A (11)
where

M * ≡ M ∞ [ V Si ], N * ≡ N ∞ [ V Si ]
. This is thus, in general, a non linear function of the parameters G S , G A . Now, due to noise, the voltage V has fluctuations δV around the rest state,

V = V Si +δ V . Assuming that δV is small enough so that M ∞ (V ) ∼ M ∞ (V Si ) + δV M ∞ (V Si ), N ∞ (V ) ∼ N ∞ (V Si ) + δV N ∞ (V Si )
, the fluctuations are given, to the first order in δV , by the stochastic differential equation

dδV dt = - 1 τ (G S ,G A ) δV + η Cm ξ t , where the characteristic time τ ≡ τ (G S , G A ) = Cm G * (G S ,G A ) with: G * (G S , G A ) = g L +g C M * + M * ( V Si -V C ) +g K N * + N * ( V Si -V K ) +G S +G A , (12)
where M * = M ∞ (V Si ) and the same for N * . Thus, δ V is, in this first order approximation, a Ornstein-Uhlenbeck process with mean zero, covari-

ance C(t, t ) = η 2 2 C 2 m τ e -t-t τ -e -t+t τ
, for t ≥ t , and variance σ 2 (t) =

η 2 2 Cm G * 1 -e -2 t τ . When t is large enough, t τ , this variance converges to a constant σ 2 ≡ σ 2 ( G S , G A ) with: σ(G S , G A ) = η 2 C m G * (G S , G A ) . (13)
Let us now call V Sd the voltage of the saddle point in region D. A necessary condition for bursting to occur is that, at some time t,

V Si + δV (t) > V Sd corresponding to a probability 1 -Π V Sd -V Si σ(t)
, where Π(x) =

x -∞ e -y 2 2 √ 2π dy. This probability is a sigmoidal function, where the inflexion point depends on G S , G A (via G * ), as well as the slope (via σ(G S , G A )) (see Fig. 4). More generally, the density and the cumulative distribution of the first passage time can be computed exactly [START_REF] Borodin | Handbook of Brownian motion: facts and formulae[END_REF][START_REF] Alili | Representations of the first hitting time density of an ornstein-uhlenbeck process[END_REF]. In the simplest case, it is given by an erf function, and more generally, by Hermitte polynomials. However, this is not enough for what we want here. We want to compute the probability, ν(G S , G A ), that, for given values of G S , G A the calcium concentration crosses the threshold θ. In terms of V , this means we want to estimate the probability that, starting from V = V Si at time 0, V (t) crosses the Sd barrier and stays long enough beyond this barrier to let calcium rise and cross the threshold. We have not been able to compute analytically this probability. This could may be analytically done using Girsanov-Freidlin-Wentsel [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF] estimates (considering the stable limit cycle, an ω-limit set, as an equipotential set), but this is beyond our technical reach.

We made therefore numerical simulations to compute this probability in the plane {G S , G A }. We proceeded as follow. For each value of G S , G A on a grid with resolution 10 -2 , one generates N samp = 100 trajectories where, for each trajectory, one starts from the rest state, iterates the dynamics with noise over a maximal time duration T = 20 s, and seeks the occurrence of an increase in the calcium conductance above the threshold θ. The estimated probability ν(G S , G A ) is the number of occurrence of such long bursts divided by N samp .

The result is shown in Fig. 4. As expected, this is a sigmoid function, well approximated by a form:

ν(G S , G A ) ∼ f G A -G Ac (G S ) σ(G S , G A ) (14)
(Fig. 4, top) where f is a sigmoid. Here, G Ac is a threshold value of G A , depending on G S and constrained by the SN 2 and Homoclinic bifurcation line (eq. ( 10)). Likewise, σ(G S ), the inverse slope of the sigmoid function, is constrained by the Homoclinic and SN 2 bifurcation lines. σ(G S ) depends on η, the noise intensity and tends to zero as η → 0. Also, it decays to zero like the distance between the SN 2 and homoclinic bifurcation line. Thus, beyond the point where SN 2 and H c meet (at the SNH point), the transition is zero-one, corresponding actually to direct transition from region A to region C by a SNIC.

Let us make a few remarks to conclude this section. First, on the role of noise intensity. When η → 0 (deterministic case), ν → 0: there is no possibility to burst in the absence of noise. For ν > 0 the probability is a sigmoid function with a transition zone delimited by the boundaries of region D. When noise level is low, it takes a long time to the cell to jump to the stable periodic orbit, but when it is there, it stays sufficiently long to produce Ach dynamics, and, for calcium dynamics to reaches the threshold θ. When η increases, the time to jump on the sPO decreases, but the time where the trajectory stays on it decreases too. If η becomes too large the cell cannot stay long enough on the stable periodic orbit to produce enough Ach to excite its neighbours. This implies that there is an interval of η values where retinal waves can actually propagate.

Second, the probability to trigger a burst is highly sensitive to the value of G S , G A because ν is quite sharp. Consider now a SAC on its way back to rest in the cycle of Fig. 2, somewhere in region D, thus in the relative refractory phase.

Its level of G S depends on its position in the cycle, its level of G A depends on the state of its neighbours. And a tiny change in these value dramatically change the probability that this cell becomes excited. Especially, taking again into account the time scale separation between hyperpolarization and cholinergic coupling, and considering the G A dynamics for a constant G S , a tiny change in the level of excitation of its neighbours dramatically impacts the probability the SAC has to become excited. This has a strong impact on waves propagation and waves distribution.

The mechanism of waves propagation and stop

When a cell is in the excited, oscillating, state, it produces Ach, according to [START_REF] Bántay | Self-organization and anomalous diffusion[END_REF], [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF], on a time scale of order 2 seconds, increasing the Ach current received by its post-synaptic neighbors. This increases the G A conductance of its neighbours cells. We investigate what happens to a cell receiving an Ach current from its neighbours. There are actually 3 cases, depending where the cell is located in the bifurcation map.

1. The cell is already excited (region C or sPO in D). Then, it produces Ach thereby exciting the cells which are exciting it, inducing a strong non linear feedback coupling mechanism, enhancing the Ach production, until the sAHP production takes place. This is illustrated, for two cells, in Fig. 5. This mutual coupling increases the Ach conductance G A of coupled cells, corresponding to moving vertically in the bifurcation map of Fig. 1. So, a cell initially on the sPO of region D (e.g. excited by noise) is eventually lead to region C where its oscillations are now essentially insensitive to noise.

2. The cell is hyperpolarized, in the absolute refractory regime, (region B or A too far from region C to be lead in C by an increase of G A ). Here, the cell cannot be excited no matter how it is perturbed. SACs in this state constitute therefore an impassable barrier where no wave can go through.

3.

The cell is hyperpolarized, in the relative refractory regime.

Here, the cell can actually be excited thanks to different modalities, which can combine. First, because the cell feels an increasing G A it moves vertically in the bifurcation map and can be eventually lead to region C.

Second, it can be excited by noise, when it is in region D, either in the rest state or on its way back to rest in the excitation-hyperpolarisation cycle of Fig. 2. This last case is the most interesting. Indeed, as discussed in the previous section, when a cell is moving in region D the probability The blue dashed line correspond to a value of g A , computed in section 5.2.1, eq. ( 25), below which cell 1 can not trigger a burst in cell 2. These two figures illustrate the mutual non linear interaction of two bursting cells. Especially, the Ach current generated by 2 prolongs the burst of 1. Compare, indeed, the black trace for cell 1 (no coupling so that cell 2 stays at A = 0), to e.g. the blue trace of cell 1 (g A = 0.05) where the burst is longer. Both cells are then mutually coupled and this mutual coupling reinforces their synchrony. As shown on the right, the effect of this coupling on burst duration τ 1,max depends on the regime, I,II,III.

to become excited due to noise depends dramatically on its position in From this simple description we anticipate that there is a minimal value g A0 ≡ g A0 ( G S , η ) for g A , the maximum Ach conductance, below which a wave cannot propagate. This value is computed in section 5.

We have here the core mechanism of waves propagation and stop. When a SAC starts to burst it excites its post-synaptic neighbours. This corresponds, for those excited neighbours, to an increase in G A , which either increases their probability to be excited by noise in region D or, if g A is large enough, leads them to region C. This way, one can generate a chain reaction, a retinal wave. Now, due to the sAHP activation, excited SACs eventually leave the bursting phase and go to an hyperpolarized state where they cannot be excited any more during a few tens of seconds. Thus, a wave leaves, behind it, an absolute refractory region where the forthcoming waves cannot go through for a long period. This induces spatio-temporal waves interactions where each wave has to propagate into a transient landscape of excitable/absolute refractory/ relative refractory cells and can collide. Before addressing the detail of waves propagation let us finish this section with two important remarks.

Link to forest fires

This dynamics bares some analogy with forest fires [START_REF] Bak | A forest-fire model and some thoughts on turbulence[END_REF][START_REF] Chen | A deterministic critical forest fire model[END_REF][START_REF] Drossel | Self-organized critical forest-fire model[END_REF][START_REF] Grassberger | Critical behaviour of the drossel-schwabl forest fire model[END_REF], as discussed in [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF][START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF]. However, there is fundamental difference, illustrated in Fig. 5: the mutual non linear coupling between cells prolongs the bursts, a feature which is definitely not present in forest fire models. Actually, this prolongation effect has been reported in Zheng et al experiments [START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF]. It is possible to block Ach receptors with a chemical cocktail. In the presence of this cocktail the burst duration is reduced. Another consequence of this non linear coupling, observed as well in experiments, is to increase the sAHP duration, hence, the persistence of refractive zones, with strong consequences on waves propagation.

Therefore, due this non linear interaction, retinal waves do not reduce to "classical" paradigm of forest fire as proposed e.g. in [START_REF] Bak | A forest-fire model and some thoughts on turbulence[END_REF][START_REF] Drossel | Self-organized critical forest-fire model[END_REF] or in [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF][START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF] for retinal waves.

However, the so-called "fire-diffuse-fire" model, appearing for intracellular calcium waves and introduced in [START_REF] Dawson | Fire-diffuse-fire model of dynamics of intracellular calcium waves[END_REF] actually present some important features of what we observe in our model: bi-stability, continuous versus saltatory propagation depending on a threshold for a characteristic parameter. The mathematical study presented in [START_REF] Thul | A bidomain threshold model of propagating calcium waves[END_REF], exhibits conditions to have this "fire-diffuse-fire" in a spatial landscape with bidomains and, in particular, the wave speed is computed. More generally, intracellular calcium waves share many properties with retinal waves: the necessity of an initial nucleus to seed the wave, the spread of excitability and the emergence of refractory domains, the effect of stochasticity (which renders calcium waves functionally robust and adaptive to changing environmental conditions [START_REF] Thurley | Fundamental properties of ca2+ signals[END_REF]). In contrast to these works, focusing on Ca dynamics modelled by phenomenological equations of reaction-diffusion type, in our work

Ca dynamics is part of a process involving bursting at fast time scale and Ach coupling at medium scale: cholinergic waves generate calcium waves which are related to the existence of refractory sAHP domains. Therefore, a possible modelling would be to couple our Ach dynamics to calcium waves shaping the sAHP landscape. This opens up a way to interesting studies at least on the numerical side, the mathematical analysis being seemingly (at least to us) out of reach.

The effect of varying Ach coupling during development

In this section we analyse numerically waves dynamics as g A , the maximal Ach conductance, varies. The Ach conductance of a nicotinic channel is proportional to A 2 γ A +A 2 with γ A ∼ 1 nM. Even if the Ach concentration resulting from the puff has been reduced when the maximum of the induced current is reached we may assume that the remaining Ach concentration is quite larger than K d . This means that A 2 γ A +A 2 ∼ 1 (saturation). Note that receptor desensitisation is quite important during a puff, as the receptors are completely overwhelmed with transmitter, so this assumption seems reasonable. Therefore, in our model the current resulting from the puff is This gives the following table of values for g A , where we added, in the comments line, the corresponding stage of development (rabbit). Stage III corresponds to glutamatergic waves birefly commented in section 2.1. What is important here is the huge decay of g A during development (2 orders of magnitude).

I puf f = -ng A (V clamp -V A )
Beyond the important observation of a rapid decay of the cholinergic synapses during maturation, this table provides us hints to explore the central topic of this paper: how could stage II retinal waves dynamics be impacted by a continuous change in the cholinergic coupling ? In the next sections we address in more detail the mechanism of wave propagation and interaction in a random, slowly evolving, landscape of G S , G A . We also study numerically the wave sizes Date (days after birth) [START_REF] Maccione | Following the ontogeny of retinal waves: pan-retinal recordings of population dynamics in the neonatal mouse[END_REF], (see Fig. 8) .

I puf f (pA) g A (

The setting of numerical simulations

Let us start with a few operational definitions used for numerical investiga-655 tions. We consider a wave as a causal event. The bursting of a cell at time t is either spontaneous, or caused by bursting cells connected to it. As a consequence, a wave is mathematically a list of events of type (t start ,S t ). S t is a subset of bursting cells at time t defined by the following constraints:

(1) The wave starts from a cell, or a group of connected cells, S start starting to burst at time t start .

(2) For t > t start , S t must intersect S t-dt ∪ ∂S t-dt , where ∂S t-dt is the neighbourhood of S t-dt . This mathematically corresponds to causality and connectedness.

(3) The wave stops at the last time t where S t is non empty. This time is called

t end .
Therefore, two disjoint clusters S t , S t (i.e. with no common cell even at the boundary) correspond to two distinct waves. When two waves intersect we still consider them as distinct, they do not penetrate each other. Thus, from this perspective, they do not spread any more at the intersection points.

We call "size" the total number of cells involved in a wave, i.e. s = t end t=tstart #S t , where #S t is the cardinal of the set S t . The sum is discrete because time is discrete in simulations. The duration of the wave is D = t end -t start . The global activity, noted n, is the total number of active cells at a given time.

All our simulations are limited to a one dimensional space. This is due to several reasons. First, considering the huge difference between the fastest and slowest time scales in the model, simulations in 2D, for reasonable sizes (∼ 100 × 100 cells) are enormously time consuming. In addition, the definition of waves boundaries is easy in 1D (boundaries are points), while these are irregular frontiers in 2D, requiring sophisticated algorithms. As we show in the next section, it is actually important to separate waves when studying their statistics, instead of considering the total number of active cells per unit time, independently of their localisation. Finally, these simulations are intended to illustrate the main content of this paper, which is theoretical. We discuss about the effect of dimensionality and connectivity in section 5.3.3.

29

In the simulations we considered characteristic times τ S and τ R , which fix the time scale of sAHP, of order 10 s, while, in experiments, they are rather of order a minute. We did that to reduce the computational time, especially, the refractory times between successive waves.

Waves dynamics and scaling 690

We now illustrate the effect of increasing g A on waves dynamics in our model.

We essentially observe 3 regimes, already commented above and now more precisely defined and illustrated by Fig. 7,8,9. Note that we used the labelling I,II,III for these regimes, not to be confused with the retinal waves stages I,II,III.

Our analysis sticks at stage II retinal waves. Regime II. Waves are competing through the sAHP landscape. When g A 700 further increases, a bursting SAC increases enough the probability of this neighbours to be excited so that some neighbours become eventually excited. The mutual excitation drives the cells in a robust synchronisation regime, as illustrated in Fig. 5. Waves can propagate but they are blocked by the sAHP landscape left by previous waves. Here, therefore, waves size 705 is dramatically dependent on the waves history. In this regime waves characteristics (size, duration, activity) are quite sensitive to g A and their probability laws have a large tail.

Regime III. Waves are cooperating through Ach coupling. When g A is large enough, waves may eventually spread through the whole lattice. In this 710 regime, waves born at different times, from different sites, eventually collide with the effect of prolonging the burst durations. When g A increases further, one is progressively driven to a full synchrony regime where all cells fire synchronously with a period fixed by sAHP characteristic times (τ S , τ R ). These 3 regimes are first roughly indicated in Fig. 7, where we plotted, in Fig. 7(a) (left), the probability ρ that a cell is bursting at a given time, as a function of g A for N = 20, 50, 100, 200, 400. We essentially see abrupt changes in the slope of ρ, allowing to (broadly, for the moment) locate these 3 regimes (Inset of Fig. 7 (a)). The separation between regime I and II holds at the unique point where ∂ρ ∂g A is maximum, that is the point where the probability ρ that a cell bursts at a given time has the maximum sensitivity to variations of g A .

Regime III is the range of g A values where ρ increases linearly. We have also indicated, in red, 4 points corresponding to the three regimes (with one point at the transition between regime I and II), used to characterize in more details waves behaviour. In Fig. 7, center, we have plotted the average wave size (label (b)) and duration (label (c)) as a function of g A , for the same values of N .

Recall (section 4.2) that we call "size" the total number of cells involved in a wave, where a wave is composed of connected bursting cells. The global activity (total number of cells bursting at a given time) is plotted on the right of Fig. 7.

Similar regimes have been reported in the literature, especially in [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF][START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF].

What these studies neglected though is the role of closeness to bifurcations, non linear feedback, and their consequences on waves propagation.

Waves propagation in the different regimes

This propagation is illustrated in Fig. 8. From left to right each column corresponds respectively to regime I (g A = 0.04 nS); regime II (g A = 0.08 nS corresponding to the point II 1 in Fig. 7 (a) and g A = 0.1 nS, point II 2 in Fig. 7 computed as:

C (i) = 1 T T 0 C i (t) C ic (t) dt - 1 T T 0 C i (t) dt 1 T T 0 C ic (t) dt, ( 15
)
where T is the total duration of the simulation, i c is the label of the center cell, and i the label of a cell in the lattice. In Fig. 8, bottom row, we have displayed the normalized correlation C (i) C (ic) so that the maximum is 1.

In regime I, (left column) the entrance of a SAC in the bursting regime initiates a propagation which ends very soon. This is because, for this regime of parameters, the Ach current, even at its maximum, is not sufficient to cross the S d barrier. Only noise affords it, with a low probability, depending on noise intensity and on g A . As a consequence, ρ, s , D , n grow slowly with g A . In regime I the calcium correlation decays rapidly with the distance to the center cell: cells are essentially decorrelated.

On the opposite, in regime III, (right column), top, we see waves spreading through the whole lattice. Waves appear almost periodically with a frequency ∼ 35 s, essentially constrained by sAHP dynamics (i.e. the return time to region D dependent on the times τ S , τ R ) and the time necessary to activate a cell and, thereby, restart a wave. In this regime ρ grows slowly, and s , D , n saturate. The saturation of s , D is due to our modelling choice: when two waves intersect we consider them as distinct. Thus, they do not spread any more at the intersection points. As a consequence, these quantities do not depend on the lattice size. The plateau is actually determined by the probability to start a wave and the wave speed. In this regime, n, the total activity, behaves in a different way. One expects the total number of active cells to increase with the lattice size, but, as we see in figure 8, all cells in the lattice are not necessarily simultaneously active. As g A further grows, however, one achieves a full synchronisation where all cells are active. These full synchrony regime is periodic and alternate of full lattice refractoriness. Finally, the instantaneous calcium correlation is strictly positive throughout the lattice, confirming that, in this regime, the correlation length is the lattice size.

The intermediate regime II is shown in Fig. 8, middle columns. On the top row, we see how waves propagation has to cope with the sAHP landscape left by the previous waves. Spreading waves are stopped by the sAHP trace left by previous waves. Therefore, waves have a widespread distribution (see also Fig. 9).

Calcium correlation also shows an interesting profile. It reveals that correlation alternate from positive to negative as the distance to the center cell increases.

Positive correlations correspond to cells in the same state as the central cell i c (i.e. bursting when i c is bursting; hyperpolarized when i c is hyperpolarized), whereas negative correlations correspond to cells in the opposite states. Thus, calcium correlations is another way to reveal the existence of refractory/active domains shaping waves propagation and interactions. Increasing g A from regime I to regime III via regime II leads to an increase of the correlation length, which has the size of the lattice in region III.

In these figures one clearly see what has been announced in the beginning of this section. In regime II waves interact via their sAHP landscape, thereby obstructing each other. In contrast, in regime III they have a direct interaction and the mutual Ach feedback prolongs the cells activity. The transition between these two regimes is reminiscent of a percolation transition, as anticipated in [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF]. However, the transition we observe holds at a specific g A value and is therefore not robust to changes of g A during development. This aspect is further discussed in section 5.3.3.

Waves statistics

A synthesis of the observed waves statistics is shown in Fig. 9, where we plot the probability distribution of s (top row), D (middle row), n (bottom row),

for N = 50, 100, 200, 400, in the 3 regimes. As before, each column corresponds to a regime (I, II 1 , II 2 , III from left to right) The plots are in log-log scale to investigate potential power laws. We observe first that the distribution of s and D do not depend on N , as expected from waves interactions. In contrast, P (n) does. One also remarks that the shape of these probabilities evolves with g A . Let us interpret what we observe thanks to the previous dynamical analysis and especially Fig. 6.

In regime I, the bursting of cells is only induced by noise. When bursting, a cell produces Ach and thereby increases the probability that a neighbour cell bursts, due to noise. As a consequence, there is an important probability to observe waves of size 1 and quite a smaller probability to observe waves of larger size. The probability distribution of waves duration D is more interesting. After a flat maximum, for D ≤ 1.2 s, this probability decreases to a local minimum, at D ∼ 1.8 s before increasing to a secondary maximum, at D ∼ 2 s. The first maximum corresponds to waves of size 1. Their duration probability decreases as the duration increases. However, waves of size 1 with a sufficiently large duration produce more Ach (in a time scale of the order 1-2 s, a few µ -1 A = 0.53 s) thereby increasing the probability to have a neighbouring cell that starts to burst, thanks to noise. But, from Fig. 5 we know this neighbouring bursting cell will increase the burst duration of the initial cell. This explains the second peak. In other words, the local minimum in the distribution marks a threshold above which the mutual coupling between two bursting neighbours prolongs the (size 2) waves duration. Finally, the probability distribution of the global activity has a bell shape with a maximum ∝ N s . Indeed, the probability of this variable is constrained, on one hand, by the number of coexisting waves and on the other hand by the joint probability that each of these waves has a certain size and that these sizes sum up to n. This explains the bimodal shape of the distribution 2 . Actually, this bimodal shape is observed in each regime 3 .

In regime II, as illustrated in Fig. 6, there is a small probability for a cell to penetrate in region C where it bursts without requiring a noise kick. Waves propagate until they rich a sAHP barrier, limiting their size and generating a long tail in the distribution. There are still waves of size 1 though, whose probability decreases with g A . For the duration we still see the local minimum and second maximum, reminiscent of regime I, but beyond the second maximum waves duration probability has, like sizes probability, a longer and longer tail as g A increases. The global activity probability still has a bell shape, but it looks as a mix of two regimes where it growths linearly (in log-log scale) up to the maximum and decreases also linearly. The maximum appears therefore as a crossover point between two power laws regimes.

In regime III, the probability to have waves of size 1 is very small. A bursting cell triggers a wave which propagates through the lattice until it meets another wave. The distribution of size is bell-shaped, with a maximum corresponding to our definition of waves size: although, in regime III, waves cooperate, when two waves meet their size does not increase. The durations distribution has a similar bell shape, starting, for g A = 2, around 2 s with a maximum around 3 s (the average duration being 5.6 s). 2 The probability of n is related to the probability of sizes by P (n) = N m=1 s 1 +s 2 +•••+sm=n P (s 1 , . . . , sm) where the sum s 1 +s 2 +•••+sm=n means the sum over all possible configurations with m sets of connected active cells (waves), each of size

s i , i = 1 . . . m with s 1 + s 2 + • • • + sm = n,
and with the constraint that the total number of cells is N . As, in regime I, waves are small and not interacting, P (n) = N m=1 s 1 +s 2 +•••+sm=n P (s 1 ) . . . P (sm).Assuming, for the simplicity of the argumentation, that P (s) ∼ Kγ -s we obtain

P (n) = γ -n N m=1 K m S N (m, n)
where S N (m, n) is the number of possibilities to make, over the set {1 . . . N }, m connected group with s i elements and with m i=1 S i = n. We therefore see that P (n) is a combination of a term decaying with n multiplied by a combinatorial term S N (m, n) which first increases with n, then decays. 3 In the other regimes the argument developed in the footnote 2, extends as we essentially need that the joint probability P (s 1 , . . . , sm) decays sufficiently fast (power low or exponential)

with n.

There is no clear evidence that the probabilities of connected waves observables (s and D) have a power law tail, although this tail is large, especially for duration which goes from 1 to 20 s. Actually, probability distributions might depend on the dimension and the connectivity, as developed in section 5.3.3.

Transport

Finally, we analysed the wave speed as a function of g A . First, we remarked that, due to waves interactions and sAHP lanscape, waves are not necessarily propagating ballistically, that is r(t) = c.t z , where r(t) is the wave radius, z the anomalous diffusion exponent and c the wave "speed" (which actually is the classical speed for z = 1, whereas, it is e.g. the diffusion coefficient if z = 1 2 ). Anomalous diffusion is further analysed in section 5.

We show a plot of c, z as a function of g A in Fig. 10. Note that the estimation of these quanties is difficult in regime I as waves are very small. We observe that z is (on average, with large error bars) smaller than 1 in regime I and III corresponding to anomalous transport. This is due to different effects though.

In regime I, waves are small and propagation is due to the effect of noise for cells essentially located in region D (see next section and Fig. 6). Thus, wave propagation is slower than ballistic. In contrast, in regime III, waves interact together. When a wave intersects another wave its radius does not evolve any more whereas its duration is incremented. Thus, the effective exponent z is lower than 1. Although it might look as an artefact of our wave definition, this criterion actually clearly shows when direct wave interactions becomes significant. To illustrate this point we have added in the figure the lines separating the 3 regimes where we clearly see how propagation is slowed down in regime I and III.

We also plot, at the bottom of Fig. 10, the speed c, together with the theoretical computation, eq. ( 26) developed in section 5. Considering that SACs are spaced by a 50 µm our speed would correspond to real waves speed in the range [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF]200 ] µm/s in agreement with experimentally observed ranges [START_REF] Singer | Potentiation of L-Type Calcium Channels Reveals Nonsynaptic Mechanisms that Correlate Spontaneous Activity in the Developing Mammalian Retina[END_REF][START_REF] Maccione | Following the ontogeny of retinal waves: pan-retinal recordings of population dynamics in the neonatal mouse[END_REF][START_REF] Sernagor | Differential effects of acetylcholine and glutamate blockade on the spa-tiotemporal dynamics of retinal waves[END_REF]. Here, transport is ballistic too but the radius of a wave stops to grow, in our definition, when two waves collide, whereas duration is incremented. This explains why z < 1. Bottom. "Speed" c i.e.

proportionality coefficient in the relation x(t) = c t z (see text). The orange line corresponds to the theoretical derivation [START_REF] Dawson | Fire-diffuse-fire model of dynamics of intracellular calcium waves[END_REF]. The deviation from the theoretical curve in regime I and III is explained in the text.

Transport equations

Up to now, we have shown how the non linear coupling between bursting cells and the recurrent return in region D give rise to a rather rich dynamics with different regimes of waves. Still, the numerical illustrations have been done 875 in one spatial dimension, with nearest neighbours interactions. So, one may ask whether a larger dimensionality or a different connectivity could notably impact the global dynamical picture. To avoid multiplying tedious numerical simulations with a larger range of possibilities when considering connectivity patterns 4 , an alternative idea is to derive canonical transport equations, holding in 1 and 2 dimensions.

In this section, we derive such a transport equation characterising the retinal waves propagation in a landscape of sAHP. These equations are constructed in the spirit of the bifurcation diagram exposed in section 2.3, with a main idea: waves transport Ach and sAHP conductances, i.e. a wave is mechanism where G S and G A , considered as fields, propagate, obeying a transport equation that we derive. This is done taking into account the time scales separation between G S and G A . From this approach we are able to compute a few important characteristics of waves dynamics, like speed in the ballistic regime. We also compare our transport equation to known equations in the non linear and statistical physics literature.

Ach conductance transport

Approximations

The equation for Ach production of neuron i is given by ( 8). We introduce the variable:

Γ i = k∈Bi U (A k ) ( 16
)
with

U (A) = A 2 γ A +A 2 . Then, the Ach conductance of neuron i at time t is G A (i, t) = g A Γ i (t).
Differentiating Γ i with respect to time gives:

dΓ i dt = k∈Bi U (A k ) [ -µ A A k + β A T A (V k ) ] .
The goal here is to obtain a closed-form equations for the time evolution of Γ i , thus for G A (i, t). For this we are going to make several approximation based on the preliminary analysis of dynamics.

(i) Piecewise linear approximation for U . We remark that U (A) is a sigmoid that can be approximated by the piecewise linear function (Fig.

12, left, top):

U (A) ∼          A 2 √ γ A , 0 ≤ A ≤ 2 √ γ A ; 1, A ≥ 2 √ γ A . (17)
In this approximation a cell k such that A k > 2 √ γ A does not contribute to the variation of Γ i because U ∼ 0. We observe that in our simulations

A < 2 √ γ A = 2.
Note that the form of T is inspired from biophysics but its shape (sigmoid) matters more than its detailed mathematical form.

Under the approximation ( 17)

U (A) ∼ U (A)A so that dΓi dt = -µ A Γ i + β A 2 √ γ A k∈Bi T A (V k ).
(ii) Mean-field approximation for T A . The time scale of Ach evolution (∼ 2s) is quite fast compared to the time scale of evolution of V (a few ms). As we are interested in transport when G S , G A vary it is relevant to replace T A (V ) by its time average on the fast time scale,

T A (V ) = lim T →+∞ 1 T T 0 T A (V (s)
) ds where G S , G A are kept constant. We then make the assumption that the time average of Ach production can be approximated as T A (V ) ≡ Ω(G S , G A ) ∈ [ 0, 1 ]. On the time scales of G S , G A Ω(G S , G A ) evolves along the pathway described by G S , G A in the bifurcation map (Fig. 6). The function Ω(G S , G A ) is shown in Fig. 12 right top. It has a sigmoidal shape, with a maximal slope in region D.

In this figure, bottom, we also show a plot of the gradient ∂Ω ∂G S , ∂Ω ∂G A and Hessian of Ω.

(iii) Continuous space limit. Cells are located at the nodes of a Z d lattice with small lattice spacing a ∼ 50 µm, d = 1, 2 (see section 2.1). We now consider G S , G A as continuous fields in space and time, G A ≡ G A (x, y, t).

Then, using the smallness of a (compared to the waves characteristic size):

∂Γ ∂t = -µ A Γ + d β A √ γ A Ω(G S , G A ) + D A ∆ Ω(G S , G A ).
where we have set:

D A = a 2 β A 2 √ γ A , (18)
the diffusion coefficient (∼ 62.500 µm 2 s -1 ), ∆ being the Laplacian operator. The first term of the right hand side comes from Ach degradation, the second comes from the local Ach production, and the third corresponds to Ach transport by cells interactions, as we develop now. From now on we will write Ω instead of Ω(G S , G A ) to alleviate notations whenever it makes no confusion.

Non linear diffusion

Using G A = g A Γ, for obtain a transport equation for Ach conductance:

∂G A ∂t = g A D A ∆ Ω + S. (19)
with:

S = -µ A G A + g A d β A √ γ A Ω, (20)
The transport equation is a non linear diffusion equation, because Ω depends non linearly on the fields G S , G A via eq. ( 14). Note that this is not a reactiondiffusion equation, in contrast e.g. to Lansdell et al approach [START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF]. Indeed, in their case Ach freely diffuse in the medium and non linearity is taken into account in a reaction term whereas, in our case, non linearity is in the diffusion coefficient, characterizing how waves non linearly interact with the sAHP landscape. The equation we obtain is also different from transport equation in calcium fire diffuse models [START_REF] Dawson | Fire-diffuse-fire model of dynamics of intracellular calcium waves[END_REF][START_REF] Thul | A bidomain threshold model of propagating calcium waves[END_REF] although propagation and underlying mechanisms are very similar (see section 3.3).

The equation ( 19) contains also a "source" containing both the degradation term and a local production term. It can therefore be either positive or negative and depends on space-time via G S , G A . This source term is shown in Fig. 11.

Although it has a very different shape in the different regimes, it is slightly negative on average so that dissipation dominates. The source term is purely deterministic here, although it should contain as well a stochastic term that mimics noise induced bursting, at the time scale of Ach evolution. This is briefly discussed in section 5. the local production coming from SACs neighbours. In this regime dynamics is ruled by fluctuations around the threshold, not by the leading source term.

Let us comment more on the effects of this non linear diffusion. It is actually useful to let the gradients of G A and G S appear in eq. [START_REF] Butts | Retinal waves are governed by collective network properties[END_REF]. We also introduce the Hessian of Ω:

H Ω =   ∂ 2 Ω ∂G 2 S ∂ 2 Ω ∂G S ∂G A ∂ 2 Ω ∂G A ∂G S ∂ 2 Ω ∂G 2 A   . (21)
This symmetric matrix characterizes the local convexity of Ω. It also defines a metric. Indeed the transport equation( 19) can now be written:

∂G A ∂t = g A D A   ( ∇G S , ∇G A ) H Ω   ∇G S ∇G A   + ∂Ω ∂G S ∆ G S + ∂Ω ∂G A ∆ G A   +S. ( 22
)
This a Laplace-Beltrami diffusion on a manifold whose metric is constrained by [START_REF] Chalupa | The visual neurosciences[END_REF]. In other words, the diagram shown in Fig. 6 has now to be thought on top of the Ω manifold, Fig. 12 top, right, whose gradient and curvature, shown in figure 12 middle and bottom row, constrain the diffusion. We see that the transport of G A depends on the spatial profile of sAHP. In particular, the extra-diagonal terms of H Ω couples ∇G S and ∇G A , the gradients of sAHP conductance and Ach conductance. As expected, the first and second partial derivatives are non vanishing essentially near the bifurcation lines delimiting region D. Interestingly, ∂Ω ∂G S < 0 in this region corresponding to an anti diffusion (focalisation) mechanism whereas ∂Ω ∂G A > 0. In view of the shape of Ω and, more specifically, on its derivatives, one sees that there are 3 main phases for cells, corresponding to specific regions in the plane {G S , G A }. At a given time, the cells are in one of these 3 phases, with spatial continuity due to the diffusion term.

(a) Quiescent phase. In regions where Ω is small and its derivative vanish [START_REF] Butts | Retinal waves are governed by collective network properties[END_REF]

reduces to ∂G A ∂t = -µ A G A + g A d β A √ γ A Ω, so that G A ∼ g A d β A µ A √ γ A Ω is
small. These cells do not participate to propagation.

(b) Bursting phase. For cells which are located in region C, Ω is large but its derivatives are small thus the Laplacian term in (19) vanishes too. There is stationary solution of the same form as above, but

G A = g A d β A √ γ A µ A Ω is
now large. These cells are the bursting cells participating to waves.

(c) Intermediate phase, for cells in the relative refractory period, in region D starting to burst or close to burst. Here, Ω depends sharply on G S , G A and its variations impact dramatically the wave propagation, according to eq. ( 19). These cells constitute waves boundaries (fronts).

The physical picture which emerges from eq. ( 22) and Fig. 12 is actually very close to the description made by Bantay and Janosi [START_REF] Bántay | Self-organization and anomalous diffusion[END_REF], in a beautiful paper about Self-Organized Criticality that we essentially rephrase, although the situation we have here is a bit more complex due to the dependence in sAHP. At a given time, one may distinguish three different spatial regions corresponding to cells in one of the 3 phases above. In the quiescent and bursting phase the diffusion is very small, while at the interface between these two phases, the diffusion grows steeply so that there is a large flux which enters the intermediate level of activity of quiescent cells until they become excitable. As a result, the boundary of the quiescent phases moves towards the bursting region, until this later is completely absorbed by the intermediate phase. The behaviour of the other boundary is completely similar, the only difference is that the direction of the flux is opposite. Now, the motion of boundaries is, in our case, strongly constrained by sAHP.

sAHP conductance

Equation ( 22) illustrates how dynamics depends on the dynamics of the sAHP conductance G S , a very slow variable (of order a minute). Using slowfast analysis one can see that V, N, C and A evolve, at the time scale of R, S, along a slow manifold parametrized by R, S. In particular, at this time scale, the calcium concentration is a function C(R, S) of R, S, whose shape can be guessed from Fig. 2, with 4 phases. It results that the coupled dynamics of S, R at the slow time scale, t s , has the form dS dts = α S C(R, S)(1 -S) -S; dR dts = α R S(1 -R) -R and is deterministic. This generates an invariant closed curve represented in Fig. 13, left, while, on the right, we see the evolution of G S as a function of t. This is, actually, a deterministic curve of the form G S (t) = u(t-t b ) where t b is the time when a cell starts to burst. This means that we may write G S (X, t) ≡ u ( t -t last (X, t) ) where t last (X, t) is the last bursting time anterior to t, for the cell located at X. Thus, t last (X, t) is a stochastic process, whose probability depends on G S , G A and therefore evolves along the waves propagation landscape. Fig. 13 is only illustrative though as the shape of G S (t) can be quite more complex depending on parameters such as g sAHP . This complex dependence on G S , which makes retinal waves dynamics extraordinary rich, lifts the general analysis of transport equations out of our technical reach and we are considering simplified situations in the next sections.

G S depends therefore on the network history via t B (x, t). This is the mathematical expression of what was anticipated earlier (section 3.2). When a retinal wave propagates, transporting Ach conductance, it leaves, behind it a refractory trace in which subsequent waves have to propagate. This induces complex interactions between waves, history depend, as depicted e.g. in Fig. 8.

Noise at the medium time scale

The transport equations are deterministic. They lack a "noise" term that mimics the random initialisation of waves. The detailed form is this "noise" is not completely clear though. It can be obtained by adding, in the equation of Ach production, a stochastic term δV to the voltage, as we did in section 3.1.2 and then integrate over time T A ( V + δV (t) ), e.g. by expanding T A in series, generating terms in powers of κ A δ V . One obtains a term proportional to

g A κ A σ β A
√ γ A which is the time integral of a Gaussian process passing through a highly non-linear function (T A ). The resulting process is therefore neither Gaussian nor Brownian. In addition it must be proportional to ν, the probability (14) that a cell burst. A possible spatio-temporal scaling for such a noise has been proposed in [START_REF] Volchenkov | Quantum field theory renormalization group approach to self-organized critical models: The case of random boundaries[END_REF]. Its main interest is to afford a renormalisation group analysis of the dynamics, which could be helpful for the characterization of retinal waves statistics in 1 and 2 D.

Consequences

Threshold value of g A

A first consequence of these approximations is that there exists a threshold value of g A , g A0 , such that for g A < g A0 no wave can propagate.

Consider the situation shown in Fig. 5. Cell 2 is at rest while cell 1 starts to burst at time t = 0. Under which condition is cell 2, initially in the quiescent phase, going to burst upon the excitation of cell 1 ? As cell 1 starts to burst at time 0, from [START_REF] Bántay | Self-organization and anomalous diffusion[END_REF], the concentration of Ach it produces is

A 1 (t) = A - 1 e -µ A t + β A t 0 T A (V 1 (s))e -µ A (t-s) ds, where A - 1 = β A T A (V Si 1 ) µ A
is the Ach concentration for cell 1 at rest (note that V Si1 depends on G S1 , G A1 via ( 11)). Now, using the mean-field approximation (ii) above we can write

t 0 T A (V 1 (s))e -µ A (t-s) ds ∼ Ω(G S 1 ,G A 1 ) µ A [ 1 -e -µ A t ] so that: A 1 (t) = A - 1 e -µ A t + β A Ω(G S1 , G A1 ) µ A 1 -e -µ A t . (23)
Note that G S1 depends on cell 1's state whereas G A1 depends on the state of its neighbours, here cell 2, which is initially in the quiescent phase, so that

G A1 ∼ 0.
From (23) A 1 increases, thereby increasing the Ach current seen by cell 2.

Without noise, the condition for cell 2 to burst is given by eq. ( 10) characterizing the crossing from region D to region C. Indeed, without noise a cell at rest in D cannot reach the limit cycle. Now, the Ach conductance G A2 seen by cell 2 is, under the approximation (i) above, G A2 (t) = g A 2 √ γ A A 1 (t). Thus, eq. ( 10) reads:

- g A 2 √ γ A A 1 (t B2 ) ( V Si2 -V A ) -G S2 ( V Si2 -V K ) = I SN
where t B2 is the time where cell 2 starts to burst by crossing the bifurcation line. Assuming that A - 1 is so small 5 that we can neglect it, gives:

t B2 = - 1 µ A log 1 + 1 g A 2 µ A √ γ A β A Ω I SN + G S2 (V Si2 -V K ) V Si2 -V A . (24)

Propagation speed

The time t B2 gives the time for a rest cell, excited by its neighbour, to burst. Now, to numerically compute the wave propagation speed, as in Fig. 10, bottom, we actually consider waves of calcium activity i.e. a cell is considered as bursting if its calcium concentration reaches a certain threshold (see section 4.2). The time to reach this threshold from the moment where cell 2 starts to burst is a time t C weakly depend on g A but that can be estimated numerically to a few hundreds of milliseconds. The inverse time t C + t B2 gives therefore the speed of propagation, in this specific geometry, and in the deterministic case:

c = a t C -1 µ A log 1 + 1 g A 2 µ A √ γ A n β A Ω I SN +G S 2 (V Si 2 -V K ) V Si 2 -V A . ( 26
)
This corresponds to the theoretical curve displayed in Fig. 10 bottom. Note that this equation is valid only in the ballistic case, because we have assumed a deterministic propagation without wave interactions. This corresponds to regime II.

Specific propagation regimes

Let us now consider several interesting limits for the transport equation.

Flat sAHP landscape and Kardar-Parisi-Zhang equation

We consider the situation where G S is homogeneous in space (flat sAHP landscape). This characterizes e.g. the propagation of a wave in a region which has not seen a previous propagation since a time larger than the characteristic time of return to rest. The terms ∇G S and ∆ G S vanish and the transport equation ( 22) reads now:

∂G A ∂t = g A D A ∂ 2 Ω ∂G 2 A ∇G A 2 + g A D A ∂Ω ∂G A ∆ G A + S(G S , G A ) (27)
Here G S acts therefore as a parameter changing the characteristics of the transport coefficients. Especially it constrains the diffusion coefficient:

D ≡ D(G S , G A ) = D A ∂Ω ∂G A ( 28
)
∂G A ∂t = g A ∂D ∂G A ∇G A 2 + g A D∆G A + S. (29)
which has the form of the Kardar-Parisi-Zhang equation (KPZ) [START_REF] Kardar | Dynamic scaling of growing interfaces[END_REF].The KPZ equation models interface growth: the first term mimics surface growth in a direction normal to the surface whereas the second term corresponds to diffusive motion along the surface. However, in contrast to KPZ the terms

D(G S , G A ), ∂D(G S ,G A ) ∂G A
are not constant. Additionally, the source term is white noise in KPZ whereas, in our case this is a function of G A , G S . This suggests interesting questions that could be addressed in further studies. Are there scale invariant solutions ? Assume that G S has such a scale invariance, what are the consequences for G A ? Is it possible to characterize (anomalous) tranport in this context ?

Transition between regime I and II

In regime I waves are weakly interacting with the sAHP landscape so that one can consider it as uniform along the path of the waves. This assumption leads to an explicit g A value for the transition between regime I and II.

Let us approximate, in region D, Ω(G S , G A ) by a linear function of

G A , Ω(G S , G A ) = Ω 0 +Ω 1 (G S ) G A with Ω 1 (G S ) ≡ ∂Ω ∂G A (G S ) > 0.
As G S is assumed spatially uniform here, the transport equation ( 19) reads, using [START_REF] Buendía | Hybridtype synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together[END_REF]:

∂G A ∂t = P -λ G A + κ ∆G A with P = g A d β A √ γ A Ω 0 , λ = µ A -g A d β A √ γ A Ω 1 , κ = g A β A Ω1 2 √ γ A a 2 .
Taking the spatial Fourier transform gives:

∂ ĜA ∂t = P δ(k) -σ(k) ĜA with σ(k) = λ + κ k 2 .
Solutions are stable if σ(k) > 0 and marginally stable if

σ(k * ) = 0.
This last condition reads:

k 2 = 2d a 2 1 - 1 g A µ A √ γ A d β A 1 Ω 1 (30)
and has a solution if

g A ≥ µ A √ γ A d β A 1 Ω1 . Let us interpret this condition. If g A < µ A √ γ A d β A 1
Ω1 then σ(k) > 0, ∀k, so that linear solutions are damped and converge to P σ(0

) = g A d β A µ A √ γ A Ω 0 .
This corresponds to a uniform quiescent phase. In this range of value of g A there is no wave propagation so that

µ A √ γ A d β A 1
Ω1 corresponds to the critical value g A0 , eq. ( 25) found in section 5.2.1 in the case of two cells.

There are slight differences though. While [START_REF] Daw | How Vision Works: The Physiological Mechanisms Behind What We See[END_REF] was holding for two cells without linear approximation on Ω, here, we find the critical value without constraint on dimension, but assuming that Ω is linear in G A .

Using this we can rewrite the positive solution [START_REF] Drossel | Complex scaling behavior of nonconserved self-organized critical systems[END_REF]

as k * ≡ k * (g A ) = √ 2d a 1 - g A 0 g A . This defines a characteristic length: Λ(g A ) = a √ 2d 1 1 - g A 0 g A , (31)
which interprets as follows. Consider cells in the quiescent phase and excite a connected cluster of cells, with radius r. If r < Λ(g A ) this perturbation is damped, whereas for r ≥ Λ(g A ) it is amplified, getting rapidly (exponentially fast) out of the range of G A where the linear approximation on Ω holds. This corresponds to drawing cells to region C. Thus, Λ(g A ) is a critical size for a cluster of excited cells to nucleate a wave.

As a consequence, the value of g A such that Λ(g A ) = a corresponds to a regime where exciting one cell is enough to propagate a wave in a uniform landscape of G S . This value is g

A = g A 0 1-1 2d
, so that, in one dimension, g A = 2 g A0 . This is the value drawn in Figs. 5, 7, 10, at the vertical dashed line separating regime I and regime II. Note that this theoretical value is obtained in a purely deterministic setting (without noise).

"Self-Organized" Criticality and beyond

Let us now explore another regime. In the zero noise limit, the transition to bursting can be only achieved when a cell switches from region D to region C. This corresponds to a sharp threshold condition where Ω(G S , G A ) takes the

form π A H [ G A -G Ac (G S ) ],
H being the Heaviside function and G Ac the bifurcation value given by eq. ( 10). Then, eq. ( 19) has the form of a continuous state sandpile, proposed by several authors [START_REF] Zhang | Scaling theory of self-organized criticality[END_REF][START_REF] Bántay | Self-organization and anomalous diffusion[END_REF] in the context of Self-Organised Criticality [START_REF] Bak | Self-organized criticality: An explanation of 1/f noise[END_REF][START_REF] Bak | How nature works: the science of self-organized criticality[END_REF]56,[START_REF] Christensen | Complexity and Criticality[END_REF][START_REF] Marković | Power laws and self-organized criticality in theory and nature[END_REF][START_REF] Buendía | Feedback mechanisms for self-organization to the edge of a phase transition[END_REF]. More generally, the presence of noise smooths the Heaviside function leading to a sigmoid of the form displayed in Fig. 12, allowing, in SOC models, power expansions and dynamical renormalisation group analysis [START_REF] Diaz-Guilera | Dynamic renormalization group approach to selforganized critical phenomena[END_REF][START_REF] Volchenkov | Quantum field theory renormalization group approach to self-organized critical models: The case of random boundaries[END_REF].

In comparison with a SOC model there are some important differences though. First, we have, in our case, an extra term S violating the local conservation required to reach a Self-Organized Critical regime [START_REF] Hwa | Dissipative transport in open systems: An investigation of self-organized criticality[END_REF]56,[START_REF] Vespignani | How self-organized criticality works: A unified mean-field picture[END_REF][START_REF] Cessac | Self-organized criticality and thermodynamic formalism[END_REF].

Actually, S behaves in a different way in the 3 regimes as shown in Fig. 11.

Second, in contrast to early SOC models the "dissipation" does not occur at the boundary but in the bulk, similarly to extended SOC models studied in [START_REF] Bak | Self-organized criticality in non-conservative models[END_REF][START_REF] Hwa | Dissipative transport in open systems: An investigation of self-organized criticality[END_REF][START_REF] Hwa | Avalanches, hydrodynamics, and discharge events in models of sandpiles[END_REF][START_REF] Grinstein | Generic scale invariance and self-organized criticality[END_REF][START_REF] Vespignani | Driving, conservation, and absorbing states in sandpiles[END_REF][START_REF] Drossel | Complex scaling behavior of nonconserved self-organized critical systems[END_REF][START_REF] Cessac | Self-organized criticality and thermodynamic formalism[END_REF][START_REF] Malcai | Dissipative sandpile models with universal exponents[END_REF]. Indeed, S is negative on average, leading to a small bulk dissipation. Finally, the wave activation in our case is not "adiabatic" as in SOC model where one waits that an avalanche stops before activating a new one, thereby preventing avalanches interactions. Thus, we don't have the infinite time scale separation necessary to reach SOC [START_REF] Grinstein | Conservation laws, anisotropy, and "self-organized criticality" in noisy nonequilibrium systems[END_REF][START_REF] Socolar | On self-organized criticality in nonconserving systems[END_REF][START_REF] Grinstein | Generic scale invariance and self-organized criticality[END_REF]. In contrast, several retinal waves can propagate simultaneously and strongly interact, as illustrated in Fig. 8. These interactions limit the size of each wave in the regime II.

In a series of recent papers, Buendia et al report on the different scenarios leading to SOC-like behaviours [START_REF] Buendía | Feedback mechanisms for self-organization to the edge of a phase transition[END_REF][START_REF] Buendía | Hybridtype synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together[END_REF]. Even if "perfect" SOC is not reached there exist observed regimes close to it, such as "self-organized bistability" (SOB), "self-organized quasi-criticality" (SOqC), "self-organized collective oscillations" (SOCO). In addition, the paper [START_REF] Buendía | Hybridtype synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together[END_REF] exhibits, in Fig. 1, a bifurcation diagram quite similar to ours although it deals with a different model (Kuramoto) and different parameters. Also, the paper [START_REF] Buendía | Feedback mechanisms for self-organization to the edge of a phase transition[END_REF] shows, in Fig. 2, a behaviour for the equivalent of our source term quite similar to ours. This suggests that our model, and especially equations ( 19), [START_REF] Chen | A deterministic critical forest fire model[END_REF], may exhibit different regimes corresponding to different modalities of transport equations reported in [START_REF] Buendía | Feedback mechanisms for self-organization to the edge of a phase transition[END_REF] especially "SOC like" behaviours. The fact that we do not observe anything striking in this direction with our simulations may just signify that one has to consider other settings. That is:

1. Dimension. We have only considered 1 dimensional simulations here.

One could happen in 2 dimensions ? We observed that the slow hyperpolarization current following the bursting phase generates a landscape in which subsequent waves have to propagate. In particular, a cell in the absolute refractory state cannot be excited by its bursting neighbours. This way, absolute refractory cells constitute temporary forbidden zone where waves cannot go through, thereby compartmentalizing space with slowly varying boundaries. Clearly, these boundaries are quite different in one and in two dimensions. In 1D, these boundaries are just points; in 2D these are lines with complex geometry. As the topological structure of active and refractory regions could play an important role in shaping the visual system, it is important to consider them in detail. This actually means, as outlined in this paper, that it is not sufficient to consider indicators for the global activity of the network, it is also essential to consider the shape and size of connected regions where activity holds. Unfortunately, this is easy to do in 1D, but is complex in 2D. (Fortunately, SACs are physiologically organized in a two dimensional layer, so, a priori, there is no need to consider 3 dimensional lattices). The detailed study of this question deserves thorough investigations beyond the scope of this paper. Yet, we want to make some comments coming out from the present analysis. phase, evolving on a very slow time scale, so an analogy with spin glasses could be more relevant [START_REF] Binder | Spin glasses: Experimental facts, theoretical concepts, and open questions[END_REF][START_REF] Mezard | Spin Glass Theory and Beyond[END_REF]. In particular, it would be interesting to investigate if evidences of ageing appear in the model when taking the thermodynamic limit. This might be a pure academic question, however, as the meaning of a "thermodynamic limit" in a retina remains to be Finally, we considered here nearest neighbours interaction. Although, this coupling shouldn't be taken "a la lettre" as a real connectivity but an effective one [START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF] one may be willing to consider more realistic connectivities. In particular, using random graphs (Erdos-Reyni) dramatically impacts the phase diagram of the Kuramoto model [START_REF] Buendía | Hybridtype synchronization transitions: Where incipient oscillations, scale-free avalanches, and bistability live together[END_REF], exhibiting a range of parameters with a Griffith phase. Griffith phase is characterized by power-laws extending over broad regions in parameters space. In mag-enhancing their activity, as shown in Fig. 5 and activity spreads through the whole lattice. Thus, in contrast to classical percolation, here cells interaction is active, potentially leading to a different "critical" behaviour (critical exponents) than in directed percolation or forest fires. Note also that the transition corresponds to a determined value of g A whereas g A is evolving during development, so that the "percolation" phase is a brief instant. Thus, regime II, extending on a whole interval (Griffith phase), might be more relevant for visual system development.

4. Self-organization. Assume that their exist some notion of criticality in retinal waves how could the retina manage to maintain a critical state as g A is evolving ? Several biological mechanisms have been proposed in the literature to explain how the brain could achieve to maintain a critical state: synaptic adaptation [START_REF] Levina | Dynamical synapses causing selforganized criticality in neural networks[END_REF][START_REF] Levina | Phase transitions towards criticality in a neural system with adaptive interactions[END_REF][START_REF] Bonachela | Selforganization without conservation: are neuronal avalanches generically critical[END_REF][START_REF] Wang | Hierarchical modular structure enhances the robustness of self-organized criticality in neural networks[END_REF][START_REF] Di Santo | Landau-ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization[END_REF][START_REF] Kinouchi | Stochastic oscillations and dragon king avalanches in self-organized quasi-critical systems[END_REF] , dynamic neuronal gains [START_REF] Brochini | Phase transitions and self-organized criticality in networks of stochastic spiking neurons[END_REF][START_REF] Kinouchi | Stochastic oscillations and dragon king avalanches in self-organized quasi-critical systems[END_REF], adaptive firing thresholds [START_REF] Girardi-Schappo | Synaptic balance due to homeostatically self-organized quasicritical dynamics[END_REF], topological self-organization [START_REF] Bornholdt | Topological evolution of dynamical networks: Global criticality from local dynamics[END_REF] (see [START_REF] Kinouchi | Mechanisms of self-organized quasicriticality in neuronal networks models[END_REF] and references therein for a recent review), homeostasy [START_REF] Menesse | Homeostatic criticality in neuronal networks[END_REF].

We first would like to remark that most of these models where made for cortical, spiking neurons. In contrast, most neurons in the retina, including SACs, do not spike6 . Also, there is no evidence of long term plasticity although short term plasticity has been reported [START_REF] Kastner | Adaptation of inhibition mediates retinal sensitization[END_REF] (dealing with non spiking neurons). It is therefore not clear to us which mechanism could apply. Note however that the most important aspect revealed by our paper, making dynamics so rich, is that region D is a recurrent region close to bifurcations. Dynamics naturally drives the cell back to this region where it is quite sensitive to perturbations thanks to basic mechanisms observed in experiments. In this way, the mechanism of "critical stability" could be close to the scenario called "Mapping Self-Organized Criticality onto Criticality" introduced by D. Sornette and co-workers in [START_REF] Sornette | Mapping self-organized criticality onto criticality[END_REF]. In the increases sharply as well as waves characteristics (size, duration) and network activity (n). The separation between regime I and II holds at the unique point where ∂ρ ∂g A is maximum, that is the point where the probability ρ that a cell bursts at a given time has the maximum sensitivity to variations of g A . The transition corresponds to a characteristic scale of wave nucleation equal to the SACs spacing. As g A further increases, in regime II, waves are essentially propagating deterministically, with a decreasing influence of noise, corresponding to cells in-between region D and C. Transport is ballistic there. What mainly limits the wave size is the interaction with the sAHP landscape left by previous retinal waves. This is therefore a regime of competition between waves. The probability ρ that a cell is bursting at time t as well as waves characteristics are increasing essentially linearly with g A , with a steep slope. In the transition between regime II and regime III one switches from a regime of competition to a regime of cooperation (non linear feedback) prolonging burst durations and sAHP. At the transition, one typically have several cooperating waves with a global activity spreading through the whole lattice. Hence, the transition between regime II and III resembles a percolation transition, as formerly pointed out by several authors [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF][START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF]. However, in contrast to classical percolation, non linear coupling cells enhances cooperative effects. In regime III, waves size and duration saturate and ρ increases linearly, with a small slope, in the range of g A values we investigated. For larger g A we expect ρ to saturate and to converge to the fraction of time a cell can be excited, dictated by sAHP parameters. We didn't explore this regime though as we were interested in ranges of g A values compatible with retinal development. In this perspective, actually, g A is decreasing, presumably leading from regime III to regime I, which, according to our estimation of table 1 could coincide with the switch from cholinergic-SAC (stage II) with glutamatergic-AII (stage III) retinal waves.

In this perspective, transition between regimes such as the percolation-like transition from regime III to regime II is a brief instant in the evolution and one does not see why it should play an important role. Nevertheless, the question of critical behaviour in retinal waves is important. Since the seminal work of Beggs and Plen [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF] -reporting that neocortical activity in rat slices occur in the form of neural avalanches with power law distributions close to a critical branching process -there have been numerous papers suggesting that the brain as a dynamical system fluctuates around a critical point. Scale-free neural avalanches have been found to occur in a wide range of neural tissues and species [START_REF] Shew | The functional benefits of criticality in the cortex[END_REF]. From a theoretical point of view, it has been suggested [START_REF] Haldeman | Critical branching captures activity in living neural networks and maximizes the number of metastable states[END_REF][START_REF] Shew | The functional benefits of criticality in the cortex[END_REF] that such a scale-free organisation could favor information storage and transfer, improvement of the computational capabilities [START_REF] Bertschinger | Real-time computation at the edge of chaos in recurrent neural networks[END_REF], information transmission [START_REF] Beggs | Neuronal avalanches in neocortical circuits[END_REF][START_REF] Bertschinger | Real-time computation at the edge of chaos in recurrent neural networks[END_REF][START_REF] Shew | Information capacity and transmission are maximized in balanced cortical networks with neuronal avalanches[END_REF], sensitivity to stimuli and enlargement of dynamic range [START_REF] Kinouchi | Optimal dynamical range of excitable networks at criticality[END_REF][START_REF] Shew | Neuronal avalanches imply maximum dynamic range in cortical networks at criticality[END_REF][START_REF] Gautam | Maximizing sensory dynamic range by tuning the cortical state to criticality[END_REF][START_REF] Girardi-Schappo | Griffiths phase and long-range correlations in a biologically motivated visual cortex model[END_REF].

It has been proposed by Hennig et al [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF] that such a critical regime could also be helpful for the visual system developmental phase triggered by stage II retinal waves. In such a regime the retina would maximize randomness and variability in the resulting activity patterns. However, they do not propose a mechanism explaining how this regime can be robust to parameters variations during development, especially cholinergic coupling. As we discussed in the paper, such a structural stability could correspond to a Griffith phase (also observed in V1 [START_REF] Girardi-Schappo | Griffiths phase and long-range correlations in a biologically motivated visual cortex model[END_REF]) instead of a classical critical regime. This interesting possibility would however deserve extensive simulations in 2D, with a realistic connectivity, and a thorough analysis data to cleanly identify a critical signature. There is a wide literature dealing with the precise characterization of a critical state, where scaling exponents obey relations fixed by renormalization group theory [START_REF] Ma | Modern theory of critical phenomena[END_REF][START_REF] Hohenberg | Theory of dynamic critical phenomena[END_REF][START_REF] Kadanoff | Scaling and universality in avalanches[END_REF][START_REF] Sethna | Hysteresis, avalanches, and disorder induced critical scaling: A renormalization group approach[END_REF]. In the field of neuroscience this has been widely discussed by authors such as Touboul and Destexhe [START_REF] Touboul | Can power-law scaling and neuronal avalanches arise from stochastic dynamics[END_REF][START_REF] Touboul | Power-law statistics and universal scaling in the absence of criticality[END_REF] and many authors have looked for solid evidences of criticality [START_REF] Girardi-Schappo | Critical avalanches and subsampling in map-based neural networks coupled with noisy synapses[END_REF][START_REF] Moretti | Griffiths phases and the stretching of criticality in brain networks[END_REF][START_REF] Girardi-Schappo | Griffiths phase and long-range correlations in a biologically motivated visual cortex model[END_REF]. In particular, it is believed that the so-called crackling noise relationship (from Sethna's terminology [START_REF] Sethna | Hysteresis, avalanches, and disorder induced critical scaling: A renormalization group approach[END_REF]) can be used to discriminate between true critical systems and "fake" ones. Recently, Kanders et al [START_REF] Kanders | Fingerprints of a second order critical line in developing neural networks[END_REF] have shown that developing neural cultures in specific conditions satisfy the crackling noise relationship. The same type of investigations could be done with our model.

We would like now to address a few further open questions related to this work.

Spontaneous Bursting

It is possible to move regions C, D downward or upward by a slight change in the leak potential V L . Experimentally, is it possible to modify V L by changing the outside potassium concentration. Note, however, that the change in a reversal potential is logarithmic in the concentration. Nevertheless, in our case, a change of V L from -72 mV, the value used in the paper, to V L = -70 mV has the effect to move region D to negative values for G A . This means that, in this case, when the cell returns back from its sAHP excursion it eventually penetrates in region C, instead of region D. Then, it starts to oscillate again (stable limit cycle in C). Therefore, in this situation, no noise is required to restart bursting. It occurs spontaneously.

In this scenario, cells burst therefore periodically, with a frequency controlled by τ R , τ S . An example is shown in Fig. 14. Here, cells have a tendency to synchronise. However, introducing a bit of randomness in initial conditions one can obtain, depending on the maximal conductances g S , g A , a regime with spatio-temporal disorder where waves interacts and leads to random regions similar to the Noise Induced Bursting regime [START_REF] Karvouniari | Mathematical and experimental studies on retinal waves[END_REF].

The cholinergic interactions of SACs in a quenched disordered sAHP landscape breaks down the spontaneous bursting periodicity and leads to a spatiotemporal chaotic dynamics quite similar to the noise induced bursting for waves dynamics [START_REF] Karvouniari | Mathematical and experimental studies on retinal waves[END_REF]. We didn't address in detail the spontaneous bursting regime in this paper though leaving it for further studies.

Understanding the functional role of retinal waves

It would be interesting to explore, at a modelling level, what could be the functional impact of the different waves regimes we observe in the development of the retino-thalamico-cortical pathways, elucidating why early neural networks would choose to maintain their activity in such regimes. As we have seen in section 4.1 experimental evidences show that g A decays during development so that, according to our results, the retina ought to go successively through regimes III (extended waves), II (wide distribution of waves and strong interaction with the sAHP landscape), and I (localized waves) although this step in evolution coincide with the arousal of stage III waves. More particularly, we reported the apparition of sAHP patterns that spatially bound the propagating waves. Those localized bounding patches are found to persist more than a sAHP timescale, indicating that in the network scale, they introduce a type of spatial memory of previous activity. In other words, they create an heterogeneous landscape where waves could propagate, inducing an effect, first on waves characteristics, and second on the selective synaptic shaping of the network through this spatial bias 7 For example, late waves (stage III) appear to be spatially bounded and more localised than stage II waves [START_REF] Maccione | Following the ontogeny of retinal waves: pan-retinal recordings of population dynamics in the neonatal mouse[END_REF], suggesting that probably late stage II waves could prepare the landscape for the transition to more spatially confined activity. Hence, those localized activity patterns, could potentially have a link with how receptive fields are formed before vision becomes functional. A next step would therefore be to link the retinal waves model presented here to a cortical V1 model such as the one developed in [START_REF] Souihel | Anticipation in the retina and the primary visual cortex : towards an integrated retino-cortical model for motion processing[END_REF] so as to investigate the role of synaptic shaping under retinal wave stimulation to drive the cortical response during development.

This work has also possible future outcomes with respect to retina therapy as well. Understanding how retinal waves are initiated and propagate in the retina could enable one to define protocols to trigger such retinal waves in the in vivo adult retina. Inducing such waves is expected to reintroduce some plasticity in the retinal tissue and the projections in the brain. This induced plasticity could have important therapeutic applications to treat patients or stimulate regeneration of retinal ganglion cell axons following optic nerve crush.

Parameters Value

Parameter

Figure 1 :

 1 Figure 1: 1. Top. Bifurcation diagram of the fast dynamics (1)-(3) in the plane G S , G A .

Figure 2 :

 2 Figure 2: Typical course of the variables C, R, S for a cell. in the presence of noise.

Figure 3 : 2 .

 32 Figure 3: Noise induced bursting. Both figures show the join evolution of voltage (upper trace) and calcium concentration (lower trace) for η = 10 pA ms 1 2 . Left. Region D. Right. region A. In region D the trajectory stays long enough onto the stable periodic orbit to trigger a high increase in calcium concentration. In region A noise can make the trajectory go in the region delimited by the saddle and its unstable manifold, leading to a short spike in V , not long enough to let the calcium concentration grow above the threshold θ (black line).

Figure 4 : 2 .

 42 Figure 4: Probability that a cell enters in the bursting regime where calcium exceeds the threshold 4 C (0) , in the space G S , G A , for η = 10 pA ms 1 2 . Left. 3 D view of the probability with color map. Right. Color map of the bursting probability, plotted together with the Homoclinic bifurcation line (dashed orange) and the SN 2 bifurcation line (dashed black). This is a zoom of region D (which extends further in the directions G S and G A ). Black corresponds to a 0 probability while yellow corresponds to a probabilty 1. The intermediate scales are shown in the color bar.

Figure 5 :

 5 Figure 5: . Left. Evolution of Ach concentration for 2 coupled cells, during a burst, for different values of g A in nS. Row 1 corresponds to cell 1 and row 2 to cell 2. Here, cell 1 bursts first due to a small pulse in V while cell 2 is driven by cell 1: it bursts due to the Ach current generated by 1. Right. Burst duration, τ 1,max , and Ach peak, A 1,max , of cell 1 under the feedback of cell 2. Region I,II,III are explained later in the paper.

540D

  which is moving since the value of G S , G A slowly evolve in time, with different time scales. Wave propagation in this regime is essentially driven by noise.

3. 4 .CFigure 6 :

 46 Figure 6: Top. Example of trajectory of a bursting cell in the plane {G S , G A } for, from left to right, g A = 0.04, 0.08, 0.1, 0.2 nS. Note the oscillations in G A coming from the bursting of neighbours even if the cell itself is not bursting (e.g. region A). Green and orange lines respectively correspond to SN and Hc bifurcation. Bottom. Probability density of cells state in the plane {G S , G A } for the same g A values as in Fig. 8. Color bar on top indicates the corresponding value of probabilities (estimated on boxes with size 1 pS 2 in the plane {G S , G A }). Note that the scales are different in the different figures. Green and orange lines respectively correspond to SN and Hc bifurcation.

  production is high enough to provoke cells synchronization and waves. As g A grows, cells have a decreasing probability to be in region D and wave propagation becomes less sensitive to noise. Finally, the right most column (regime III) 610 corresponds to cells mostly located in region C, with high Ach production, inducing a high global activity.

615 4 . 1 .

 41 Fig 3B in Zheng et al, [109]. They puff Ach with a concentration 1 mM while they block synaptic transmission with C 2+ d . They monitor the induced Ach current and plot its maximum (arising after ∼ 2s) as a function of the developmental time. They observe a clear decay. From this figure we extrapolated the table 1 below and the plausible behaviour of g A .

  where n is the number of activated cholinergic synapses connecting a SAC, V clamp = -70 mV (patch clamp in Zheng et al experiment) and V A = 0 mV. Thus g A = -I puf f n V clamp .The number n is difficult to estimate from experiments. From[START_REF] Zheng | A transient network of intrinsically bursting starburst cells underlies the generation of retinal waves[END_REF] we may fix n ∼ 25 -30 so that g A ∼ I puf f 2000 , which gives us hints on the conductance per synapse (assuming actually that n is constant during development).

695 4 . 3 . 1 .

 431 Quantitative description of the regimesRegime I. Localized waves. When g A is small, cells can be excited by noise but the Ach they produce does not increase sufficiently the probability of excitation of their neighbours to induce a large scale propagation.

  (a)); regime III (g A = 0.2 nS). The top row shows the time evolution (vertical axis) of calcium concentration in the lattice (horizontal axis) in these different regimes. We can see the initiation and propagation of waves. The bottom row shows the normalized correlation of calcium activity of the center cell (index i c ) with another cell in the lattice (index i). The correlation is numerically

Figure 7 :

 7 Figure 7: Left (a). Probability ρ that a cell is bursting at a given time as a function of g A , for N = 20, 50, 100, 200, 400. The inset shows the numerical derivative of ρ for N = 400. The maximum of the derivative defines the transition between regime I and II (leftmost dashed line). This transition is further explained in section 5. The second dashed line, roughly defines the transition between regime II and III and arises at g A ∼ 0.14 nS. Red points with labels correspond to specific values of g A : one in regime I, two in regime II (II 1 at the transition between I and II and II 2 ) and one in regime III, used for representative plots in the subsequent figures. Center. Average wave size, s , label (b), and duration, D , label (c), in seconds, as a function of g A , for the same values of N (in logscale). Vertical black lines correspond to the separation in 3 regimes. Right (d). Mean activity, n , (see text for a definition) in the same conditions. Note that size and duration do not depend on N , the lattice size, in contrast to n, which, actually grows further for larger values of g A (not shown).

Figure 8 :

 8 Figure 8: Spatial characterization of the 3 wave regimes. Top. Propagation of waves in the presence of noise (η = 6 pA ms -1 2 ) for 200 cells illustrated by the evolution of calcium concentration. From left to right: g A = 0.04 nS (regime I); g A = 0.08 and g A = 0.1 nS (regime II, point II 1 , II 2 in Fig. 7); g A = 0.2 (regime III). Hot zones (orange, yellow) are active (bursting) zones; black zones are absolute refractory zones; purple zones correspond to SACs in the rest state. Propagating waves are stopped by absolute refractory regions. Bottom. Normalized calcium correlation C (i) C (ic) between the central cell and other cells for the same values of g A . There are zones of positive correlation and negative correlations (see text).

Figure 9 :Figure 10 :

 910 Figure 9: Probability distributions of waves quantities for different values of N in the 3 regime corresponding to points I, IIa, IIb, III in Fig. 7 (from left to right). Top. Sizes distribution; Middle. Durations distributions; Bottom. Activities distributions.

Figure 11 :

 11 Figure 11: Histogram of the source term, S, in equation (19) in the 3 regimes (from left to right).

  region. In the intermediate region the diffusion is very large, thus any spatial inhomogeneity disappears quite fast. Consequently, the inflow flux is transported through the intermediate region to the quiescent phase boundary increasing the

Figure 12 : 2 . 1 ..> 0 .

 12210 Figure 12: Top, left. Piecewise linear approximation for U with γ A = 1 nM 1 2 . The dashed line corresponds to 2 √ γ A . Top, right. Average value Ω ≡ T A as a function of G S , G A . Note the level lines, at the bottom, and compare to the bifurcation diagram Fig. 1. Middle row, left. Gradient of Ω, ∂Ω ∂G S . Middle row, right. Gradient of Ω, ∂Ω ∂G A. Gradients are essentially non

Figure 13 :

 13 Figure 13: Left. Trajectory of R, S during the cycle bursting, hyperpolarization, repolarisation. "FP" means "Fixed Point". Right. Trajectory of G S (τ S , τ R = 8.3 s) in the same conditions. Note that the shape can be more complex depending on parameters such as g sAHP .

  clarified. Another question coming out from the analogy with an Ising model is "what could play the role of the external field". Photoreceptors become more and more active during the first weeks of development and it could be interesting to investigate if they could play a role, in stage II, by acting on Bipolar cells and thereby on Starburst Amacrine cells.

Figure 14 :

 14 Figure 14: Spontaneous bursting (SB) regime. Voltage Rhythmic bursting for one isolated SAC as a function of time. The interburst depends on τ S , τ R .

  [START_REF] Abel | Relationships between intracellularneocortical pyramidal neurons calcium and afterhyperpolarizations in neocortical pyramidal neurons[END_REF]. More details are given in the legend.

Additional information, not necessary for the understanding of the paper, are available on the web page https://team.inria.fr/biovision/bifurcations-map-ofthe-retinal-waves-model/. Especially, the region in the interval

∼ ( G S , G A ) ∈ [2.2, 3.2[×[1.3, 2.7[,

containing the SNH, BT, and SNH points, is not discussed in the paper, for a simple reason: the dynamics of our model, with its parameters grounded on experiments, never enters in this region (see

Fig 6)

. A zoom on this region, with the corresponding bifurcation diagram, can nevertheless, be found at the aforementioned web page. Let us also mention that, in this region, dynamics ought to be very close to the model studied by Moreno-Spiegelberg et al

Table 1 :

 1 Values of g A versus perinatal days, extrapolated from Fig.3Bin[START_REF] Zheng | A developmental switch in the excitability and function of the starburst network in the mammalian retina[END_REF].

	nS)	Comments

and durations distribution. As we show, it strongly depends on the maximal intensity of Ach coupling g A . The values of g A displayed in Table

1

are actually only indicative on the modelling side, as we are considered, in simulations, a nearest neighbours topology different from the real connectivity. Nevertheless, 650 we tune model parameters to have stage II waves in a similar range of values for g A . Especially, the wave speed we observe are in agreement with experimentally measured values

  1.4. The almost-vanishing of the source term, in regime I, expresses that the local degradation of acetylcholine is compensed by

		0.6					0.15					0.1
					0.01 nS					0.09 nS					0.16 nS
		0.5			0.02 nS 0.03 nS					0.1 nS 0.11 nS					0.2 nS 0.24 nS
					0.04 nS					0.12 nS					0.28 nS
		0.4			0.05 nS 0.06 nS		0.1			0.15 nS
	P(S)	0.3			0.07 nS 0.08 nS	P(S)					P(S)	0.05
		0.2					0.05
		0.1
		0					0					0
		-0.002	0	0.002	0.004	0.006	-0.002	0	0.002	0.004	0.006	-0.002	0	0.002	0.004	0.006
				S (nM/s)					S (nM/s)					S (nM/s)

  First, a large part of the mechanism of retinal waves can be explained by local dynamics and does not depend on dimensionality. Yet, the synchronisation mechanism is a non local effect, depending, as we saw, on the non linear coupling intensity g A , sAHP conductance g S and noise η, fixing the characteristic size of active and refractory regions. It is also constrained

	dritic tree is rather short compared to the size of the retina (2 mm in
	the mouse, compared to 50 µm for the dentritic tree radius of a SAC),
	so SACs interaction can be considered as local. Yet, Ising model, taught
	us that local interactions can induce long range correlations at a phase
	transition point. However, phase transitions in Ising model, and more
	generally, models with short range interactions do not appear in 1D. In
	this sense too, the behaviour of retinal waves could be quite different in
	1D and 2D.
	by boundary conditions. Especially, zero boundary conditions imposes
	quiescent cells at the boundaries, thereby constraining the way how waves
	propagate and interact and their characteristic size. Here, we can make
	a formal analogy with statistical mechanics and Ising model. SACs den-

2. Potential phases transition and connectivity. The formal analogy with Ising model actually raises further comments and questions. First, dynamics here is quite more complex than an Ising model (or classical forest fires) as couplings are slowly evolving in time, and depend non linearly in dynamics. Refractory domains constitute somewhat a glassy

Table 2 :

 2 Typical parameters values used in the model.

Recent experimental works by E. Sernagor lab, suggest that some specific cells, whose type

Even if the real SACs connectivity is known at maturation[START_REF] Tauchi | The shape and arrangement of the cholinergic neurons in the rabbit retina[END_REF], we don't know about papers reporting how connectivity evolves during development.

√ γ A β A Ω µ A ( V Si 2 -V A )

Neural spikes are produced by sodium and potassium fluxes, whereas it was shown by Zhou et al[START_REF] Zheng | A developmental switch in the excitability and function of the starburst network in the mammalian retina[END_REF] that sodium plays no role in SACs bursting.

This spatial heterogeneity might also have an impact on blood vessels organisation with, thus, a very long term influence.

Acknowledgements

We thank the reviewers for insightful remarks and constructive criticism. D.M.K. was supported by a doctoral fellowship from Ecole Doctorale des Sciences et Technologies de l'Information et de la Communication de Nice-Sophia-Antipolis (EDSTIC) and by the RENVISION European Union Project No. 600847. This work also benefited from the support of the Neuromod institute of the University Côte d'Azur and the Doeblin federation. We warmly acknowledge Lionel Gil, Matthias Hennig, Evgenia Kartsaki, Evelyne Sernagor, Olivier Marre and Serge Picaud for their help.

This equation is sufficient to characterize the wave propagation in the geometry we are considering in this paper: one dimensional and nearest neighbours interactions. Note however that is generalizes to a cell at rest surrounded by n bursting neighbours whose Ach concentration increases according to [START_REF] Choi | Intrinsic bursting of aii amacrine cells underlies oscillations in the rd1 mouse retina[END_REF] provided these neighbours started to burst at the same time (synchrony). Let us now comment on this result.

1. This equation is actually implicit. Indeed, Ω depends on G S1 , G A1 which depend on t. Thus, Ω actually depends on t B2 via G S1 , G A1 . Likewise, V Si2 depends on G A1 . However, using the methodology developed in this paper, we may consider that t B2 is a function of the parameters G S1 , G A1 which slowly move with time in the bifurcation map as shown in Fig. 6 (growing of Ach conductance upon bursting of neighbours).

2. Likewise, t B2 depends on G S2 , characterizing the level of hyperpolarization of cell 2. The higher G S2 the longer it is for cell 2 to eventually burst, until a value of G S2 where t B2 becomes infinite.

3. More precisely, for eq. ( 24) to have a solution with a positive t B2 one needs that -1

As

such that the cell 2 is excited by cell 1 only if g A > g A0 . This value depends on G A1 , G S1 via Ω, and also on G S2 . The most important here is that g A0 depends on G S2 , the level of hyperpolarisation of cell 2. As G S2 increases g A0 increases so that the conductance g A necessary to trigger excitation propagation to a quiescent cell depends of its level of hyper-polarisation, as expected. For a cell at rest, we found, with our parameters value (Table 2) g A0 ∼ 0.04 nS. This corresponds to the blue dashed line in Figs. corresponding to a percolation phase transition", where waves of activityoften referred to as "avalanches"-are distributed according to power laws (see Fig. 4 of [START_REF] Hennig | Early-stage waves in the retinal network emerge close to a critical state transition between local and global functional connectivity[END_REF], although it is not completely clear if this is a power law). They interpret this regime as an indication that "early spontaneous activity in the developing retina is regulated according to the following principle: maximize randomness and variability in the resulting activity patterns". A similar point of view is defended in [START_REF] Lansdell | A reaction-diffusion model of cholinergic retinal waves[END_REF]. The equivalent of a percolation transition in our model is at the transition between regime II and regime III. One switches from a competitive regime where waves interact with the sAHP landscape, preventing them to propagate, to a cooperative regime where waves collides, with the effect of non linearly past, this line of thought has collected a broad consensus and resulted in numerous interesting results [START_REF] Sornette | Mapping self-organized criticality onto criticality[END_REF][START_REF] Levina | Dynamical synapses causing selforganized criticality in neural networks[END_REF][START_REF] Levina | Phase transitions towards criticality in a neural system with adaptive interactions[END_REF][START_REF] Di Santo | Landau-ginzburg theory of cortex dynamics: Scale-free avalanches emerge at the edge of synchronization[END_REF][START_REF] Rubinov | Neurobiologically realistic determinants of self-organized criticality in networks of spiking neurons[END_REF][START_REF] Shew | Adaptation to sensory input tunes visual cortex to criticality[END_REF]]. Yet, the recurrence to region D is not sufficient to ensure criticality, requiring also long range correlations via Ach coupling. A thorough investigation of criticality in our model would reveal whether the dynamics we propose is sufficient or if we need extra mechanisms such as the one proposed in [START_REF] Kossio | Growing critical: Self-organized criticality in a developing neural system[END_REF] where the feedback of the control parameter to the order parameter is carried out through the dynamics of the synaptic tree radius.

Discussion

In this paper, following [START_REF] Karvouniari | A biophysical model explains the oscillatory behaviour of immature starburst amacrine cells[END_REF], we answered several questions with respect to our understanding of retinal waves. Which is the biophysical mechanism that generates sustained periodically bursting in immature SACs and which are the parameters that control it? Is there a potential, universal, mechanism of waves generation yet accounting for the observed waves variability ? How do SACs synchronize in order to produce propagating waves? How do waves propagate ?

How do they stop? How do the characteristics of waves depend on biophysical parameters ?

The

Table 3: Typical parameters values related to the acetylcholine neurotransmitter (Ach).