Propagating Waves Can Explain Irregular Neural Dynamics

Dynamic wave patterns emerging from balanced, spatially extended networks

The balance of a neuron often refers to the approximately equivalent amounts of inhibitory and excitatory synaptic inputs that it receives (van Vreeswijk and Sompolinsky, 1996, 1998; Salinas and Sejnowski, 2000; Xue et al., 2014). This can be measured by the ratio β = <c̅i̅j̅I̅> / <c̅i̅j̅E̅> where excitatory current c_ij^E(t) = g_ij^E(t)|V_ij(t) − V_E|, inhibitory current c_ij^I(t) = g_ij^I(t)|V_ij(t) − V_I|, <> indicates averaging over time excluding refractory periods, and the overbar represents averaging across neurons. When β ≈ 1 (i.e., excitatory and inhibitory inputs are approximately equilibrated), we refer to this as the balanced condition. When β ≈ 1, the randomly coupled network (see Materials and Methods) generates asynchronous dynamics without any structured patterns, as found in previous balanced networks (van Vreeswijk and Sompolinsky, 1996; Renart et al., 2010). In this balanced state, individual neurons fire irregularly (Fig. 1a), and the distribution of membrane potential is Gaussian.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

The formation of propagating wave patterns with complex dynamics depends upon a balance of excitation and inhibition, and upon spatially extended coupling. Colors in the images represent membrane potential between inhibitory reversal V_I = −80 mV (blue) and threshold V_T = −55 mV (red), as indicated by the color bar. Each image is of a network of 300 × 300 neurons, but only a 100 × 100 region is shown for clarity. a, Snapshot of spatiotemporal activity in the randomly connected, balanced network. For this case, there are no structured patterns. b–d, Snapshots of spatiotemporal patterns in the balanced, spatially extended network. Patterns take the form of multiple, localized, crescent-shaped propagating waves and wandering patches, which are outlined by the green boxes. The patterns may interact in a variety of ways, such as mutual or partial annihilations, in which one or both patterns disappear; deflection, in which patterns move in opposite directions; and fission, in which a wave splits into two. e, A spiral wave formed in the unbalanced network in which excitation dominates; the center of the spiral wave drifts across the network. f–h, Snapshot of spatiotemporal activity when the balanced condition in a network is broken so that excitation dominates. The separation between the waves decreases, and then nearby waves merge; this effect persists until plane waves form. i, Average width of the waves as a function of time. The balanced condition is broken at 3 s, and for this case the excitation is dominant. j, Spatiotemporal activity when balance in a network is broken so that inhibition dominates; patterns tend to dissipate quickly such that they cannot propagate over long distances. k, Propagation distances of the wave patterns as a function of time. The balanced condition is broken at 3 s, and inhibition dominates in this case. l, Snapshot of the spatiotemporal activity of the spatially extended network when noisy inputs are included. Despite the presence of random inputs, propagating wave patterns are still able to form from the network.

In contrast, the balanced, spatially extended network generates spatiotemporal patterns with complex dynamics. As shown in Figure 1b, c, d, these patterns include several localized, crescent-shaped waves that move rapidly, and localized, patchy patterns that slowly wander around (outlined by the green boxes). The initiation sites and subsequent trajectories of both types of pattern are seemingly random, and the patterns interact when they approach each other; these interactions include repulsive collisions, and partial and full annihilations. Because both types of pattern can propagate over the network, we generally refer to them as propagating wave patterns. Notably, propagating wave patterns with random initiation sites and propagation directions have been observed in the spontaneous activity of visual cortex of rats (Han et al., 2008, their Fig. 2b) and in the sensorimotor cortex of behaving mice (Ferezou et al., 2007). Multiple wave patterns that sweep across the cortex have been found in the visual cortex of cats (Arieli et al., 1995, 1996), and the coexistence of waves and localized patchy patterns has been observed in the barrel cortex of rodents (Petersen et al., 2003).

When the network input is unbalanced (β ≤ 0.95) and disinhibited, wave patterns with more regular, global properties emerge. In particular, global ring or target waves occur periodically, with frequency proportional to the external input. Alternatively, depending on the initial conditions, the network sometimes develops plane waves that propagate in a preferred direction, or spiral waves, the centers of which can spread across the entire lattice (Fig. 1e). This is largely consistent with recordings in cortical slices and intact neocortex of rats: when neurons were disinhibited using bicuculline to block GABA receptors, alternative patterns of ring, plane, and spiral waves were observed (Huang et al., 2004, 2010).

