Gap junctions set the speed and nucleation rate of stage I retinal waves

Spontaneous waves in the developing retina are essential in the formation of the retinotopic mapping in the visual system. From experiments in rabbits, it is known that the earliest type of retinal waves (stage I) is nucleated spontaneously, propagates at a speed of 451±91 μm/sec and relies on gap junction coupling between ganglion cells. Because gap junctions (electrical synapses) have short integration times, it has been argued that they cannot set the low speed of stage I retinal waves. Here, we present a theoretical study of a two-dimensional neural network of the ganglion cell layer with gap junction coupling and intrinsic noise. We demonstrate that this model can explain observed nucleation rates as well as the comparatively slow propagation speed of the waves. From the interaction between two coupled neurons, we estimate the wave speed in the model network. Furthermore, using simulations of small networks of neurons (N≤260), we estimate the nucleation rate in the form of an Arrhenius escape rate. These results allow for informed simulations of a realistically sized network, yielding values of the gap junction coupling and the intrinsic noise level that are in a physiologically plausible range.

Eq. (1-3). We are specifically interested in the propagation mechanism, and thus in the 10 time the voltage needs to go from the vicinity of the resting potential to the peak 11 potential that marks the occurrence of the first spike. Consequently, the reset 12 mechanism Eq. (3) can be discarded. Until this first spike time, the gating variable u 13 can roughly be regarded as constant, cf. the phase space trajectory displayed in 14 Fig. 1(b). Accordingly, we replace u i in Eq. (1) by the constant u r = −19.2 mV., i.e. 15 the resting value of the gating variable. The resulting expression is still a system of 16 coupled differential equations, due to the interaction with neighboring neurons in the 17 network, that is given by the GJ current Eq. (5). To circumvent this problem we simply 18 replace the membrane potential of a bursting neighboring neuron by the constant value 19 V b , an estimate of the temporal average of a neuron's voltage during the burst, which is 20 described in detail below. The voltage of the silent neighbors is replaced by the resting 21 potential. For the discussion of the propagation mechanism these replacements are fairly 22 reasonable, because the propagation process can be thought of as a certain number of 23 synchronously bursting neurons (typically, one or two) exciting a connected neighbor 24 with additional quiescent neighbors. 25 Under these conditions, the dynamics Eq. (1) can generally be recast into the form This equation can be solved by separating the variables and integrating. It is however 27 necessary to distinguish the two cases of V 1 and V 2 being real valued (corresponding to 28 the situation with a stable resting potential) or complex valued (unstable situation). We 29 will refer to real valued (V 1 , V 2 ) as (V r , V c ) with V r < V c . Real values are e.g. obtained 30 for the case of an isolated neuron, or equivalently G = 0, and are related to our original 31 parameters as follows leading to (V r , V c ) = (−64 mV, −60 mV) for our parameter choice. In this case, the 33 integration of Eq. A yields April 12, 2019 1/7 which can be inverted and thus leads to the explicit solution for the voltage trajectory: 35 Turning to the case of complex valued (V 1 , V 2 ), we will refer to them as The explicit expression of V m and γ depend on the specific setup of the neighbors (see To analytically understand the propagation of the burst from cell to cell, we 41 approximate the time-dependent driving voltage of the bursting cell by its time average, 42 V b , resulting in the GJ current: This replacement is justified because during the period of interest (the time needed for 44 the driven cell to generate a spike), the bursting neuron generates at least a few spikes, 45 i.e. changes rapidly compared to the voltage of the driven neuron.

46
The mean voltage during the burst can be estimated by integrating the voltage over 47 one inter-spike interval T ISI , for simplicity considered for an isolated cell (G = 0) that is 48 started at V (t = 0) = V reset , u(t = 0) = u r , where T ISI = t(V peak ) − t(V reset ). Using Eq. C, we find Using Eq. D and Eq. H in Eq. G, we can calculate the integral and further simplify the 51 resulting expression For our standard parameters this givesV b ≈ −34 mV (T ISI ≈ 73 msec).

53
Propagation Speed in the Three-Neurons Approximation

54
Let us now turn to the situation, where an initially quiescent neuron is driven to a burst 55 due to the GJ current that results from one bursting neighbor in a one-dimensional 56 chain. To this end, we consider the three neurons i − 1 (bursting), i (to be excited), and 57 i + 1 (quiescent). Neuron i is brought from the resting potential V r to the peak 58 potential V peak in a period T ISI + T B that approximately consists of one inter-spike 59 interval T ISI (the time needed by neuron i − 1 to go from reset to peak potential for the 60 first time) and the BOTD T B that determines the speed of the wave via During the period t ∈ [−T ISI , T B ] (light blue shaded area in Fig. A(b)), we approximate This equation can be recast into the form of Eq. A with complex V 1,2 = V m ± iγ in the 65 case of propagating bursts, where For the calculation of T B we can then employ Eq. E, yielding A comparison between calculated burst onset times according to Eq. O and T B obtained 68 from the simulation of a one-dimensional chain is shown in Fig. A(c). The 69 approximation of T B (Eq. O, solid line) shows reasonable agreement with the simulation 70 results (symbols), in particular for small values of G. This is so, because our averaging 71 argument is most reasonable at weak coupling between the neurons (similar to standard 72 approximations for weakly coupled oscillators [S1]).. As stated in the main text, we 73 approximate the propagation speed of waves passing through a two-dimensional network 74 of neurons by Eq. 9. This corresponds to the assumption that the wave's front is 75 reasonably well approximated by a planar shape, if far enough from its origin. As 76 displayed in Fig. 3(b), the propagation mechanism can then be mimicked by a 77 one-dimensional situation with rescaled distance and coupling strength.

