Model for the single retinal ganglion cell

We use the phenomenological Izhikevich neuron model [31, 32], known for displaying biologically plausible dynamics. Due to its discontinuous fire and reset mechanism, it is a computationally efficient model of a bursting neuron. Comparable dynamics can be obtained from two-dimensional excitable models such at the Morris-Lecar model, under incorporation of an additional third dimension, e.g. a calcium-dependent potassium current, cf. Sec. 5.2 in [33]

The model can be regarded as a quadratic integrate-and-fire neuron for the membrane voltage V_i(t) of the ith neuron with an additional slow recovery variable u_i(t), also referred to as gating variable (cf. Fig 1(a) for the nullclines of the system): (1) (2) (3) The membrane recovery variable provides negative feedback to the voltage (cf. Fig 1(b) and 1(c) top). The parameters a, b, d as well as V_rest, V_crit, V_reset, and V_peak determine the spiking regime of the neuron, with V_rest < V_crit < V_peak. The time-scales of the voltage and gating variable are defined by τ_V and τ_u, respectively. For u(t) ≡ 0 and I(t) ≡ 0, V_rest and V_crit are the stable and the unstable fixed points of the dynamics, respectively. If V_i ≥ V_peak, the membrane potential is reset to V_reset, the kth spike time, t_i,k, is registered, and the recovery variable is increased by the constant value d. We choose parameters such that the burst characteristics of our model neuron illustrated in Fig 1 roughly agree with experimental measurements from Syed et al. [11]. Specifically, we aim at a burst duration of about 1 − 2 seconds (cf. Fig 1(c) bottom) and a spike frequency during bursts of about 5 − 15 Hz. We find those characteristics reasonably met for: a = 0.1, b = 0.3, d = 1.2, τ_V = 100 msec, τ_u = 0.0003⁻¹ msec, V_rest = −76 mV, V_crit = −48 mV, V_peak = 30 mV, V_reset = −50 mV. The bursting mechanism is illustrated in Fig 1. The chosen time-scale of the gating variable u is comparatively large, but not uncommon for cortical neurons [34].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Burst mechanism of the single neuron model.

(a) shows the nullclines of the Izhikevich neuron model in phase space (V, u) without current, RI = 0. The green dashed line shows the voltage nullcline and the blue dashed line shows the gating variable nullcline, respectively. Intersections of these two lines are fixed points of the system. The lower fixed point, indicated in red, is stable and represents the resting state of the neuron at (V, u) = (V_r, u_r) = (−64mV, −19.4mV). The gray vertical lines indicate the peak voltage V_peak and the reset voltage V_reset. (b) shows the path in phase space of a neuron that is initially in the resting position, but exposed to an external current with RI = 2 mV from t = 0. The temporal evolution of the separate components u and V is illustrated in (c).

https://doi.org/10.1371/journal.pcbi.1006355.g001

The total current RI_i = R[I_gap,i + I_noise,i] is a superposition of the intrinsic noise current and GJ currents from neighboring cells (see below). The intrinsic noise originates from fluctuations of the various channel populations (sodium, calcium, and different potassium channels, see e.g. [35]) and is approximated by white Gaussian noise: (4) with 〈ξ_i(t)〉 = 0 and 〈ξ_i(t)ξ_j(t′)〉 = δ_ij δ(t − t′) and D is the noise intensity. We perform simulations at discrete times with a time step of Δt = 0.1 msec according to an Euler-Maruyama integration scheme, see supporting information S1 Text.

Retinal network

Ganglion cells are distributed within the ganglion cell layer with a decreasing density towards the outer regions of the retina. For instance, the density in rabbits covers a range from 5000 cells/mm² down to 200 cells/mm² (the mean value is 800) [36]. In a previous study of retinal waves observed in rats, Butts et al. [4] used a ganglion cell density of ∼ 4000 cells/mm². In their simulations they placed neurons in a regular triangular lattice for which the given density translates to a lattice spacing of 17 μm. Because we focus on the rabbit retina, we assume a triangular lattice with a different lattice spacing of 38 μm, reflecting the lower cell density (800 cells/mm²) for this system. The reported experimental observations on characteristics of stage I retinal wave were obtained from retina patches of roughly 3 × 5 mm. A mean cell density of 800 cells/mm² translates to a total cell number estimate of 12,000 cells in the studied system. For comparability, we use a similar number of cells for simulations (i.e. 12,100 = 110 × 110). The triangular lattice structure can be seen in Fig 2(a).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Wave propagation in the deterministic system (D = 0).