To further illustrate that balanced excitation and inhibition are needed for the formation of propagating wave patterns with complex dynamics, we break the balanced condition in the network to study the changes in these patterns. For this purpose, the model is run for the usual transient time and then after 3 s, the inhibitory connection strength W_I is altered to obtain either an excitation or an inhibition dominated case. For the excitation dominated case, the typical changes of these patterns after breaking the balance are shown in Figure 1f, g, h: the separation between these localized wave patterns decreases, and then nearby waves merge with each other; this process continues until eventually global plane waves are formed (Fig. 1h). To quantify such changes, we calculate the average width of the wave patterns identified by our algorithm (see Materials and Methods) as a function of time. Figure 1i shows that the size of the waves increases rapidly after breaking the balanced condition and continues to increase until all the patterns are merged into a plane wave. On the other hand, after breaking the balanced condition to obtain an inhibition dominated case, the wave patterns change from ones with relatively long range propagation to ones tending to disappear shortly after they are formed (i.e., they are short-lived due to the strong inhibition) (Fig. 1j). Based on this observation, we calculate the average propagation length (500 ms time windows) of each pattern as a function of time. Figure 1k shows that, as expected, the propagation length of the waves rapidly decreases after the balance is broken until it reaches a certain value that is quite small relative to the size of the network.

In our network model, the balanced state and the resultant propagating wave patterns emerge without fine tuning of parameters. For instance, changing the values of F^E and F^I still leads to the balanced state as long as F^E is high enough relative to F^I such that individual neurons can fire even without recurrent synaptic inputs. The population-averaged firing rates increase nearly linearly as a function of external input strength (F^E), as has been found for randomly connected, balanced networks (van Vreeswijk and Sompolinsky, 1996). Additionally, when the neurons in our network receive stochastic inputs (see Materials and Methods), they still settle to a balanced state and the network is able to form propagating wave patterns (Fig. 1l), as found in the network with a constant external input. Moreover, the strengths of recurrent excitation (W_E) and inhibition (W_I) have a relatively broad range over which the balanced condition is preserved; for instance, given W_E = 0.23, the network is balanced for 0.23 ≤ W_I ≤ 0.35. In a recent study of neural field models (Rosenbaum and Doiron, 2014), it was found that, in the continuum limit, balanced firing rate solutions require that the spatial spread of external inputs be broader than that of recurrent excitation, which in turn should be greater than or equal to that of recurrent inhibition; however, this condition is not satisfied for finite-size networks, and crucially, propagating wave activity cannot be captured by the mean field models, as pointed out by Rosenbaum and Doiron (2014). In our study of the finite-size spiking neural network with emergent propagating wave patterns, the balanced condition with the ratio of average excitatory current and inhibitory current close to 1.0 is quite robust, and these wave patterns underlie variable firing activity of neurons as shown in the following section.

Spikes of neurons exhibit a doubly stochastic process

We next characterize the dynamics of the propagating wave patterns and show that these dynamics cause the doubly stochastic process of neural firing activity, with variability in spike timing and slow fluctuations in firing rates (Churchland et al., 2010, 2011). First, the center of mass of each pattern is calculated and tracked over time (see Materials and Methods). Figure 2a, b shows the typical trajectories of crescent-shaped waves and patchy patterns, respectively; both patterns appear to have random initiation sites and trajectories. They have distinct dynamics, with crescent-shaped waves propagating over long distances and patchy patterns wandering around a particular area. By calculating the displacement of centers of mass for each pattern, we find that the crescent-shaped waves have a mean speed of 2.0 ± 0.6 gridpoints/ms (Fig. 2c). We also change the connectivity density (i.e., the number of connections each neuron receives over a given distance) and find that the average speed of the waves increases monotonically as the density increases.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Dynamic properties of propagating wave patterns. a, A sample of the trajectories of crescent-shaped waves. Both the initiation sites (black circles) and propagation paths (gray lines) of these waves appear to be random. The trajectories are terminated as a result of interactions with other waves or as a result of the wave fissioning. b, Trajectory of a typical patchy pattern with initiation site depicted by the circle. The pattern tends to “wander” around a particular area and covers a smaller region than the crescent-shaped waves (note the smaller 20 × 20 scale of this figure). c, Distribution of the speeds of the crescent-shaped waves, which has a mean of 2.0 gridpoints/ms. d, Mean-squared displacement MSD(τ) as a function of time increment τ. Red line indicates the MSD for crescent-shaped waves; a fit to the straight-line part in the log-log plot (indicated by the black line above) gives MSD(τ) ∝ τ^α with α = 1.9 (nonlinear least-squares fitting), indicating the superdiffusive dynamics of crescent-shaped waves. To confirm that the crescent waves propagate in a superdiffusive way rather than a ballistic way (α = 2), we compare our fit against a quadratic fit by calculating the Bayes factor of B = P(D | M)/P(D | M₀) where P(D | M) is the likelihood of the power law fit M with α = 1.9 given the data D, and P(D | M₀) is the likelihood of the quadratic fit M₀; we obtain B > 1000, indicating that the superdiffusive case is much more likely than the ballistic case. Blue line indicates the MSD for patchy patterns: for t > 10 ms, the slope of the line of best fit (the black line above) is α = 1.0, which indicates Brownian dynamics over long time scales.