78
In order to see this, note that neurons of the same color in Fig. 3(b) that are part of 79 a perfectly planar wave front share exactly the same state (V, u). If the voltage of all 80 horizontal neighbors is identical, links between these neurons can be discarded, because 81 the GJ current is zero. We set the first spike time of the bursts of all red neurons as 82 time origin. Now, every single blue neuron feels an excitatory current from two bursting 83 neurons (connected via the links indicated in red/blue). The leak current of this one 84 blue neuron is affected by the two links connecting this neuron to two yellow neurons 85 (indicated by blue/yellow lines in Fig. 3(b)), which are to a good approximation at rest 86 at this instant of time. Doubling the excitatory current and the additional leak current 87 via GJs can be expressed by doubling the GJ conductance parameter G in Eq. O. Last 88 but not least, the propagation of the wave within one burst onset time difference is not 89 in direction of the link, but we have to consider the reduced distance 3/4 · . Putting 90 everything together, we obtain Eq. 9. To discuss the BOTD in an alternative approach, we would like to calculate the velocity 94 of burst propagation using the one-dimensional version of our system in the continuum 95 limit. Following the literature (e.g. [S2]), we construct a traveling front applying 96 singular perturbation theory and numerical shooting.

97
To this end, we discard the dynamics of the slower gating variable u, constraining it 98 to its resting value u r and use the recast expressions introduced above, cf. Eq. K. In the 99 continuum limit, the GJ coupling can be interpreted as a diffusion term, resulting in: April 12, 2019 4/7 For similar reaction-diffusion systems, the velocity of the traveling wave has often been 101 determined using the following ansatz: The parameter s represents the velocity of the traveling wave, such that ξ is a 104 coordinate in a frame moving with velocity s. This second order ODE is transformed 105 into two first-order ODEs, here introducing W = dV dξ , such that Obtaining the front solution (V, W ) T (ξ), and thereby s, requires the definition of 107 boundary conditions. The resulting boundary value problem can be solved by numerical 108 shooting. The idea of this method is to solve a boundary value problem by treating it as 109 an initial value problem (which can simply be evaluated by numerical integration) and 110 change the parameters of the system, until the boundary conditions are fulfilled.

111
The "left" boundary condition in our case is given by the resting state of the neuron 112 and therefore This fixed point is a saddle and, for the numerical shooting method, the system is 114 initialized in its vicinity on the unstable manifold, cf. Sec. 6.2.1 in [S2]. Next, the 115 parameter s is iteratively changed, until the "right" boundary condition is reasonably 116 well recovered. We assume this right boundary condition to be where ξ right denotes the value of ξ at the right boundary, which is not infinite because it 118 is reached after finite time. Note that this does not represent a fixed point as common 119 in other applications of numerical shooting. In the original discrete and discontinuous 120 model Eq. (1-3), V = V peak is the effective threshold upon reaching of which the 121 membrane potential is reset to V reset . There is no such reset in the one-dimensional 122 continuum approximation. To mimic this threshold behavior, we choose V peak as the 123 end point of the numerical shooting.

124
Further, we chose lim ξ→∞ W (ξ) = 0, which may seem rather arbitrary. There is 125 actually a range of values for the parameter s, that lead to lim ξ→∞ V (ξ) = V peak . We 126 observe that there is a lower bound for the wave speed |s|, for which the trajectory of interesting to see this strikingly good agreement and could argue the applicability of 139 numerical shooting methods for studying discrete and discontinuous bursting neuron 140 models, such as ours.

141
Numerical Simulation Methods

142
The numerical simulations of our system were performed according to the following 143 Euler-Maruyama integration scheme:

157
In Fig. 4, we find a finite size effect that vanishes for larger system sizes. The 158 boundary conditions introduce a measurable effect on the spontaneous nucleation rate 159 of retinal waves for smaller system sizes. This effect is due to correlations of distant 160 neighbors in the system. To gain a better understanding of the magnitude of correlation 161 on voltage fluctuations, we simulated a one-dimensional chain of neurons. Because we 162 are interested in the sub-threshold voltage fluctuations we simplified the neuron model 163 to a version that does not generate spikes, i.e. with a linearized deterministic part of the 164 dynamics (GJ and noise current unchanged) at the stable fixed point: Simulation results for the Pearson correlation coefficient of the membrane potential as 166 function of distance for a ring of 15 coupled neurons are illustrated in Fig. B(a). With 167 increasing coupling strength, there is a non-zero correlation between the voltages, even 168 for neurons as far as 4 space units apart. This leads to a measurable increase of the 169 overall membrane fluctuations of a ring with chain length up to 5 neurons for G = 0.5, 170 cf. Fig. B(b). For larger chain lengths, the periodicity has no effect on the voltage