Lattice structure of the network shown in (a) for the one-dimensional (1D) and two-dimensional (2D) simulations. Voltage traces of five model neurons (vertically shifted for better visibility), coupled in a one-dimensional chain with G = 0.1 (b) and G = 0.5 (c). The respective first neuron (bottom trace) was initialized in the bursting regime, i.e. (u(t = 0), V(t = 0)) = (u_rest, V_reset). Snapshots of waves on a two-dimensional triangular lattice (voltage and recovery variable in top and bottom panels, respectively) with G = 0.1 at different time instances as indicated (d).

https://doi.org/10.1371/journal.pcbi.1006355.g002

Here, we ignore for simplicity the inhomogeneous and irregular structure of the ganglion cell layer. We place N = n × n single ganglion cells in a rectangular domain on a triangular lattice such that every cell is connected with GJs to six nearest neighbors, (the lattice structure is illustrated in Fig 2(a)). For illustrative purposes, we will also consider a one-dimensional chain, in which each neuron has only two neighbors. Because we are interested only in stage I waves, prior to synaptogenesis, these cells are not connected to any other cells, i.e. bipolar and amacrine cells are not part of our model. We choose a common approach (e.g. [21]) to model the GJ current as diffusive and instantaneous coupling by (5) where G is the rescaled dimensionless GJ coupling, i.e. G = R/R_gap. The membrane resistance R of retinal ganglion cells can experimentally be measured and is in the range of 100-500 MΩ, e.g. [37]. R_gap is the GJ resistance between neighboring ganglion cells in the retina, which depends on the connexin type and the transjunctional voltage difference and is roughly R_gap ≈ 1GΩ [38, 39]. The values of R and R_gap imply a physiological range for our parameter of G ∈ [0.1, 0.5]. Because the time course of the action potential produced by our neuron model is only a coarse approximation of the electrophysiological shape of a spike, the GJ coupling may be stronger or weaker than assumed here. This gives additional justification for choosing a wider range of G.

For the two-dimensional setup, we apply two different boundary conditions. For estimating the noise dependence of propagation velocities and nucleation rates, we perform small system simulations (N∼50-260) with periodic boundary conditions in both directions (system on a torus) in order to avoid strong finite-size effects. Simulations of the full system with N∼12,000 are carried out with two additional layers of neurons on the boundary, that are not exposed to intrinsic noise (cells on the system boundary have fewer neighbors, between 2 and 5 instead of 6). Neurons in the two outer layers of the large simulations are discarded from all statistical evaluations.

Single propagating waves running through the network can be captured by the population activity [40] (6) where the index k runs over the spike times of the ith neuron. Hence, A(t) is the firing rate, averaged over the network and the time bin Δt_A. We use Δt_A = 0.5 seconds, which is comparatively large and covers multiple spikes when the cells are bursting.

Wave propagation

If we couple cells in a chain (as indicated in Fig 2(a)—1D) and initiate a burst in one of them, we see a propagation of the burst along the chain (cf. Fig 2(b)); similar voltage traces have also been seen in simulation of computational models of cortex slices, e.g. [41]. A higher propagation speed can be achieved by increasing the GJ conductance parameter G Fig 2(c). The picture is similar in our two-dimensional setup, for which snapshots are shown in Fig 2(d). In this case, the wave has been evoked by enforcing a burst in the lower left corner. It propagates as a circularly shaped wave front, which is a consequence of the regularity and rotational symmetry of the system. The gating variable u (lower row in Fig 2(d)) can be associated with the experimentally accessible calcium dynamics and resembles calcium fluorescences images [11]. Compared to the membrane potential (top row), the wavefront of the gating variable lags behind, as it slowly builds up during the burst.