To further quantify the dynamics of these patterns, we calculate their MSD (see Materials and Methods). For the crescent-shaped waves, MSD appears to be linear on a log-log scale (Fig. 2d, red line), which implies that it is a power function of time increment τ such that MSD(τ) ∝ τ^α. By using nonlinear least-squares fitting via the Levenberg-Marquardt algorithm (Press et al., 2007), we obtain an exponent of α = 1.9. For Brownian motion, the MSD is a linear function of time (α = 1); when α ≠ 1, a diffusion process has a nonlinear relationship to time and is called anomalous diffusion. In particular, if 1 < α < 2, it is a superdiffusive process; and if 0 < α < 1, it is a subdiffusive process (Metzler and Klafter, 2000). This result therefore indicates that the crescent-shaped waves move in a superdiffusive manner. Anomalous superdiffusion indicates that the motion of these waves is not fully random as is the case for Brownian motion and that, instead, there are long-range correlations in their propagations (Metzler and Klafter, 2000). For the patchy patterns that wander around, MSD(τ) ∝ τ^α and α = 1.0 (using nonlinear least-squares fitting) for τ > 10 ms, indicating that these patterns move in a normally diffusive manner over long time scales and that they move much slower than crescent-shaped waves (Fig. 2d, blue line).

Figure 3a shows that neurons within the network exhibit dynamic transitions between low-activity and high-activity states. These transitions appear to coexist with variable spike timing of individual neurons, which is a characteristic feature of a doubly stochastic process (Churchland et al., 2011; Churchland and Abbott, 2012; Litwin-Kumar and Doiron, 2012). To quantify the variability of spike timing, we first calculate the CV of ISIs as outlined in Materials and Methods. The average value of CVs across all trials is 1.1 ± 0.1, which is in agreement with experimental studies (Softky and Koch, 1993; Holt et al., 1996). However, a high CV may reflect fluctuations in firing rates as well as variability in the timing of individual spikes. To quantify such fluctuations in firing rates, we calculate the spike count Fano factor (see Materials and Methods) for each neuron by using a fixed time window of Δt = 100 ms, which yields a mean value of 1.4 ± 0.5. Fano factors with values >1.0 have been observed in many cortical systems (Tolhurst et al., 1983; Britten et al., 1993; Churchland et al., 2011) and are evidence for the fluctuations of neural firing rates. The fluctuations of firing rates over long time scales are further demonstrated by calculating Fano factors as a function of time window. As shown in Figure 3b, Fano factor increases linearly as Δt increases; this is because of the transitions between high and low-activity states that happen over long time scales (Churchland et al., 2011). A recent modeling study of balanced networks with clustered connections has found similar doubly stochastic firing activity and transitions between high and low firing activity states (Litwin-Kumar and Doiron, 2012), in which such transitions are caused by spontaneous switching dynamics between different attractors. However, for our case, the dynamics of the propagating wave patterns cause these transitions; this mechanism is illustrated below.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Properties of the spiking activity of network neurons. a, Raster plot showing the spike times of a subpopulation of 50 excitatory neurons. Firing activity with variable spike timing is punctuated by dynamic transitions between high and low firing activity states. b, For the balanced, spatially extended network, Fano factors increase as a function of length of time window (blue line). The neurons in the balanced, randomly connected network (red line) have a high Fano factor, but it is independent of time window. For the disinhibited network (green line), the Fano factor is independent of time window and has a value of almost zero. c, Histogram of the number of neurons in a 5 × 5 subgrid that would fire when there is a crescent-shaped wave passing by the subgrid. d, Percentage of patchy patterns as a function of random rewiring probability P. In the balanced, spatially extended network, the percentage of patchy patterns increases with P. e, Average Fano factor as a function of random rewiring probability P. In the balanced, spatially extended network, the Fano factor increases with P. f, Raster plot showing the spike times of a subpopulation of 50 excitatory neurons in the spatially extended network with 10% of the connections randomly rewired.

Propagating wave patterns generate irregular spiking activity

