Abstract
A novel, biophysically realistic model for earlystage, acetylcholinemediated retinal waves is presented. In this model, neural excitability is regulated through a slow afterhyperpolarization (sAHP) operating on two different temporal scales. As a result, the simulated network exhibits competition between a desynchronizing effect of spontaneous, cellintrinsic bursts, and the synchronizing effect of synaptic transmission during retinal waves. Cellintrinsic bursts decouple the retinal network through activation of the sAHP current, and we show that the network is capable of operating at a transition point between purely local and global functional connectedness, which corresponds to a percolation phase transition. Multielectrode array recordings show that, at this point, the properties of retinal waves are reliably predicted by the model. These results indicate that early spontaneous activity in the developing retina is regulated according to a very specific principle, which maximizes randomness and variability in the resulting activity patterns.
Introduction
Retinal waves are a well studied example of spontaneous, patterned activity during CNS development and have been characterized in many vertebrate species (Wong, 1999). They appear before the retina is capable of responding to light and are observed as rhythmic bursts of action potentials in retinal ganglion cells (RGCs) that propagate in wavelike events across the surface of the retina. An involvement of retinal waves in visual system development has long been hypothesized, and recent studies demonstrate that manipulating wave activity can indeed affect the maturation of the projections of RGC axons (Firth et al., 2005) and the development of the retinal circuitry (Sernagor and Grzywacz, 1996; Sernagor et al., 2001; Mehta and Sernagor, 2006). It is also increasingly clear that in some cases specific properties of retinal waves are required for normal development (MuirRobinson et al., 2002; McLaughlin et al., 2003; Cang et al., 2008; Shah and Crair, 2008; Sun et al., 2008). This suggests that the adequate regulation of spontaneous retinal activity is an important prerequisite for normal visual system development.
Important questions remain unresolved in terms of understanding the physiological processes underlying the typical spatiotemporal properties of retinal waves and how they are maintained and regulated during development. Although some wave features, such as their frequency of occurrence or propagation velocity, are highly variable across species, there are also important commonalities. In particular, in all species investigated so far, retinal waves display a high degree of randomness. Successive retinal waves occur at variable intervals, and exhibit a high degree of variability and randomness in terms of initiation point, size and trajectory (Feller et al., 1997; Stellwagen et al., 1999). Similar forms of spontaneous activity have also been found in other developing neural structures including the hippocampus and spinal cord (Feller, 1999; O'Donovan, 1999).
Here, by combining biophysical modeling and multielectrode array (MEA) recordings, we show that retinal waves arise at a very specific network state which can be interpreted in terms of a percolation model (Essam, 1980), a classical model in statistical mechanics. In particular, we demonstrate that a slow afterhyperpolarization (sAHP) current that has recently been characterized in starburst amacrine cells (SACs) (Zheng et al., 2006) acts on two different time scales, providing a mechanism to desynchronize neurons and dynamically regulate the network connectedness. A mean field approximation shows that this network can undergo a percolation phase transition, separating the states of purely local and global functional connectedness. A comparison of simulated activity with experimental MEA recordings shows that retinal waves emerge close to this transition point, which is manifested as powerlaw distributions of wave sizes and durations.
Materials and Methods
Computational model.
The model implements a regular, twodimensional, hexagonal lattice of starburst amacrine cells. Starburst amacrine cells were modeled using the membrane equation, with a membrane capacity C = 160 pF, resistance R = 180 MΩ and resting potential V_{r} = −65 mV. Synaptic inputs and membrane channels have conductances g_{i}(t) and associated reversal potentials E_{r}^{i}. A fast voltagedependent calcium current was modeled using an adiabatic approximation of its fast dynamics (Morris and Lecar, 1981): The halfactivation point was set to V_{Hca} = −20 mV, the slope V_{SCa} = 10 mV, the peak conductance g_{m}^{Ca} = 10 nS, and the reversal potential E_{r}^{Ca} = 50 mV. An sAHP current was modeled as a calciumdependent potassium conductance with four gating particles r(t), similar to I_{sAHP} in cortical neurons (Abel et al., 2004). This current depends on the intracellular calcium concentration [Ca(t)], and a second, slower calciumdependent process s(t): The peak conductance was g_{m}^{sAHP} = 30 nS, the reversal potential E_{r}^{sAHP} = −90 mV, the activation constant α = 2400, and the decay time τ_{r} = 5 s. A cooperative activation of the slow calciumdependent process was assumed, with an activation constant β = 30, Michaelis constant γ = 10^{−8} and decay time τ_{s} = 50 s. The calcium concentration was calculated as follows: with scaling factor δ = 10^{8} and decay time τ_{Ca} = 50 ms. Neurons have an intrinsic source of noise, which was implemented as filtered shot noise: where P(t,r(V)) is a Poisson process with a voltagedependent rate r(V), event size α = 200 pS and reversal potential E_{r}^{n} = 50 mV. The source of cellintrinsic noise in SACs is not known at present, and may be the result of stochastic ion channel gating, so a voltagedependent ratemodulated Poisson process was used with λ(t) = k_{0}p(V)(1 − p(V)), where p(V) was the calcium channel activation function in Equation 1, and k_{0} = 1400/s. These activation variables were chosen as simulations suggest that low calcium channel densities in small, distal SAC dendrites can induce spontaneous firing in SACs (M. H. Hennig, unpublished data). Simulations have, however, also shown that the same results are obtained with other noise models, such as an unmodulated Poisson process (data not shown).
SACs receive excitatory input from all neighboring cells j within in a circular region of a radius of three cells with a Gaussian distancedependent synaptic weight with SD σ_{w} = 1. The threshold for transmitter release was V_{T} = −60 mV, and synaptic currents were lowpass filtered with a time constant τ_{syn} = 0.2 s to account for the slow cholinergic transmission in immature SACs (Zheng et al., 2004). The total synaptic conductance is then g_{syn} = g_{syn}^{max}Σ_{j}g_{syn}^{j} with g_{syn}^{max} = 320 nS. A forward Euler scheme was used for integration, with a time step of 1 ms. Unless stated otherwise, statistical measures were obtained from 4000 s simulated activity on a lattice of 56 × 56 SACs, where the central 50 × 50 neurons were used for analysis to reduce boundary effects. The simulation software was implemented in C++. The source code is available upon request.
MEA recordings.
This study was done using C57b1/6 neonatal mice and turtle embryos from the species Pseudemys Scripta Elegans. All animal procedures were conducted under the United Kingdom Home Office, Animals (Scientific procedures) Act 1986. Mouse pups were killed by cervical dislocation and enucleated before retinal isolation. Turtle embryonic ages were determined according to specific staging criteria (Yntema, 1968). The embryos were anesthetized by hypothermia, decapitated and enucleated before retinal isolation. The isolated retina was then transferred to the experimental chamber and placed, RGC layer facing down, onto MEAs consisting of 60 titanium nitride electrodes (30 μm diameter, 200 μm spacing) arranged in an 8 × 8 grid on indium tin oxide substrate (Multi Channel Systems). Better coupling between the tissue and the electrodes was achieved by placing a small piece of polyester membrane filter (5 μm pores) (Sterlitech) on the retina followed by a slice anchor holder (Warner Instruments). For mice, the retina was kept at 32°C and continuously perfused (2–5 ml/min) with artificial CSF containing the following (in mm):118 NaCl, 25 NaHCO_{3}, 1 NaH_{2}PO_{4}, 3 KCl, 1 MgCl_{2}, 2 CaCl_{2}, and 10 glucose, equilibrated with 95% O_{2} and 5% CO_{2}. For turtles, the retina was kept at 29°C and perfused with Ringer's solution containing the following (in mm): 96.5 NaCl, 2.6 KCl, 2 MgCl_{2}, 31.5 NaHCO_{3}, 10 glucose, 10 HEPES and 4 CaCl_{2}, equilibrated with 95% O_{2} and 5% CO_{2}. Signals were amplified (gain 1200) and acquired using a 128channel analog to digital converter (Multi Channel Systems MC_Card). Signals were digitized at 25 kHz and acquired without filtering using the software MC_Rack (Multi Channel Systems). The time of occurrence of spontaneous spikes was thresholddetected with MC_Rack (the typical threshold was at signal amplitude that is 3× below the baseline noise). Hence, the firing rate on each electrode reflects the overall activity level generated by all RGCs on that same electrode (typically 4–5 cells). Using the software MC_Data Tools (Multi Channel Systems), the times of spike occurrence were converted into text files for further analysis. Retinal waves were then detected using custom software written in Matlab (Version 7.21; Mathworks).
Analysis of MEA recordings.
Spike trains of MEA recordings of retinal waves typically show occasional bursts of spikes at high rates, which are embedded in lowfrequency background activity. For a precise and reliable burst detection, first the ranks R(t) for the interspike intervals and the probability distribution P(C_{s}) for the spike count in a fixed time window were obtained (mouse, 1 s; turtle, 4–6 s; these values reflect the typical burst duration). Then, each spike train was analyzed spike by spike, and bursts were detected using a combination of two threshold criteria based on the mean firing rate on each electrode. Specifically, the onset of a burst was defined as an event where R(t) < θ_{r}, where the rank threshold θ_{r} was manually adjusted (typical value for mouse data were θ_{r} = 0.2 and for turtle data θ_{r} = 0.1), and where the spike count C_{s} in the fixed window beginning at the current spike exceeded a threshold. The threshold was taken as the value of C_{s} where P(C_{s}) = θ_{s} (mouse, θ_{s} = 0.05; turtle, θ_{s} = 0.3), and calculated individually for each electrode. The end of a burst was determined as the time of the spike where the spike count first dropped below half of the onset threshold. This method was insensitive to the large differences in average and peak firing rates between different electrodes.
Waves were detected as temporally overlapping groups of bursts. The neighborhood relationship between electrodes was not taken into account, as occasionally multiple waves can be initiated at the same time on the area covered by the MEA (compare Fig. 4A). To avoid the overlap of temporally separated waves, a maximum duration was imposed on bursts (2.5 s or longer). This value was chosen because activity propagation between neighboring electrodes (with a distance of 200 μm, or 280 μm for the diagonal) is expected to take 2 s or less, since typical wave propagation velocities are >100 μm/s (Wong et al., 1993; Feller et al., 1997). Visual inspection confirmed that this method was suitable for detection of coherent wave events, while shorter maximal burst durations could lead to splitting of waves. Not imposing a maximal duration had only minor effects for the mouse data, but was necessary for the turtle recordings, where bursts often lasted 10 s or longer.
Estimation of powerlaw exponents.
To obtain reliable estimates of the parameters of the powerlaw distributions of wave sizes and lifetimes, and to test the quality of these fits, we used the methods outlined in detail by Clauset et al. (2007). In the following, we provide a brief summary of the way these techniques were applied here. A powerlaw probability distribution has the form with x ≥ x_{m} > 0 for all x and α > 1. x_{m} is the lower bound for the powerlaw distribution and C the normalization constant. For realvalued x, we have C = (α − 1)x_{m}^{α − 1}, and if x is discrete valued, C = ζ (α, x_{m})^{−1}, where ζ (α, x_{m}) is the generalized zeta function.
To estimate the powerlaw exponents for the distributions of retinal wave sizes and lifetimes, we used a combination of minimization of the Kolmogorov–Smirnov (KS) statistic to estimate the lower bound x_{m} and a maximum likelihood estimator (MLE) to estimate the exponent. For the discretevalued wavesize distributions with n data points, α was estimated by maximizing the loglikelihood function: For the continuous lifetime distributions, a closedform estimator for the exponent β is given by the following: Since the lower end of the distributions of the simulated and experimental data are affected by aliasing, we estimated powerlaw exponents and calculated the KS statistic for different values of x_{m} up to a prespecified limit x_{m}^{+}. The value of x_{m} which minimized the KS statistic, and the corresponding exponent, were then taken as the best estimates. The fixed upper limits for the experimental data were x_{m}^{+} = 6 for the wavesize distributions and x_{m}^{+} = 5 s for the lifetime distributions, to prevent spurious effects due to the small sample sizes.
The goodnessoffit was estimated by comparing the KS statistic for the real data to that obtained from synthetic data sampled from powerlaw distributions with the same parameters. To this end, 2000 synthetic data sets, each with the same number of data points as in the original set, and sampled from a powerlaw distribution with the same parameters, were generated, and their parameters were estimated. Then, the KS statistic was calculated for each data set, and the p value was determined as the fraction of times the KS statistic for the synthetic data sets exceeded that for the original data. For experimental data sets, the estimates for the exponents and x_{m} were obtained as specified above. To reduce the time needed to compute the p values for the discrete simulated wavesize distributions, which typically contained >1000 events, a continuous approximation for the MLE was used:
Results
Characterization of simulated SAC activity
Retinal waves early in development depend on cholinergic synaptic transmission between SACs (Feller et al., 1996; Zhou and Zhao, 2000), spontaneous, calciumchannel mediated bursts in SACs, and a sAHP conductance which regulates their intrinsic excitability (Zheng et al., 2006). These mechanisms have been implemented with a recurrent network of conductancebased model neurons, which reproduce the physiological properties of immature SACs.
An isolated simulated SAC produces infrequent, spontaneous bursts caused by cellintrinsic noise (Fig. 1A, black traces). Calcium influx during a burst activates the sAHP current, which induces a strong, longlasting hyperpolarization preventing further spontaneous bursts. Experimental manipulations have shown that longer or more intensive depolarizations increase the sAHP duration (Zheng et al., 2006). This was modeled by combining direct, calciumdependent current activation with an additional slow, presumably secondmessenger mediated process. The relaxation time of the current then depends on the amount of calcium influx and is longer after stronger or more intensive depolarization. This is illustrated in Figure 1A, which compares the effects of brief current injections with that of spontaneous bursts. In either case, calcium influx activates the sAHP current, which leads to a temporary dampening in cellular excitability and suppression of subsequent spontaneous bursts. Increasing the current intensity leads to a stronger activation of the slow component of the sAHP current, which in turn results in a longer suppression of spontaneous bursts. Therefore, the average latency of the first spontaneous burst after a current injection increases monotonously with the current magnitude, until it saturates at a value determined by the decay time of the slow sAHP component (Fig. 1B).
This differential activation of the sAHP current becomes important when SACs are synaptically coupled. Spontaneous bursts can then occur in isolation (Fig. 1C, asterisks), and coincident spontaneous bursts in nearby neurons trigger propagating waves of activity (Fig. 1C, diamonds). During a wave, strong synaptic drive from multiple, simultaneously active SACs leads to a stronger depolarization than intrinsic bursts, which in turn causes a longerlasting reduction in excitability of the participating SACs. Since waves uniformly reduce the excitability of participating neurons, wave boundaries are determined by recently active retinal areas, as shown in previous studies (Feller et al., 1996, 1997).
Intrinsic bursts also temporarily reduce the excitability of single SACs, which can prevent their participation in a subsequent wave (Fig. 1C, arrow). Since the sAHP current induced by intrinsic bursts decays faster than following waves, this results in competition between a desynchronizing effect of intrinsic bursting, and the synchronizing effect of synaptic transmission. This requires that synaptic connections between SACs are sufficiently weak, such that an isolated intrinsic burst cannot trigger a propagating wave (Fig. 1C, asterisks). Under these conditions, propagating activity patterns emerge, exhibiting a mixture of randomness and spatial coherence as in experimentally recorded retinal waves. Two examples of simulated waves are shown in Figure 2, A and B, illustrating that a wave only recruits neurons that are sufficiently depolarized to generate bursting activity, but stops when it encounters neurons with a reduced excitability due to recent activity (supplemental Movie, available at www.jneurosci.org as supplemental material).
Simulated retinal waves
An important prerequisite for the appearance of the characteristic random, noncyclic activity patterns of retinal waves is therefore the coexistence of synchronizing, propagating activity and desynchronizing, intrinsic bursting in the network. This condition is met for a range of physiologically plausible parameter combinations in the model. To determine when the typical properties of retinal waves are best reproduced, we investigated the properties of simulated waves for different intensities of the cellintrinsic noise (parameter k_{0} in Eq. 6; note that this parameter reflects the peak noise rate, and that the actual rate is typically much lower as it also depends on the membrane potential). Changing the noise intensity directly affects the rate of spontaneous, intrinsic bursting, and is hence best suited to determine its role in network activity desynchronization.
The simulations show that a gradual increase in noise intensity causes a transition between two qualitatively different types of spatiotemporal network activity (see also supplemental Fig. S1, available at www.jneurosci.org as supplemental material). For weak noise, the network produces a mixture of large and small waves at a low frequency, and the distributions of wave sizes and lifetimes have a bimodal form (Fig. 2C–E, blue curves). This reflects strongly correlated network activity, as evident in the slow decay of the spatial correlations (Fig. 2F, left), and occurs because SACs produce few intrinsic bursts that hinder the propagation of activity. The lowwave frequency is due to the low probability of simultaneous bursting of nearby neurons, which is required for wave initiation.
Increasing the noise intensity leads to the appearance of more frequent waves of intermediate size and duration (Fig. 2C–E, red curves). At a certain noise intensity (k_{0} ≈1.4 kHz; medium noise), the wavesize distribution changes from a bimodal to a monotonously decreasing function, and follows a power law with an exponent of α = 1.5. Rescaling the lattice of SACs by combining multiple neurons, and declaring each group of neurons as active only when more than half of the neurons are active (majority rule) yields the same powerlaw distributions (supplemental Fig. S2, available at www.jneurosci.org as supplemental material), confirming that the activity indeed has scalefree properties. At the same time, the wave lifetime distribution takes the form of an exponentially truncated power law with an exponent of β = 2, which is also preserved after rescaling. In this case, the stronger desynchronizing effect of the more frequent intrinsic bursting activity leads to a faster decay of the spatial correlations than at low noise intensities (Fig. 2F, middle).
At noise intensities above this point, waves become even more frequent, but the largest waves are no longer observed, and the wavesize and lifetime distributions now have the form of exponentially truncated powerlaw functions (Fig, 2C–E, green curves). The correlations now decay rapidly, consistent with strong desynchronization of the network activity through very frequent intrinsic bursting (Fig. 2F, right). These simulated waves strongly resemble the patchy, nonpropagating activity observed in more mature animals (Sernagor et al., 2003; Syed et al., 2004).
These observations show that this network undergoes a transition from a highly synchronized to a strongly desynchronized state, and produces activity patterns clearly different from oscillatory, phasedlocked behavior, where more compact event size and interval distributions are expected. The observed noncyclic behavior also manifest itself through the random, unbiased distribution of wave initiation points (Fig. 2G), a central feature of retinal waves (Feller et al., 1996). Yet, subsequent wave trajectories are not entirely random (Fig. 2H) but consistent with a historydependent regulation of the neural excitability (Feller et al., 1997).
Simulated of pharmacological manipulations
So far, the transition from the synchronized to the desynchronized network state has been characterized by varying the noise intensity, but the same behavior can also be observed by systematically changing other, experimentally more accessible, parameters. In particular, certain pharmacological manipulations can either increase or decrease synchronization, and hence induce transitions between the two regimens. To directly reproduce and interpret these effects with the model, we initially choose parameters such that powerlaw distributions are obtained to reproduce normal retinal waves (see below), and then simulated experimental manipulations by changing appropriate parameters. First, we investigate the effects of partially blocking cholinergic synaptic transmission, which has been shown to reduce wave frequency and spatial (Bansal et al., 2000; Sernagor et al., 2000, 2003). To simulate this, the peak synaptic conductance was uniformly reduced in the model. Consistent with the experimental results, this leads to a decrease in mean wave frequency and wave size as well as duration (Fig. 3A). These effects occur because more coincident bursts are required to trigger a wave when synapses are weak. This causes not only a decrease in wave frequency, but also a relative increase of the rate of intrinsic, desynchronizing bursts. As a result, the network enters a more desynchronized state, which leads to the disappearance of the largest waves (Fig. 3A, left). The opposite effect, an increase of wave size, duration and frequency, is observed when the synaptic strength is increased in the model. This mimics the experimentally observed increased frequency of calcium transients in the chick retina during blockade of ACh breakdown (Catsicas et al., 1998). However, a general elevation of the synaptic drive with nicotinic ACh receptor (nAChR) agonists abolishes wave activity (Zhou and Zhao, 2000), perhaps due to receptor desensitization, ion channel inactivation or permanent activation of potassium conductances.
Second, an increase of cAMP levels, which reduces the strength and affect the kinetics of the sAHP current in SACs (Zheng et al., 2006), has been shown to increase wave size and frequency (Stellwagen et al., 1999). This could be reproduced by reducing either the sAHP peak conductance (Fig. 3B), or the activation variable of the slow activation component of the sAHP current (Fig. 3C). Both manipulations cause the cellular excitability to be restored more quickly after activity, and lead to more frequent waves. Additionally, the desynchronizing effect of intrinsic bursts, mediated by the sAHP current, is reduced. This increase the relative frequency of the largest waves, and results in bimodal wavesize distributions as has been observed in calciumimaging experiments (Stellwagen et al., 1999). Enhancing the sAHP current or the slow activation component, however, causes the opposite effect, with a decrease in average wave size, duration and frequency. This manipulation has not been tested experimentally so far.
Collectively, these results show that qualitatively different manipulations of synaptic efficacy, cellular excitability or intrinsic noise intensity (compare Fig. 2C–E) all have similar effects on the synchronization of the activity of the network. Each of these variables can be used to induce transitions between the synchronized and desynchronized network state.
Analysis of MEA recordings
Our results raise the question which of these two network states can best reproduce retinal waves. Calciumimaging studies indicate that wavesize distributions are decreasing, highly skewed functions (Feller et al., 1997; Stellwagen et al., 1999; Bansal et al., 2000), similar to the powerlaw distributions predicted by the model at the transition point between the synchronized and desynchronized state. To test this experimentally, MEA recordings from mouse and turtle retinas were performed at early developmental stages, where waves primarily depend on cholinergic synaptic transmission (Feller et al., 1996; Sernagor and Grzywacz, 1999; Sernagor et al., 2003; Syed et al., 2004). If a behavior similar to that found in imaging experiments is observed with MEA recordings covering larger retinal areas, this would provide strong support for powerlaw distributions since only a powerlaw function retains its shape during spatial rescaling. However, due to the spacing between electrodes, MEA recordings also undersample the RGC mosaic, which will affect the measurement of wave sizes. Simulations show that spatial undersampling distorts powerlaw distributions in related models (Priesemann et al., 2007; V. Priesemann, M. H. J. Munk, and M. Wibral, unpublished observations), and increases the estimated powerlaw exponents in our simulations (supplemental Fig. S3, available at www.jneurosci.org as supplemental material). Hence, we would generally expect to find more rapidly decaying functions than predicted by the model.
Waves were detected as single or multiple coincident bursts of spikes on multiple electrodes (Fig. 4A) (see Materials and Methods). The normalized histograms in Figure 4, B and C, shows that the wavesize and lifetime distributions can, with one exception (Fig. 4E), be fitted by power functions (see Materials and Methods for details on the fitting procedure). The estimated exponents are close to, but always larger that the predicted α = 1.5 for wave sizes and β = 2 for lifetimes (ranges: α = 1.55–1.93; β = 2.32–2.84). These deviations are expected for a moderate undersampling of the RGC mosaics (supplemental Fig. S3, available at www.jneurosci.org as supplemental material). We also estimated the goodnessoffit by comparing the Kolmogorov–Smirnov statistic with that of synthetic data sets, which yields p values to quantify the support for a powerlaw hypothesis (see Materials and Methods). However, simulations show that spatial undersampling also affects these p values, in particular for the wavesize distributions (data not shown). Still, in one case, weak support for a powerlaw was found for the wavesize distribution (Fig. 4D, P2 mouse), and in four cases, this test supports the powerlaw hypothesis for the lifetime distributions (Fig. 4B–D,G). The analysis of two additional shorter dataset from postnatal day 5 (P5) mouse retinas produced similar results (data not illustrated; retina 1: 700 s activity, 34 waves, α = 1.78, p_{α} = 0.23, β = 2.43, p_{β} = 0.11; retina 2: 750 s activity, 25 waves, α = 1.73, p_{α} = 0.0,β = 2.25, p_{β} = 0.08; here, the small sample size may, however, introduced a bias). Together, these results are consistent with the model prediction that retinal waves exhibit statistical properties best described by powerlaw functions.
One retina from a P5 mouse did show bimodal size and lifetime distributions (Fig. 4E), and similar results were obtained in older animals (data not shown). Experimental observations suggest that around day P5 in mouse, GABAergic transmission begins to influence RGC spiking, increases spatial correlations and has subtle effects on wave spatial extent and frequency (Wang et al., 2007; E. Sernagor, unpublished data).GABA has a depolarizing effect at this developmental stage due to elevated [Cl^{−}]_{i}, so it is possible that maturation of the GABAergic drive causes an effective synaptic strengthening. In this case, the models predicts, as observed, bimodal distributions and increased spatial correlations (compare Fig. 3A).
Percolation in the retinal network
Powerlaw event size distributions are also called scalefree distributions, because their shape is unaffected by spatial scaling, such as a multiplication by a constant factor. Then, waves of all possible sizes are produced by the retina, ranging from very frequent waves involving just a few SACs to the very rare events of waves covering the whole retina. In other words, at this point, the properties of the desynchronized and the synchronized network states, characterized by small and large waves respectively, coexist. In a more abstract sense, this constitutes the typical signature of a continuous phase transition between an ordered and disordered network state. Similar phase transitions have been previously studied in related lattice models (Bak et al., 1987; Hopfield and Herz, 1995; Vespignani and Zapperi, 1998), and were suggested as an explanation for powerlaw event size distributions measured in cortical slices and cultures (Beggs and Plenz, 2003).
As outlined above, the sAHPmediated desynchronization is a key mechanism to control the network state. It is activated by neural activity and temporally reduces the neural excitability. The size of a wave is, therefore, always determined by the number of excitable SACs available in the network, which in turn depends on the rate of waves and intrinsic bursting in SACs, and the size of previous waves. Therefore, on average, the central parameter that determines the wave properties is the density of excitable neurons. A system with these properties can be analyzed with methods from percolation theory, which was developed to describe the spreading of activity through a medium, where motion is solely determined by the stochastic properties of the medium (Essam, 1980). The corresponding mechanisms in our model are the propagation of neural activity through synapses, where the state of the neurons, excitable or refractory, dictates where activity can spread and where it cannot.
Importantly, percolation models display a continuous phase transition. When the density of excitable neurons is low, activity cannot spread far and waves are small. Increasing the density then also increases the mean wave size. At a particular density (the percolation threshold), the system undergoes a transition from a state of purely local connectedness to a state, where connections between neurons across the whole network can exist. This point is characterized by a diverging characteristic length scale, and as a result, the mean size of the clusters of connected neurons diverges. Then, propagating activity will assume powerlaw distributions. In particular, as observed here, the expected mean field exponents are α = 1.5 for wave sizes and β = 2 for their lifetimes (Muñoz et al., 1999).
To describe our model in terms of a percolation model, we define three activity states of SACs: stable (S) for a refractory neurons, critical (C) for a nonrefractory neuron, and active (A) for a bursting neuron, with associated densities ρ_{s} and ρ_{c}. The following state transitions are possible (Fig. 5A): Stable to critical with a probability p (S: p → C), and critical to active either with probability fρ_{c}, if the neuron participates in a wave, or with probability r in isolation, if it produces an intrinsic burst (C: fρ_{c} → A and C: r → A). Each firing event is always followed by a transition into the stable state, and firing events are brief and rare (p ≪ 1); hence, we can write C: fρ_{c} → S and C: r → S. Assuming local homogeneity and ignoring higher orders in ρ_{a}, the temporal evolution of ρ_{c} is then given by the following:
with ρ_{s} + ρ_{c} = 1. Here, p is not necessarily a fixed quantity and may depend on the previous activity. This allows for an estimation of the density of excitable neurons (the critical state density) through quantities that can be measured in simulations. To this end, we consider what happens after a burstor wave has occurred. Then, we obtain ρ_{c}^{+w} = ρ_{c} − fρ _{c}^{2} after a wave, and ρ_{c}^{+b} = ρ_{c} − rρ_{c} after an intrinsic burst. Assuming stationarity, i.e., ρ_{c}^{+w} ≈ ρ_{c}^{+b}, yields the estimate ρ_{c}^{+} ≈
The estimated density of excitable neurons ρ_{c}^{+} as a function of the noise intensity is shown in Figure 5B. For weak noise, ρ_{c}^{+} is small, indicating that most cells simultaneously become refractory due to the predominantly very large waves. ρ_{c}^{+} increases with increasing noise strength, consistent with a higher rate of intrinsic bursting and a more variable wave size. At the points where the wavesize distribution best approximates a power law with an exponent of α = 1.5 (Fig. 5C), the density of critical states is ρ_{c}^{+} ≈ 0.7. This value is very close to the threshold for site percolation on a hexagonal lattice (p_{c} ≈ 0.6970) (Suding and Ziff, 1999). Hence, wavesize powerlaw distributions are found at the point where our approximation predicts a percolation phase transition. When the noise intensity is further increased, ρ_{c}^{+} also increased and eventually exceeds unity. Then, the approximation is invalid because the strong noise causes bursting in refractory neurons that cannot participate in waves.
In conclusion, these results show that the predicted and experimentally observed powerlaw statistics are consistent with a network state close to a percolation phase transition. Here, the dynamic regulation of the network excitability maintains a state at the boundary between a purely local and global connectedness of the SAC network.
Discussion
Returning to the issues raised at the beginning of this study, our results show that earlystage retinal waves indeed reflect a very specific network state characterized by highly variable spatiotemporal patterns, the absence of cyclic behavior, and powerlaw distributed wave sizes and lifetimes. Theoretical arguments show that these properties result from a dynamic regulation of the neural excitability, which maintain the SAC network close of a percolation phase transition. At this point, the network is at the transition point between a state of purely local and global functional connectedness. Local functional connectedness yields strongly desynchronized activity and small waves, and in the globally connected state, activity is highly synchronous with proportionally more large waves events. At the transition point, the neural activity simultaneously reflects both states, and therefore exhibits maximal variability with respect to wave sizes and durations. This is a defining feature of retinal waves (Feller et al., 1997), and we confirm the predicted powerlaw behavior with MEA recordings of retinal waves. Moreover, we demonstrate with the model how different pharmacological interventions can specifically shift the network away from this state.
Developmental significance
There are at least two reasons why such activity patterns may be relevance to neural development. The first is based on the consideration of what may constitute a stable state in the highly recurrent network in the immature retina and how it can be maintained and controlled by activitydependent mechanisms. In particular, there may be two competing demands. On the one hand, the network has to maintain a certain level of activity, to ensure neural survival, to form and stabilize synapses, and to trigger cellintrinsic activitydependent developmental molecular pathways. Then, the synchronized, globally connected phase appears to be a suitable regimen. Yet, high synchrony also renders the network very unstable to small perturbations and leads to epileptiform activity. Therefore, it is at the same time necessary to limit the amount of synchrony. These two requirements are exactly balanced at the point of the phase transition. A similar effect has also been demonstrated recently in a randomly driven network model, where the stability and sensitivity to input was maximized at a transition from a quiescent regimen to a state with propagating, selfsustained activity (Ribeiro and Copelli, 2008). Since such behavior is also seen in cultured networks (see below), this may be a state which neural networks without an external drive generally adopt during development.
Second, it is now well established that spontaneous neural activity during development plays an important role in sculpting out the mature neural circuitry in various neural systems (Moody and Bosma, 2005). Here, we show that earlystage retinal waves exhibit maximal variability and containevents on all length scales due to their scalefree character. This randomness ensures that activity patterns are not biased toward a typical scale or sequence of events. This may be a general advantage while the retinal circuitry develops and retinal projections are established in visual centers of the brain. Indeed, the different correlation structure in earlystage waves in mice lacking the β2AChR subunit (Sun et al., 2008), for example, may cause a sustained lack of refinement of retinocollicular projections (McLaughlin et al., 2003) and prevent the development of eyespecific thalamic layers (MuirRobinson et al., 2002). Since our results suggest that earlystage retinal waves reflect a rather precisely regulated network state, it would be interesting to assess how sensitive higher visual system development is to modest changes of the spatiotemporal wave properties. In this context, it is interesting to note that the spatial frequency content in natural images also decays with a power law (Field, 1987), so retinal waves might present the visual system with an early opportunity to adapt to such patterns.
Clearly, this network state is not maintained during development. Changes toward a stronger emphasis on feedforward processing, as in the mature retina, might begin as soon as GABAergic synapses first appear, where we observed first deviations from the powerlaw behavior. At later stages, when glutamatergic transmission drives retinal waves, waves are small and patchy (Zhou and Zhao, 2000; Sernagor et al., 2003), which in our model would correspond to a state well below the percolation threshold for lateral activity propagation.
Relations to other systems and previous models
Powerlaw distributed event sizes and lifetimes with the same exponents (α = 1.5 and β = 2) were also found in cortical slices and cultures (Beggs and Plenz, 2003, 2004), and recently in the developing cortex in vivo (Gireesh and Plenz, 2008). These studies have proposed that this phenomenon could be explained in terms of a branching process, where a critical state exists when each spike in one neuron causes on average exactly one spike in another neuron. A critical branching process predicts powerlaw event size and lifetime distributions with the observed exponents, and is directly related to a large class of models that display the same type of phase transition (Zapperi et al., 1995). In models of neural networks, the same behavior has been observed in slowly driven recurrent network with nonleaky or nonrefractory neurons (Hopfield and Herz, 1995; Eurich et al., 2002; Buice and Cowan, 2007), but theoretical work also shows that it may be difficult to reconcile this behavior with the biophysical properties of neurons (Dickman et al., 2000) (but see Levina et al., 2007). Here, we have identified a biologically plausible model that reproduces these properties. It is worth noting that the models mentioned above are generally different to the one presented here. Mathematical analysis has shown that the lattice all belong to the class of directed percolation models, where activity propagates forward in time (Vespignani and Zapperi, 1998; Buice and Cowan, 2007). In the present model, in contrast, the refractory mechanism prevents spreading of activity into previously active regions, and it is therefore related to dynamic isotropic percolation.
A refractory mechanism as the basis for the regulation of the spatiotemporal properties of retinal waves was first suggested in a model by Feller et al. (1997) and recently confirmed experimentally (Zheng et al., 2006). As in the present study, Feller et al. (1997) concluded that wave trajectories depend on the density of nonrefractory amacrine cells, which is determined by their activation history due to previous waves. This earlier model, however, required an additional filtering of amacrine cell activity through spatiotemporal integration and thresholding in RGCs to obtain compact waves. Here, we show that, in line with recent experimental and theoretical results (Zheng et al., 2006; Godfrey and Swindale, 2007), the properties of retinal waves are already present at the level of the SAC network. In addition, the model by Feller et al. (1997) required a certain degree of variability in the decay times of the refractory mechanism to prevent synchronization (Butts et al., 1999). A similar solution was recently suggested in a modeling study by Godfrey and Swindale (2007), where a highly nonlinear mechanism was used to amplify difference in refractory times between cells in the center and near the borders of waves. In this lattice model, the resulting network exhibits chaotic behavior, which may be a consequence of these severe nonlinearities.
Here, we reconcile and extend these previous approaches by suggesting a physiologically plausible mechanism for network desynchronization, which is based on differences between intrinsic bursting and depolarizations in SACs (Zheng et al., 2006), and allows us to provide a general theoretical framework to explain the origins of retinal waves. Given its generality, this principle might also apply to other development neural systems.
Footnotes

This work was supported by the Edinburgh Compute and Data Facility [partially supported by the eDIKT (eScience Data, Information and Knowledge Transformation) initiative]. We thank members of the CARMEN (Code Analysis, Repository & Modelling for ENeuroscience) consortium for providing infrastructure for data sharing. M.H.H. was supported by a United Kingdom Medical Research Council Fellowship (G0501327), and E.S. was supported by grants from the Engineering and Physical Sciences Research Council (EP/E002331/1) and the Centre of Excellence for Life Sciences Limited (OneNorthEast). We thank Drs J. Cortes, S. Eglen, M. Herrmann, M. Muñoz, and M. van Rossum for discussions.
 Correspondence should be addressed to Matthias H. Hennig, Institute for Adaptive and Neural Computation, School of Informatics, University of Edinburgh, 10 Crichton Street, Edinburgh EH8 9AB, UK. m.hennig{at}ed.ac.uk