In both, one-dimensional and two-dimensional simulations in Fig 2, we have set the intrinsic noise intensity to zero in order to illustrate that wave propagation does not hinge on the presence of fluctuations. We note already here, that the propagation speed in the two-dimensional system matches the order of magnitude of biologically observed values. To determine the speed of the waves from simulation such as shown in Fig 2(c), we approximate the wave’s shape as circular with a fixed center. We define a wavefront as the group of neurons that spike within the same time bin of Δt = 0.1 seconds (see left illustration in Fig 3(a)) and measure the front’s mean distance from the center and its mean time instance of occurrence. From the differences of these distances and times, we determine the mean velocity, which we find to be weakly distance dependent, but saturating at about 350 μm from the origin of the wave, cf. Fig 3(b). In the following, all velocity values are averaged over measurements for the range of distances 350 − 650 μm (shaded area in Fig 3(b)) from the point of initiation and we refer to this measuring method as concentric method. The velocities are shown in Fig 3(c) as a function of the GJ parameter for the physiologically relevant range of G (see Methods). We obtain velocities that are in the range of values observed in the rabbit retina [11], cf. the shaded area in Fig 3(b). The experimental mean value of about 450 μm/sec is attained for G ≈ 0.4.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Wave speed: Measurement and dependence on GJ coupling.

Neural groups with simultaneous burst onset of an exemplary simulation (time resolution Δt = 0.1 seconds) are shown in the panel (a) left, for three consecutive time bins in different colors. At large distances from the origin, the shape of a wavefront can be approximated as planar, cf. (a) middle. The mechanism of burst propagation can then be mimicked by a one-dimensional situation. Therefore in our theoretical derivations, the distance and coupling strength has to be modified, cf. (a) right and details in the main text. Squares in (b) represent the speed of the concentric wave (G = 0.4) as a function of the distance from the wave’s origin (lower left corner of the simulation domain), measured as described in the text. Alternatively, the speed can be assessed by measuring burst onset times along different fixed directions of the network, i.e. at blue and green sites shown in the inset of (b). The resulting wave speeds as functions of distance (blue and green lines) agree closely with the concentric method (squares in (b)). The speed shown in (c) is the mean value of simulation data (symbols) of the shaded area in (a) as a function of the GJ coupling G. Simulation results are compared to v_2D(G), Eq (9). The vertical and horizontal shadings indicate the physiological range of G (see Methods) and the observed velocities in the rabbit retina [11], respectively.

https://doi.org/10.1371/journal.pcbi.1006355.g003

The propagation and its speed can be theoretically understood as follows. Assuming a steep wave profile, the speed of the wave is given by the inverse of the time it takes a bursting neuron to excite its neighbors, times the displacement of the corresponding wave fronts. We refer to this time as burst onset time difference (BOTD). For simplicity, we neglect noise and consider in the following a one-dimensional setup consisting of three neurons: one initially quiescent neuron (i) is connected to a bursting neuron (i − 1) on one side and to a quiescent neuron (i + 1) on the other side. They are separated by the lattice spacing ℓ = 38 μm, hence the velocity is defined as v_1D = ℓ/T_B. Therein, T_B denotes the analytical approximation of the BOTD for this one-dimensional case.

The approximation T_B for the BOTD between neighboring neurons can be derived using three assumptions (details in S1 Text). First, we assume a constant gating variable (u(t) ≈ u_r = const), which is reasonable on a short time-scale, because τ_u ≫ τ_V. Second, we replace the voltage variable of the bursting neuron V_i−1(t) by its temporal average , that can be analytically calculated (see S1 Text) and for our standard parameters is mV. Third, we replace the voltage of the quiescent neuron that is not directly connected to the bursting neuron by the resting potential, V_i+1 = V_r. Consequently, the GJ current seen by the driven neuron reads , and the resulting dynamics until the voltage V_i reaches the peak potential for the first time is effectively one-dimensional and can be recast to the form (cf. details in S1 Text): (7) This first order ordinary differential equation can be solved via separation of variables to find t(V). We obtain it by first calculating the difference of the times from the voltage being at its peak potential and its resting potential. However, the driven neuron is already exposed to the driving GJ current while the voltage of the bursting neuron travels to its first spike time (cf. Fig. A of S1 Text). Therefore, for simplicity we subtract the first inter-spike interval T_ISI from the beforehand calculated time difference: (8) The explicit expression is lengthy and derived in S1 Text, resulting in Eq. O of S1 Text. Comparing T_B to simulations of a one-dimensional chain shows a reasonable agreement (cf. Fig. A of S1 Text), although the theory overestimates the simulated values, in particular, for larger values of G. For comparison we also discuss a corresponding result for the wave velocity in the continuum limit in S1 Text.