We next consider the relationship between the spiking activity of individual neurons and the collective dynamics of propagating wave patterns. By tracking the motion of individual patterns, we observe that, when a crescent-shaped wave front moves toward a certain area, neurons in the area would receive synaptic inputs from neurons that are firing in the front, and their membrane potentials are depolarized such that they may subsequently fire. Because the wave moves in a superdiffusive way, it leaves that location rapidly, which means that neurons around that location fire no more than once. To verify this, we consider a 5 × 5 subgrid with 25 neurons, which is smaller than the area of a typical wave pattern and thus cannot contain more than one wave trajectory at a time. We then calculate over all trials the number of neurons within this grid that would fire when a wave passes by. The resultant distribution, shown in Figure 3c, has a mean value of 11.2, meaning that, when a wave passes by, only a fraction of neurons in the grid fire. The complex dynamics of the waves, namely, their random origins and trajectories (Fig. 2a), results in the irregular passage of these wave fronts around that location; this causes the high variability of the spike times indicated by the CV. Because the waves tend to leave each location rapidly and cause only a fraction of neurons in a grid to fire, their passage corresponds to the relatively low-activity state. This can be illustrated by choosing a neuron that has transitions between the low and high-activity states, and we find that, during the time periods when the neuron is not within a patchy pattern and crescent waves occasionally pass the location of the neuron, the mean firing rate of the neuron is 16 Hz. On the other hand, when a patchy pattern moves to this location, the pattern wanders around the location due to its normally diffusive dynamics, and drifts away slowly; this enables the neuron within that location to fire persistently, corresponding to the high-activity state as shown in Figure 3a. For the time periods when the example neuron is within a patchy pattern, we find that the mean firing rate is 137 Hz, and similar values can be found for other neurons, which are passed by patchy patterns. The low-activity state is therefore associated with crescent-shaped waves, whereas the transition to the high-activity state is caused by the approach of a patchy pattern.

Based on the above result, one would expect that the time scale of these transitions is closely related to the density of the patchy patterns. To test this prediction, we randomly rewire connections in the network with a certain probability P (see Materials and Methods). Despite the random connections, the localized patterns shown in Figure 1b–d can still emerge for P ≤ 10%. Figure 3d shows that, for the network with random connections, the proportion of patchy patterns relative to the number of crescent-shaped waves is greater than that in the original model; indeed, the proportion increases monotonically as the rewiring probability P increases. This increase in the proportion of patchy patterns results in an increase in Fano factors (Fig. 3e). This is because the greater the proportion of the patchy patterns, the greater the frequency of transitions between high and low-activity states, as shown in Figure 3f. The Fano factor asymptotically approaches a particular value, which corresponds to an absence of the crescent-shaped waves. Because the crescent-shaped waves are mostly responsible for the low-activity state, their absence means that neurons exhibit high firing rates or no firing at all; this results in greater fluctuations of firing rates.

To show that the propagating wave patterns emerging from the balanced network are essential for doubly stochastic firing activity, we calculate the Fano factors of the unbalanced network and the random network. The global ring, plane, and spiral waves (Fig. 1e,h), which occur in the unbalanced network, cannot cause variable spike timing. In this case, the Fano factors of firing activity are not only independent of time windows but also have very low values (Fig. 3b, green line). The randomly coupled network, which cannot generate any structured patterns (Fig. 1a), can produce spike count variability, but it does not exhibit slow fluctuations in firing rates; as shown in Figure 3b, whereas its Fano factor is close to 1, it is independent of time window.

Previous models of balanced networks have shown that averaged spike count correlations of pairs of cortical neurons form a wide distribution (i.e., with a high SD) with a mean value close to zero (Renart et al., 2010; Litwin-Kumar and Doiron, 2012). We calculate the correlation coefficients (see Materials and Methods) of spike counts of pairs of randomly selected neurons in the balanced, spatially extended network. Figure 4a shows a wide distribution of correlations with a mean of 0.0008 ± 0.0002, similar to that found in Litwin-Kumar and Doiron (2012). We also calculate the correlation coefficients as a function of distance (Fig. 4b): at relatively short distance, there are positive correlations that decrease as the distance increases, which is largely consistent with the distance-dependent neuronal correlations as observed in Smith and Kohn (2008) and Ecker et al. (2014). The correlations are negative at intermediate distances and approach zero when distance is large enough; this is consistent with the negative correlations found for the randomly selected neuron pairs (Fig. 4a) (Renart et al., 2010; Ecker et al., 2014). Negative correlations have also been observed by Vaadia et al. (1995).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Spike count correlations between pairs of neurons in the balanced, spatially extended network. a, Histogram of the correlation coefficients of random pairs of neurons; it is wide and centered around zero. b, Average correlation coefficients, computed over all neuron pairs, as a function of their distance. c, Histogram of spike count correlations for nearby neurons (separated by 5 gridpoints), which has a mean value of 0.66, as indicated by the arrow. d, Histogram of spike count correlations for nearby neurons (separated by 5 gridpoints) when they are in the high-activity state; the mean value is 0.22, as indicated by the arrow. e, Histogram of spike count correlations for nearby neurons (separated by 5 gridpoints) when they are in the low-activity state; the mean value is 0.15, as indicated by the arrow.

As illustrated above, there are transitions between low- and high-activity states, with the latter being produced by the slowly propagating patchy patterns. That is, when a patchy pattern moves to a location in the network, it drifts away from the location slowly. This makes the neurons within the location maintain sustained firing states for a relatively long period of time, resembling the “up state” found in the cortex. To show whether the correlations depend on the states of neuron firing activity for nearby neurons, as observed in Ecker et al. (2014), we also calculate the distribution of spike count correlations when patchy patterns are over a neuron pair (high-activity state), and when they are not (low-activity state); these are then compared with those obtained over all trials. To do this, we classify the high firing state as those Δt = 50 ms intervals with spike counts corresponding to a firing rate of >50 Hz. The mean values of the spike count correlations are 0.66, 0.22, and 0.15 for all trials (Fig. 4c), for the case when both neurons are in a high firing state (Fig. 4d), and for the case when they are not (Fig. 4e), respectively. The mean value of the correlation coefficients of the high-activity state is significantly larger than that of the low-activity state (p < 0.001), and both are significantly smaller than that for all trials (p < 0.001). These results therefore indicate that, in our model, the spike count correlations depend on the states of neural firing activity.

In total, the results we have obtained show that the balanced, spatially extended network with emergent propagating wave patterns is able to capture salient characteristics of variable firing activity of cortical neurons, including their variable spike timing, slow fluctuations of firing rates, small spike count correlations of randomly selected neuron pairs, and state-dependent spike count correlations.

Synchronized synaptic inputs to individual neurons

In the input synchrony model proposed for explaining variable neural dynamics, synaptic inputs to individual neurons consist of large, brief, and occasional inputs or “bumpy” inputs, indicating the synchronous arrival of many spikes of presynaptic neurons (Stevens and Zador, 1998; Okun and Lampl, 2008). These kinds of synchronized bumpy synaptic inputs have been found in the cortex (Hasenstaub et al., 2005; DeWeese and Zador, 2006; Okun and Lampl, 2008; Poulet and Petersen, 2008; Tan et al., 2014), but the mechanisms that give rise to such inputs have not been studied. We now show that, in our balanced, spatially extended network, the emergent propagating wave patterns provide a network mechanism to generate such inputs.

Figure 5a shows the dynamics of the excitatory synaptic conductance received by individual neurons in our model, which are characterized by quiescent periods punctuated by large, brief, and occasional excursions in amplitude, with variable magnitudes and shapes. Figure 5b shows the distribution of the excitatory synaptic conductance g_ij^E(t) received by each neuron over the whole time course of one trial; the distribution has a heavy tail. We use a linear-log plot because we can use it to directly compare with empirical observations in DeWeese and Zador (2006) and Hromádka et al. (2013), but it is evident that the distribution is heavy-tailed if the scale of both axes is linear. It is interesting to note that heterogeneous synaptic inputs with heavy-tail distributions have been found in the somatosensory cortex (Lefort et al., 2009). To quantify these non-Gaussian dynamics, we calculate the kurtosis for the distribution of the excitatory synaptic conductance of each neuron within the network (see Materials and Methods). We find that these values display a wide range of activity with a maximum kurtosis of 17.6 and a minimum kurtosis of −1.1. However, negative values of kurtosis are rare; most values are positive, with a population mean of 5.9 (SD 1.6). These kurtosis values are comparable with those implied by whole-cell recordings of auditory cortex in awake rats (Hromádka et al., 2013). Because the kurtosis is a measure of the proportion of samples within the peak and tail of a distribution compared with a normal distribution, the kurtosis is positive for heavy-tail distributions (see Materials and Methods). The negative values of kurtosis indicate that distributions of the inputs may occasionally be bimodal; such distributions may result from periods of quiescence punctuated by periods of sustained input.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Properties of synaptic inputs to neurons in the balanced, spatially extended network. a, Time series of excitatory synaptic conductance g_ij^E(t) received by a neuron over 1 s; it consists of quiescent periods punctuated by brief excursions in conductance, indicating the arrival of highly synchronized inputs of many neurons. The time courses of g_ij^E(t) have variable shape and duration; for instance, at ∼200 ms, the bump in synaptic conductance has a relatively longer duration. b, Histogram of excitatory synaptic conductance g_ij^E(t) received by each neuron over the whole time course of one trial (i.e., the trace of the conductance as shown in a). The bumpy synaptic input results in a heavy-tail distribution of conductances. Note the logarithmic scale of the coordinate; this is used for direct comparison with relevant empirical data, but the heavy tail is still evident for a plot with both axes linearly scaled. c, Excitatory (solid line) and inhibitory (dashed line) inputs received by a test neuron in the model of a wave front of width N_w ≈ 19 approaching a test neuron head-on. d, The conductances received by two nearby neurons in the network, shown as red and blue lines, as a wave passes by. Solid lines indicate excitatory conductance. Dashed lines indicate inhibitory conductance. e, Average cross-correlation of synaptic inputs received by neuron pairs separated by an interneuron distance of 5. There are high cross-correlations of 0.69 (SD 0.04) for excitatory-excitatory (E-E), 0.56 (SD 0.05) for inhibitory-excitatory (I-E), and 0.91 (SD 0.02) for inhibitory-inhibitory (I-I) conductances, respectively. Red bars represent maximum cross-correlation. Blue bars represent cross-correlation at zero time lag. f, The conductances received by two nearby neurons in the network as a wave approaches them from the opposite direction to that for d. g, Cross-correlation between the inhibitory input and the excitatory input found for nearby neurons in the network (black line) and the cross-correlation calculated for the model of a wave approaching a test neuron (red line).