In the two-dimensional setup at larger times, the wave attains a planar shape as indicated in Fig 3(a), where red circles represent bursting neurons and blue and yellow circles represent driven and quiescent neurons. Now, we assume that the wave front is perfectly flat and all neurons shown in the same color share an identical voltage. In that case, the propagation mechanism simplifies to two bursting neurons exciting one quiescent neuron, whose membrane potential is further affected by two quiescent neurons. Hence, we can mimic the quasi one-dimensional situation by doubling the value of G and additionally taking into account the modification of the effective length, i.e. ℓ_eff = (3/4)^1/2 ℓ, see Fig 3(a). Consequently, we can approximate the velocity in the two-dimensional system as (9) Calculated velocities v_2D(G) are shown in Fig 3(c) by the blue line, underestimating the true velocity (circles) but providing a correct order-of-magnitude estimate. Note that so far we restricted the considerations to a purely deterministic setup. Our simulations with noise indicate that moderate fluctuations have only little impact on the mean velocities.

Wave nucleation

In the stochastic version of our system, we observe spontaneous waves that resemble those found in experiments [11]. Experimentally, it was observed by Syed et al. [11] that the spontaneously nucleated waves appear with a mean inter-wave interval T_IWI of 36 seconds. In our model, waves are initiated by noise, since neurons are set in the excitable regime and cannot generate periodic spiking or bursting without external input. We expect that the nucleation rate per neuron depends strongly on the noise intensity D. To characterize this dependence, we simulate small systems (N∼50-260, see Methods) with periodic boundary conditions for two different values of GJ coupling and different noise intensities, cf. Fig 4.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Arrhenius plot.

Spontaneous nucleation rate as function of the inverse noise intensity obtained from four two-dimensional systems with different system sizes as indicated and periodic boundary conditions. From the linear fit of these data, an effective potential barrier ΔU and a rate prefactor r₀ can be estimated (dependence of ΔU on system size shown in inset).

https://doi.org/10.1371/journal.pcbi.1006355.g004

With the understanding that every neuron has the same chance to trigger a wave, the global nucleation rate should be linear with N to a first approximation. Thus we measure the nucleation rate per neuron as r = 1/(T_IWI N). As demonstrated in Fig 4 by the linear dependence of the rate’s logarithm on the inverse noise intensity, we obtain an Arrhenius rate (10) The effective potential barrier ΔU depends on G and the system size N and saturates for sufficiently large systems (inset) for both values of G.

The increase of the potential barrier with G can be understood to first approximation by the effective change of the current-voltage relation in the single neuron. The GJ coupling term Eq (5) leads to an effective increase in the leak current that stabilizes the resting potential and makes it harder to initiate a burst. This mechanism is dominant in comparison to the influence of other coupling effects and the stochasticity of the neighbors on the nucleation rate (supported by additional simulations, see Fig. B of S1 Text).

The more subtle dependence of ΔU on the system size can be explained as follows: Coupling stochastic neurons in small systems with periodic boundary conditions leads to spatial correlations and thus effectively to stronger noise. This effect can be neglected for large system sizes or weak coupling, but has a measurable effect otherwise (cf. Fig 4 and Fig 4 inset).

Discussion of large-scale simulation results

Our results so far can be used to predict the mean inter-wave interval and the propagation speed of retinal waves for a system size N = 12,100 that roughly corresponds to the experimentally studied patch size in Ref. [11]. Vice versa, we can infer an approximate value of the noise intensity D that leads to the experimentally observed value of T_IWI = 36 seconds and test this in numerical simulations of the full system.