We next illustrate that the dynamics of propagating wave patterns in our network provide a mechanistic explanation for the bumpy, synchronized synaptic inputs, as shown in Figure 5a. It is apparent that a wave front, which includes multiple spiking neurons, provides a source of synchronized input to any neurons that it approaches. However, the superdiffusive dynamics of the crescent-shaped waves means that such synchronized input is transient, as the wave will pass by and move out of range of the neuron rapidly. We can quantify the number of spikes participating in a bumpy synaptic input by considering a planar wave front approaching a “test” neuron. The maximum distance of excitatory connections is D^E = 10; above this distance, the test neuron will not receive any recurrent input (see Materials and Methods). Therefore, we assume the initial displacement between the wave and the test neuron to be d = −D^E. If the wave approaches the test neuron head-on, the neuron will continue to receive input until the wave passes out of range at d = D^E. Because d is in units of neuron separation (i.e., gridpoints), the test neuron receives excitatory inputs from the wave front (i.e., at 2D^E + 1 locations); this is the range over which the wave provides input to our test neuron as it approaches and then moves away; in other words, it is the length of the wave trajectory that affects the test neuron. The width of a typical wave front consists of N_w ≈ 19 (Fig. 1i) depolarized neurons. The number of neurons that could participate in a bump of synaptic input is the width multiplied by the length of the trajectory (i.e., N_w(2D^E + 1) − 1), where the minus one excludes the test neuron from receiving input from itself. Excluding neurons, which are separated from the test neuron by a distance greater than D^E, we find that 315 neurons can participate in the bump. This is similar to the number of synthetic excitatory postsynaptic potentials used to generate synchronized inputs to cortical neurons in brain slices (Stevens and Zador, 1998), but it is approximately an order of magnitude less than that found in vivo (DeWeese and Zador, 2006). This discrepancy is likely due to the scaling of our network (see Materials and Methods).

The model of a typical wave front approaching a test neuron can also be used to explain the shape of the bumps of synaptic conductance. Consider the wave front moving at the average speed of 2 gridpoints/ms (Fig. 2c) head-on toward a test neuron and then away from it. As excitatory coupling strength has an exponential dependence on distance (Eq. 4, see Materials and Methods), the amplitude of the excitatory inputs first increases and then decreases. The resultant time course of the input produced in this model is similar to the time course (full width at half-maximum is 4.4 ms) and amplitude (∼2000 μS) of the bumps of largest amplitude in the network model (Fig. 5c, solid line). Bumpy synaptic inputs in our network usually have diverse profiles, but such variations can arise from differences in the sizes of firing wave fronts and the complex trajectories of these waves, which often do not approach a neuron head-on but move within its proximity. In summary, the crescent-shaped waves with superdiffusive dynamics can produce bumpy, synchronized inputs with short durations, to individual neurons.

Bumpy synaptic inputs in the balanced, spatially extended network occasionally have longer durations. For instance, the bumpy synaptic conductance centered around 200 ms in Figure 5a has a relatively longer duration, during which the neuron receives synchronized inputs for a sustained period (i.e., up states), similar to intracellular recordings in the auditory cortex of awake rats (Hromádka et al., 2013). These inputs are caused by patchy patterns because, when these patterns move toward a particular neuron, they wander around the neuron and drift away from it slowly in a normally diffusive manner, providing relatively sustained synchronized inputs to the neuron. To further quantify the bumpy, synchronized synaptic inputs, we study their durations, which we define as the time intervals during which excitatory conductance exceeds some threshold. We choose this threshold to be the leak conductance (i.e., g_ij^E(t) > g_L), but as long as the threshold level is not too small, its precise value is not critical; choices near this value yield similar results. Using this definition, we classify bumpy inputs with long duration as those of duration >500 ms, as was used by Hromádka et al. (2013). We find that bumpy inputs of long duration are rare, comprising only 0.05% of all bumps. Furthermore, most time spent away from rest is spent in short bumps; only 1.2% of total simulation time is spent in bumpy inputs with duration >500 ms. Indeed, bumpy inputs of long duration constitute 3.6% of the duration of all bumpy, synchronized inputs, despite, on average, being approximately an order of magnitude longer than bumpy inputs of short duration (Fig. 5a). This is in agreement with experimental results, which suggest that sustained bumps (up states) in inputs are rare in the awake cortex (Hromádka et al., 2013).

Synaptic inputs are correlated and inhibition lags behind excitation

We now demonstrate that our model is also able to capture the cross-correlations of synaptic inputs to pairs of nearby neurons, as measured by using dual-cell recording in Okun and Lampl (2008). Figure 5d shows the typical excitatory (solid lines) and inhibitory (dashed lines) conductances received by two neurons separated by a small distance of 5 gridpoints, as a wave passes by; conductance is positive (Eq. 2 in our model). We calculate the cross-correlations of these inputs over pairs of neurons and over trials and find that there are large peak values of cross-correlation with 0.69 (SD 0.04) for excitatory-excitatory, 0.56 (SD 0.05) for inhibitory-excitatory, and 0.91 (SD 0.02) for inhibitory-inhibitory conductances, respectively (Fig. 5e). This is comparable with the cross-correlations of excitatory and inhibitory inputs found in nearby cortical neurons (Okun and Lampl, 2008; Gentet et al., 2010).

As shown in Figure 5d, excitatory inputs to the nearby neurons reach their peak values at almost the same time; this is also the case for the inhibitory inputs. In contrast, inhibition lags excitation significantly (p < 0.001). As in Okun and Lampl (2008), such lags between inhibition and excitation are then measured based on the peak values of their cross-correlations averaged over neuron pairs and trials; these cross-correlations attain their peak value at non-zero lag, which is 2.3 ± 0.4 ms. This lag is comparable to the lag of 3.8 ± 4.9 ms as found in Okun and Lampl (2008). We also calculate the cross-correlation of excitatory and inhibitory conductances arriving at the same neuron and obtain a lag of 2.2 ± 0.5 ms, comparable with the lag of 2.4 ± 3.6 ms reported by Wehr and Zador (2003) for such inputs received by the same neuron. Our model is therefore able to reproduce a characteristic feature observed in cortical neurons, namely, that inhibition lags behind excitation by several milliseconds. This supports the hypothesis that inhibition controls the integration time window of excitation, enabling neurons to act as coincidence detectors.

In our model, the bumpy excitatory synaptic inputs are always ahead of inhibitory inputs, as shown in Figure 5d, regardless of which direction the approaching wave moves; Figure 5f shows such inputs when there is a wave moving in the opposite direction to that in Figure 5d. This effect is a result of the distance-dependent coupling and the distinct synaptic time scales of excitatory and inhibitory dynamics (τ_d^E = 2.0 ms and τ_d^I = 7.0 ms in Eq. 3). To illustrate this, we consider the model mentioned earlier, of a propagating wave front consisting of N_w = 19 firing excitatory neurons, and extend it to include inhibitory neurons. Because inhibitory neurons only constitute 20% of the neurons in the network, there are half as many inhibitory neurons as there are excitatory neurons along both dimensions. Therefore, the firing wave front consists of ∼N_w/2 inhibitory neurons, which occur at d = −D^I−D^I −2, …, −1, 1, …, D^I −2, D^I, where D^I = 15 is the spatial range of inhibitory neurons. Using this extended model, we obtain the time course of inhibitory input received by our test neuron, shown by Figure 5c (dashed line); this is similar to what is observed in our balanced network (Fig. 5d,f). When we calculate the cross-correlation between excitatory and inhibitory bumps, we find that it obtains its peak value at 3.9 ms (Fig. 5g, red line), similar to the lag found when a wave passes a neuron in the balanced, spatially extended network (Fig. 5g, black line). Based on the model of a wave front approaching a test neuron, we also find that the time lag remains the same for a wave moving in an opposite direction (i.e., starting at d = D^I and ending at d = −D^I). However, if we run the model with the same synaptic time scales (τ_d^E = 2.0 ms and τ_d^I = 2.0 ms), no lag exists in the cross-correlation function; if we run the model with the same distance-dependent coupling for excitatory and inhibitory neurons (K_ij,i′j′^E = K_ij,i′j′^I in Eq. 4), the lag is only 2.3 ms. These results therefore suggest that the lags between excitation and inhibition are dependent on the interplay of the synaptic time scales and the distance-dependent coupling but are independent of the propagation direction of the wave patterns. We also find that this lag is not evident for the balanced, randomly connected network, in which no wave patterns can form.