For our estimation of the rough value of the noise intensity in a large system, we have to take into account that the single neuron undergoes a substantial refractory period of T_ref ≈ 14 seconds after bursting (estimated from small-system simulations investigating the minimal mean inter-wave interval for various noise intensities). The mean inter-wave interval is then given by T_IWI = T_ref + 1/[N ⋅ r(D)] and the estimated value of the noise intensity follows from the Arrhenius law, Eq (10), as (11) (for G = 0.4, and r₀ = 6 and ΔU = 0.71, fit parameters from Fig 4, solid line with N = 256).

The estimated parameters, G = 0.4 and D = 0.050, can now be used in a large-scale simulation. In Fig 5(a), we show snapshots of the full system’s gating variable (a proxy for the experimentally accessible calcium concentration). The wave front seen in the experimentally observable area (box in Fig 5(a)) looks similar to experimental measurements, cf. Ref. [11]. From Fig 5(b), it becomes evident that the mean inter-wave interval becomes much shorter for a slight increase in D. The mean inter-wave interval at these parameter values is not exactly 36 seconds, but somewhat larger: these statistics depend very sensitively on the value of the noise intensity (i.e. on the second leading digit, cf. Fig 5(c) middle). This is seen in the global population activity, that reveals a wave going through the system as a single peak vs. time.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Large-scale simulations.

A network of 12,100 GJ coupled and noisy Izhikevich neurons display spontaneously nucleated waves that propagate with velocities that are comparable to experimentally observed values. (a) Snapshots of the gating variable (associated to a proxy for the calcium concentration) at different time instances during one wave running through the system (G = 0.4, D = 0.05). The small rectangle indicates the dimensions of the experimentally accessible observation area [11]. (b) Population activity A(t) (with Δt_A = 0.5 seconds, see Eq (6)) of the entire system over a larger time window for different noise levels. One wave, as shown in (a) collapses here into a single peak; time differences between adjacent peaks are the inter-wave intervals T_IWI,i (one indicated by an arrow). (c) Mean velocity, mean inter-wave interval and standard deviation of the subthreshold membrane voltage as a function of the noise intensity for a small range around the estimated value D_target = 0.05. Dashed lines indicate experimental mean values from Ref. [11], solid red line shows the wave speed for D = 0, extracted from the circle at G = 0.4 in Fig 3.

https://doi.org/10.1371/journal.pcbi.1006355.g005

The dependence of crucial neural statistics on the noise intensity is illustrated in Fig 5(c). In contrast to the mean inter-wave interval, the mean velocity of the wave does not depend strongly on the noise (Fig 5(c), top) but stays close to the experimentally observed mean value (dashed line). This is due to the fact, that the wave, once it is initiated, is largely determined by the deterministic propagation mechanism explained above. The fine tuning of the noise intensity shows that the experimental value of 〈T_IWI,exp〉 = 36 seconds is attained for a noise level of D = 0.052, slightly larger than D* (estimated above). How realistic is this noise level? To address this question, we show at the bottom of Fig 5(c) the standard deviation of the subthreshold voltage fluctuations, σ_V, as a function of the noise intensity D. σ_V increases only slightly with D and attains values around 1.6 mV.

To our knowledge, there are no detailed investigations of intrinsic noise sources in retinal ganglion cells at embryonic age. Because in this developmental stage there are no chemical synapses present [42], the synaptic background fluctuations can be excluded for our system: in the recurrent networks of the cortex, fluctuations stem mainly from the many synaptic interactions among the neurons, resulting in the famous asynchronous irregular state [43] that can be highly variable [44–46]. In our system, one likely source of variability is channel noise that typically leads to small membrane potential fluctuations with a standard deviation σ_V below 0.6 mV [47, 48]. The noise intensity that is required for the experimentally observed inter-wave interval results in sub-threshold voltage fluctuations that are three times bigger, cf. Fig 5(c) bottom, suggesting that besides ion channel noise there are additional sources of fluctuations present. These could result from stochasticity of GJs itself but also indirectly from GJs via differences in individual resting potentials (for the heterogeneity of the resting potential in similarly sized cells, pyramidal cells in the cortex, see [49]). In any case, the apparent voltage fluctuations of about 1.6 mV are well within the range of experimentally observed voltage noise in embryonic ganglion cells (cf. Fig. 1 in Ref. [11]).