Membrane potential dynamics are non-Gaussian

In previous models with balanced excitation and inhibition, membrane potential follows a random walk and its distribution is Gaussian (Shadlen and Newsome, 1998); this is also the case for our randomly connected network, in which the membrane potential has a Gaussian distribution with a kurtosis near zero. However, for the balanced, spatially extended network, the synchronized bumpy inputs cause occasional large excursions of membrane potential from rest; this results in a heavy-tail distribution of membrane potential (Fig. 6a). This non-Gaussian dynamic is confirmed by calculating the kurtosis, which has a mean value over all neurons of 5.9 (SD 1.2), close to the kurtosis found in DeWeese and Zador (2006). The distribution of kurtosis values for the membrane potential of each neuron shows that there is a great diversity in membrane potential dynamics, with a maximum kurtosis of 11.8 and a minimum kurtosis of −0.3 (Fig. 6b). We find similar distributions for the synaptic inputs (Fig. 5b), which indicates that there are bumpy dynamics for both synaptic input and membrane potential. This coincides with the passing of waves, which causes highly synchronized synaptic inputs and a resultant sudden increase in membrane potential.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

The dynamics of membrane potential in the balanced, spatially extended network. a, Distribution of the membrane potential of all neurons over the time course of one trial is heavily tailed. Note the logarithmic scale of the coordinate; this is used for direct comparison with empirical data, but the heavy tail is still evident if the scale of both axes is linear. b, Distribution of kurtosis values of membrane potentials of each neuron for all trials. c, Average auto-correlation of membrane potential of neurons within one trial. It has a frequency of 33 Hz, and the envelope of the auto-correlation decays as time lag increases. d, Average auto-correlation of excitatory synaptic inputs received by neurons within one trial. It has a frequency of 33 Hz, and the envelope of the auto-correlation decays as time lag increases. e, Average cross-correlation of membrane potential of pairs of nearby neurons within one trial. It also has a frequency of 33 Hz, and it decays with time lag but has a smaller amplitude (maximum cross-correlation is 0.5).

Finally, we consider the correlation properties of membrane potential. We calculate the auto-correlation of the membrane potential of each neuron and then calculate the average. As shown in Figure 6c, there are oscillations in auto-correlation as a function of time lag with a frequency of 33 Hz. Applying the same procedure to the excitatory synaptic inputs received by each neuron (Fig. 6d) yields a similar auto-correlation, with a characteristic frequency of 33 Hz. We also calculate the cross-correlation of the membrane potentials of nearby neurons. The averaged cross-correlation has the same frequency as the averaged auto-correlation but has a smaller amplitude of 0.5, as shown in Figure 6e; this is due to the dependence of excitatory connection strength upon distance (Eq. 4, see Materials and Methods). These correlated functions of membrane potential produced in our model are similar to those found in previous experiments of awake animals (Poulet and Petersen, 2008; Gentet et al., 2010). The frequency is, however, an order of magnitude greater than that found experimentally (Poulet and Petersen, 2008; Gentet et al., 2010), a discrepancy that is likely due to the scale of our network (see Materials and Methods). Suppose that the speed of the wave patterns is v and the spatial separation between them is l; then on average, the interval between waves approaching a certain area would be l/v. The corresponding frequency of synaptic inputs to the neurons around the area would therefore be v/l, meaning that increasing wave separations can potentially decrease the frequency. Indeed, when the network is scaled up to 600 × 600 with larger synaptic coupling ranges (D^E = 20 and D^I = 30 in Eq. 4), we find that the average frequency decreases from 33 Hz to 21 Hz.