Summary

The investigations presented in this paper propose a GJ-based model of stage I waves in the developing retina. Starting with a neuron model that roughly reproduces the spiking properties of a burst of one single retinal ganglion cell, we incorporated GJ coupling of physiologically plausible strength and temporally uncorrelated fluctuations. This allowed us to reproduce the characteristics of wave nucleation and slow wave propagation in the early retina. Earlier it was believed that GJs can play a role in fast neural transmissions only [5, 9], since the current in electrical synapses responds much quicker than neurotransmitters in chemical synapses. As shown in our paper, however, it is possible to obtain a limited transmission speed in a simple Ohmic model of the GJ coupling. Furthermore, although stochastic fluctuations are strong enough to ignite bursts with the correct nucleation rate, they do not distort the propagating fronts very much, i.e. the wave propagation is still a reliable process.

The reason for the slow transmission we observe can be found in the nonlinear dynamics of the single neuron. The Izhikevich model that we use for the ganglion cell is essentially a quadratic integrate-and-fire neuron model with a slow adaptation variable. This model is the normal form of a saddle-node bifurcation and has a pronounced latency if close to this bifurcation, i.e. the spike response to a current step (in our case provided by a neighboring bursting cell) is considerably delayed because the system experiences the “ghost of the former fixed point”, see Ref. [50]. The presence of weak noise modifies this picture only slightly [51].

Although our model accounts for the most important features of wave nucleation and propagation for stage I retinal waves, it cannot explain the strong variability of the experimentally measured statistics (error of velocity ±91 μm/sec [11]). This is due to a number of model simplifications, which we now concludingly discuss. Firstly, the real system is much more heterogeneous than in our model, both with respect to the lattice structure as well as with respect to the local coupling between cells; secondly, GJs may couple more than next neighbors and their conductivity may be noisy and voltage gated; thirdly, the detailed dynamics of ganglion cells is certainly more complex than can be captured by the Izhikevich model; last but not least, the white Gaussian noise in our model is a rather coarse approximation of the channel noise and other fluctuations in the system.

In our model, we arranged the neurons on a highly regular lattice with a cellular spacing according to an experimentally determined mean value of cell density, neglecting the strong heterogeneities in the distribution [36]. On this lattice, each cell is connected to exactly six nearest neighbors. Given the aforementioned heterogeneity, the numbers and distances between neighbors will be more broadly distributed than in our model. Incorporating these heterogeneous features in the simulations would likely broaden the range of observed velocities and thus better reflect the considerable variability found in experimentally measured values.

The soma size of (rabbit) retinal ganglion cells (< 30μm, e.g. Ref [36]) is smaller than our employed lattice spacing, implying GJ coupling between dendrites rather than soma-soma coupling only. The size of the dendritic arbor of retinal ganglion cells is ∼ 100 − 130μm, thus suggesting direct communication between cells that are up to the threefold of the lattice spacing apart. In our simulations with only next-neighbor coupling, we could reproduce the experimentally observed velocity with a comparatively large coupling constant of G = 0.4 (physiological range was G ∈ [0.1, 0.5], see Methods). It is conceivable, that this large G value is an effective description of a system with larger effective GJ neighborhood but with a smaller (and possibly distance-dependent) coupling value G. Put differently, we expect similar results for the wave speed in a system with extended coupling neighborhood but reduced coupling strength per connection (with the latter still being within the physiological range).

Regarding the neuron model and the incorporation of noise, we note that for developed retinal ganglion cells detailed multi-compartment conductance-based models with stochastic ion channels exist [35]. With more electrophysiological data available, it will certainly be possible to develop biophysically more realistic models of the bursting ganglion cell at the early stage. Furthermore important for our problem will be the incorporation of stochastic models of GJs [52] with voltage-dependent kinetics [53, 54] and the heterogeneity of physiological parameters such as the resting potential. Such detailed models are certainly difficult to simulate for large networks but could be employed to estimate the total noise intensity in the system and to identify the dominant noise source, cf. similar approaches in Refs. [35, 55, 56].