Between Perfectly Critical and Fully Irregular: A Reverberating Model Captures and Predicts Cortical Spike Propagation

Abstract Knowledge about the collective dynamics of cortical spiking is very informative about the underlying coding principles. However, even most basic properties are not known with certainty, because their assessment is hampered by spatial subsampling, i.e., the limitation that only a tiny fraction of all neurons can be recorded simultaneously with millisecond precision. Building on a novel, subsampling-invariant estimator, we fit and carefully validate a minimal model for cortical spike propagation. The model interpolates between two prominent states: asynchronous and critical. We find neither of them in cortical spike recordings across various species, but instead identify a narrow “reverberating” regime. This approach enables us to predict yet unknown properties from very short recordings and for every circuit individually, including responses to minimal perturbations, intrinsic network timescales, and the strength of external input compared to recurrent activation “thereby informing about the underlying coding principles for each circuit, area, state and task.


Introduction
In order to understand how each cortical circuit or network processes its input, it would be desirable to first know its basic dynamical properties. For example, knowing which impact one additional spike has on the network (London et al. 2010) would give insight into the amplification of small stimuli Suarez et al. 1995;Miller 2016). Knowing how much of cortical activity can be attributed to external activation or internal activation (Reinhold et al. 2015) would allow to gauge how much of cortical activity is actually induced by stimuli, or rather internally generated, for example in the context of predictive coding (Rao and Ballard 1999;Clark 2013). Knowing the intrinsic network timescale (Murray et al. 2014) would inform how long stimuli are maintained in the activity and can be read out for short-term memory (Buonomano and Merzenich 1995;Wang 2002;Jaeger et al. 2007;Lim and Goldman 2013). However, not even these basic properties of cortical network dynamics are generally known with certainty.
In the past, insights about these network properties have been strongly hampered by the inevitable limitations of spatial subsampling, i.e., the fact that only a tiny fraction of all neurons can be recorded experimentally with millisecond precision. Such spatial subsampling fundamentally limits virtually any recording and hinders inferences about the collective response of cortical networks (Priesemann et al. 2009(Priesemann et al. , 2014Ribeiro et al. 2010Ribeiro et al. , 2014Levina and Priesemann 2017).
To describe network responses, two contradicting hypotheses have competed for more than a decade, and are the subjects of ongoing scientific debate: one hypothesis suggests that collective dynamics are "asynchronous-irregular" (AI) (Burns

Material and Methods
We analyzed in vivo spiking activity from Macaque monkey prefrontal cortex during a short-term memory task (Pipa et al. 2009;Franke et al. 2010;Pröpper and Obermayer 2013), anesthetized cat visual cortex with no stimulus (Blanche and Swindale 2006;Blanche 2009), and rat hippocampus during a foraging task (Mizuseki et al. 2009a(Mizuseki et al. , 2009b) (Supp. 1). We compared the recordings of each experimental session to results of a minimal model of spike propagation, which is detailed in the following.

Minimal Model of Spike Propagation
To gain an intuitive understanding of our mathematical approach, make a thought experiment in your favorite spiking network: apply one additional spike to an excitatory neuron, in analogy to the approach by (London et al. 2010). How does the network respond to that perturbation? As a first order approximation, one quantifies the number of spikes that are directly triggered additionally in all postsynaptic neurons. This number may vary from trial to trial, depending on the membrane potential of the postsynaptic neurons. However, what interests us most is m, the mean number of spikes triggered by the one extra spike. Any of these triggered spikes can in turn trigger spikes in their postsynaptic neurons in a similar manner, and thereby the perturbation may cascade through the system.
In the next step, assume that perturbations are started continuously at rate h, for example through afferent input from other brain areas or sensory modalities. Together, this leads to the mathematical framework of a branching model (Harris 1963;Heathcote 1965;Pakes 1971;Beggs and Plenz 2003;Haldeman and Beggs 2005;Ribeiro et al. 2010;Priesemann et al. 2013Priesemann et al. , 2014. This framework describes the number of active neurons A t in discrete time bins of length Δt. Here, Δt should reflect the propagation time of spikes between neurons. Formally, each spike i at the time bin t excites a random number Y t i , of postsynaptic spikes, on average = m Y t i , . The activity + A t 1 , i.e., the total number of spikes in the next time bin is then defined as the sum of the postsynaptic spikes of all current spikes A t , as well as the input h t : , t This generic spiking model can generate dynamics spanning AI and critical states depending on the input , and hence is well suited to probe network dynamics in vivo (see Supp. 3 for details). Most importantly, this framework enables us to infer m and other properties from the ongoing activity proper. Mathematical approaches to infer m are long known if the full network is sampled (Heyde and Seneta 1972;Wei 1991). Under subsampling, however, it is the novel estimator described in  that for the first time allows an unbiased inference of m, even if only a tiny fraction of neurons is sampled.
A precise estimate of m is essential, because the dynamics of the model is mainly governed by m (Fig. 1a). Therefore, after inferring m, a number of quantities can be analytically derived, and others can be obtained by simulating a branching model, which is constrained by the experimentally measured m and the spike rate.

Simulation
We simulated a branching model by mapping a branching process (Eq. (1) and Supp. 3) onto a random network of = N 10,000 neurons in the annealed disorder limit (Haldeman and Beggs 2005). An active neuron activated each of its κ = 4 postsynaptic neurons with probability κ = p m/ . Here, the activated postsynaptic neurons were drawn randomly without replacement at each step, thereby avoiding that two different active neurons would both activate the same target neuron. The branching model is critical for = m 1 in the infinite-size limit, and subcritical (supercritical) for < m 1 ( > m 1). We modeled input to the network at rate h by Poisson activation of each neuron at rate h N / . Subsampling (Priesemann et al. 2009) was applied to the model by sampling the activity of n neurons only, which were selected randomly before the simulation, and neglecting the activity of all other neurons. Thereby, instead of the full activity A t , only the subsampled activity a t was considered for observation.
If not stated otherwise, simulations were run for = L 10 7 time steps (corresponding to~11 h). Confidence intervals were estimated according to  from = B 100 realizations of the model, both for simulation and experiments.
We compared the experimental recordings to three different models: AI, near-critical, and reverberating. All three models were set up to match the experiment in the number of sampled neurons n and firing rate = 〈 〉 ( ⋅Δ ) R a n t / t . The AI and nearcritical models were set up with branching ratios of = m 0 or = m 0.9999, respectively. In addition, the reverberating model matched the recording in =m m, wherem was estimated from the recording using the novel subsampling-invariant estimator (see below). For all models, we chose a full network size of = N 10 4 and then determined the appropriate input in order to match the experimental firing rate. Exemplarily for the cat recording, which happened to represent the medianm, this yieldedˆ= m 0.98, = n 50, and = R Hz 7.25 . From these numbers, = h 290, = h 5.8 and = h 0.029 followed for the AI, reverberating, and near-critical models, respectively.
In Fig. 2, the reverberating branching model was also matched to the length of the cat recording of 295 s. To test for stationarity, the cat recording and the reverberating branching model were split into 59 windows of 5 s each, before estimating m for each window. In Fig. 1c, subcritical and critical branching models with = N 10 4 and = A 100 t were simulated, and = n 100 units sampled.

Subsampling-Invariant Estimation ofm
Details on the analysis are found in Supp. 2. For each experimental recording, we collected the spike times of all recorded units (both single and multi-units) into one single train of , reflecting the propagation time of spikes from one neuron to the next.
From these experimental time series, we estimatedm using the multistep regression (MR) estimator described in all detail in . In brief, we calculated the linear regression slope Δ r k t , which describes the linear statistical dependence of + a t k upon a t , for different time lags δ = Δ t k t with = … k k 1, , max . In our branching model, these slopes are expected to follow the relation , where b is an unknown parameter that depends on the higher moments of the underlying process and the degree of subsampling. However, it can be partialled out, allowing for an estimation of m without further knowledge about b. Throughout this study, we chose = k 2500 max (corresponding to 10 s) for the rat recordings, = k 150 max (600 ms) for the cat recording, and k max = 500 (2000 ms) for the monkey recordings, assuring that Δ k t max was always in the order of multiple intrinsic network timescales. In order to test for the applicability of a MR estimation, we used a set of conservative tests . The exponential relation can be rewritten as an exponential autocorrelation function . While the precise value of m depends on the choice of the bin size Δt and should only be interpreted together with the bin size (Δ = t m s 4 throughout this study), the intrinsic network timescale is independent of Δt. Therefore, we report both values in the following.

Reverberating Spiking Activity In Vivo
We applied MR estimation to the binned population spike counts a t of the recorded neurons of each experimental session across three different species (see Methods). We identified a limited range of branching ratios in vivo: in the experimentsm ranged from 0.963 to 0.998 (median = 0.98 , for a bin size of ), which is only a narrow window in the continuum from AI ( = m 0) to critical ( = m 1). For these values of found in cortical networks, the corresponding τ are between 100 ms and 2 s (median 247 ms, Figs 1e and S1). This clearly suggests that spiking activity in vivo is neither AI-like, nor consistent with a critical state. Instead, it is poised in a regime that, unlike critical or AI, does not maximize one particular property alone but may flexibly combine features of both . Without a prominent characterizing feature, we name it the reverberating regime, stressing that activity reverberates (different from the AI state) at timescales of hundreds of milliseconds (different from a critical state, where they can persist infinitely).

Validity of the Approach
There is a straight-forward verification of the validity of our phenomenological model: it predicts an exponential autocorrelation function δ r t for the population activity a t . We found that the activity in cat visual cortex (Fig. 2a,a') is surprisingly well described by this exponential fit (Fig. 2b,b'). This validation holds to the majority of experiments investigated (14 out of 21, Fig. S1).
A second verification of our approach is based on its expected invariance under subsampling: We further subsampled the activity in cat visual cortex by only taking into account spikes recorded from a subset ′ n out of all available n single units. As predicted (Fig. 2c), the estimates ofm, or equivalently of the intrinsic network timescale τ, coincided for any subset of single units if at least about five of the available 50 single units were evaluated (Fig. 2c'). Deviations when evaluating only a small subset of units most likely reflect the heterogeneity within cortical networks. Together, these results demonstrate that our approach returns consistent results when evaluating the activity of ≥ n 5 neurons, which were available for all investigated experiments.

Origin of the Activity Fluctuations
The fluctuations found in cortical spiking activity, instead of being intrinsically generated, could in principle arise from nonstationary input, which could in turn lead to misestimation of m (Priesemann and Shriki 2018). This is unlikely for three reasons: first, the majority of experiments passed a set of conservative tests that reject recordings that show any signature of common non-stationarities, as defined in . Second, recordings in cat visual cortex were acquired in absence of any stimulation, excluding stimulusrelated non-stationarities. Third, when splitting the spike recording into short windows, the window-to-window variation ofm in the recording did not differ from that of stationary in vivo-like reverberating models ( = p 0.3, Fig. 2d,d'). For these reasons, the observed fluctuations in the estimates likely originate from the characteristic fluctuations of collective network dynamics within the reverberating regime.

Timescales of the Network and Single Units
The dynamical state described by m directly relates to an exponential autocorrelation function with an intrinsic network timescale τ = −Δt m / ln . Exemplarily for the cat recording, = m 0.98 implies an intrinsic network timescale of τ = 188 ms, with Δ = t 4 ms reflecting the spike propagation time from one neuron to the next. While the autocorrelation function of the full network activity is expected to show an exponential decay (Fig. 3a, blue), this is different for the autocorrelation of single neurons-the most extreme case of subsampling. We showed that subsampling can strongly decrease the absolute values of the autocorrelation function for any non-zero time lag (Fig. 3a, gray). This effect is typically interpreted as a lack of memory, because the autocorrelation of single neurons decays at the order of the bin size (Fig. 3a, red). However, ignoring the value at δ = t 0, the floor of the autocorrelation function still unveils the exponential relation. Remarkably, the autocorrelation function of single units in cat visual cortex displayed precisely the shape predicted under subsampling (compare Fig. 3a,b).

Established Methods are Biased to Identifying AI Dynamics
On the population level, networks with different m are clearly distinguishable (Fig. 1a). Surprisingly, single-neuron statistics, namely inter-spike interval (ISI) distributions, Fano factors, conventional estimation of m, and the autocorrelation strength δ r t , all returned signatures of AI activity regardless of the underlying network dynamics, and hence these single-neuron properties do not serve as a reliable indicator for the network's dynamical state.
First, exponential ISI distributions are considered a strong indicator of Poisson-like firing. Surprisingly, the ISIs of single neurons in the in vivo-like branching model closely followed exponential distributions as well. The ISI distributions were almost indistinguishable for reverberating and AI models (Figs 4a,a' and S2). In fact, the ISI distributions are mainly determined by the mean firing rate. This result was further supported by coefficients of variation close to unity, as expected for exponential ISI distributions and Poisson firing (Fig. S2).
Second, for both the AI and reverberating regime, the Fano factor F for single unit activity was close to unity, a hallmark feature of irregular spiking (Tolhurst et al. 1981;Vogels et al. 1989;Softky and Koch 1993;de Ruyter van Steveninck et al. 1997;Gur et al. 1997;Kara et al. 2000;Carandini 2004) (Fig. 5g, analytical result: Eq. (S9)). Hence it cannot serve to distinguish between these different dynamical states. When evaluating more units, or increasing the bin size to 4 s, the differences became more pronounced, but for experiments, the median Fano factor of single unit activity did not exceed = F 10 in any of the experiments, even in those with the longest reverberation (Figs 4b,b' and S3). In contrast, for the full network the Fano factor rose to ≈ F 10 4 for the in vivo-like branching model and diverged when approaching criticality (Fig. 5g, analytical result: Eq. (S5)).
Third, conventional regression estimators (Heyde and Seneta 1972;Wei 1991) are biased towards inferring irregular activity, as shown before. Here, conventional estimation yielded a median ofˆ= m 0.057 for single-neuron activity in cat  at lag δ = t 0 to the next lag ±Δ r t (gray solid line). We showed previously that this drop is a subsampling-induced bias. When ignoring the zero-lag value, the autocorrelation strength is decreased, but the exponential decay and even the value of the intrinsic network timescale τ of the network activity are preserved (inset). The red, dashed line shows a potential, naive exponential function, fitted to the single-neuron autocorrelation function (gray). This naive fit would return a much smaller τ. (b) The autocorrelation function of single-neuron activity recorded in cat visual cortex (gray) precisely resembles this theoretical prediction, namely a sharp drop and then an exponential decay (blue, inset), which persists over more than 100 ms. A naive exponential fit (red) to the single-neuron data would return an extremely short τ. visual cortex, in contrast toˆ= m 0.954 returned by MR estimation (Fig. S9).
Fourth, for the autocorrelation function of an experimental recording (Fig. 3b) the rapid decay of δ r t prevails, and hence single-neuron activity appears uncorrelated in time.

Cross-validation of Model Predictions
We compared the experimental results to an in vivo-like model, which was matched to each experiment only in the average firing rate, and in the inferred branching ratiom. Remarkably, this in vivo-like branching model could predict statistical properties not only of single neurons (ISI and Fano factor, see above), but also pairwise and population properties, as detailed below. This prediction capability further underlines the usefulness of this simple model to approximate the default state of cortical dynamics. First, the model predicted the activity distributions, ( ) p a t , better than AI or critical models for the majority of experiments (15 out 21, Figs 4c,d,c',d', S5, and S6), both for the exemplary bin sizes of 4 ms and 40 ms. Hence, the branching models, which were only matched in their respective first moment of the activity distributions (through the rate) and first moment of the spreading behavior (through m), in fact approximated all higher moments of the activity distributions ( ) p a t . Likewise, the model predicted the distributions of neural avalanches, i.e., spatio-temporal clusters of activity (Figs 4e,f,e', f', S7, and S8). Characterizing these distributions is a classic approach to assess criticality in neuroscience (Beggs and Plenz 2003;Priesemann et al. 2014), because avalanche size and duration distributions, ( ) p s and ( ) p d , respectively, follow power laws in critical systems. In contrast, for AI activity, they are approximately exponential (Priesemann and Shriki 2018). The matched branching models predicted neither exponential nor power law distributions for the avalanches, but very well matched the experimentally obtained distributions (compare red and blue in Figs 4e,f,e',f', S7, and S8). Indeed, model likelihood (Clauset et al. 2009) favored the in vivo-like branching model over Poisson and critical models for the majority experiments (18 out of 21, Fig. S7). Our results here are consistent with those of spiking activity in awake animals, which typically do not display power laws (Bédard et al. 2006;Ribeiro et al. 2010;Priesemann et al. 2014). In contrast, most evidence for criticality in vivo, in particular the characteristic power-law distributions, has been obtained from coarse measures of neural activity (LFP, EEG, BOLD; see Priesemann et al. 2014 and references therein). Last, the model predicted the pairwise spike count crosscorrelation r sc . In experiments, r sc is typically between 0.01 and 0.25, depending on brain area, task, and most importantly, the analysis timescale (bin size) (Cohen and Kohn 2011). For the cat recording the model even correctly predicted the bin size dependence of r sc from¯≈ r 0.004 sc at a bin size of 4 ms (analytical result: Eq. (S12)) to¯≈ r 0.3 sc at a bin size of 2 s (Fig. 4g). Comparable results were also obtained for some monkey experiments. In contrast, correlations in most monkey and rat recordings were smaller than predicted (Figs 4g' and S4). It is very surprising that the model correctly predicted the crosscorrelation even in some experiments, as m was inferred only from the temporal structure of the spiking activity alone, whereas r sc characterizes spatial dependencies.
Overall, by only estimating the effective synaptic strength m from the in vivo recordings, higher-order properties like avalanche size distributions, activity distributions and in some cases spike count cross-correlations could be closely matched using the generic branching model.

The Dynamical State Determines Responses to Small Stimuli
After validating the model using a set of statistical properties that are well accessible experimentally, we now turn to making predictions about yet unknown properties, namely network responses to small stimuli. In the line of London et al. (2010), assume that on a background of spiking activity one single extra spike is triggered. This spike may in turn trigger new spikes, leading to a cascade of additional spikes Δ t propagating through the network. A dynamical state with branching ratio m implies that on average, this perturbation decays with time constant τ = −Δt m / log . Similar to the approach in London et al. (2010), the evolution of the mean firing rate, averaged over a reasonable number of trials (here: 500) unveils the nature of the underlying spike propagation: depending on m, the rate excursions will last longer, the higher m (Figs 5a-c and S11). The perturbations are not deterministic, but show trial-to-trial variability which also increases with m (S11b).
Unless > m 1, the theory of branching models ensures that perturbations will die out eventually after a duration Δ d , having accumulated a total of = ∑ Δ Δ = s t d t 1 extra spikes in total. This perturbation size Δ s and duration Δ d follow specific distributions (Harris 1963), which are determined by m: they are power law distributed in the critical state ( = m 1), with a cutoff for any < m 1 (Figs 5f and S11c,d). These distributions imply a characteristic mean perturbation size Δ s (Fig. 5d), which diverges at the critical point. The variability of the perturbation sizes is also determined by m and also diverges at the critical point (inset of Figs 5d and S11e).
Taken together, these results imply that the closer a neuronal network is to criticality, the more sensitive it is to external perturbations, and the better it can amplify small stimuli. At the same time, these networks also show larger trial-to-trial variability. For typical cortical networks, we found that the response to one single extra spike will on average comprise between 20 and 1000 additional spikes in total (Fig. 5e).

The Dynamical State Determines Network Susceptibility and Variability
Moving beyond single spike perturbations, our model gives precise predictions for the network response to continuous stimuli. If extra action potentials are triggered at rate h in the network, the network will again amplify these external activations, depending on m. Provided an appropriate stimulation protocol, this rate response could be measured and our prediction tested in experiments (Fig. S11g). The susceptibility ∂ ∂ R h / diverges at the critical transition and is unique to a specific branching ratio m. We predict that typical cortical networks will amplify a small, but continuous input rate by about a factor fifty (Fig. S11h, red).
While the input and susceptibility determine the network's mean activity, the network still shows strong rate fluctuations around this mean value. The magnitude of these fluctuations in relation to the mean can be quantified by the network Fano factor = [ ] 〈 〉 F A A Var / t t (Fig. 5g). This quantity cannot be directly inferred from experimental recordings, because the Fano factor of subsampled populations severely underestimates the network Fano factor, as shown before. We here used our in vivolike model to obtain estimates of the network Fano factor: for a bin size of Δ = t 4 ms it is about ≈ F 40 and rises to ≈ F 4000 for bin sizes of several seconds, highlighting that network fluctuations probably are much stronger than one would naively assume from experimental, subsampled spiking activity.
Distinguishing Afferent and Recurrent Activation Last, our model gives an easily accessible approach to solving the following question: given a spiking neuronal network, which fraction of the activity 〈 〉 A is generated by recurrent activation from within the network, and which fraction can be attributed to external, afferent excitation h? The branching model readily provides an answer: the fraction of external activation is 〈 〉 = − h A m / 1 (Fig. 5h). In this framework, AI-like networks are completely driven by external input currents or noise, whereas reverberating networks generate a substantial fraction of their activity intrinsically. For the experiments investigated in this study, we inferred that between 0.1% and 7% of the activity are externally generated (median 2%, Fig. 5i).
While our model is quite simplistic given the complexity of neuronal network activity, keep in mind that "all models are wrong, but some are useful" (Box 1979). Here, the model has proven to provide a good first order approximation to a number of statistical properties of spiking activity and propagation in cortex. Hence, it promises insight into cortical function because (i) it relies on simply assessing spontaneous cortical activity, (ii) it does not require manipulation of cortex, (iii) it enables reasonable predictions about sensitivity, amplification, and internal and external activation, (iv) this analysis is possible in an area specific, task-and state-dependent manner as only short recordings are required for consistent results.

Our Results Resolve Contradictions Between AI and Critical States
Our results for spiking activity in vivo suggest that network dynamics show AI-like statistics, because under subsampling the observed correlations are underestimated. In contrast, typical experiments that assessed criticality potentially overestimated correlations by sampling from overlapping populations (LFP, EEG) and thereby hampered a fine distinction between critical and subcritical states (Pinheiro Neto and Priesemann, in preparation). By employing for the first time a consistent, quantitative estimation, we provided evidence that in vivo spiking population dynamics reflects a reverberating regime, i.e., it operates in a narrow regime around = m 0.98. This result is supported by the findings by Dahmen et al. (2016): based on distributions of covariances, they inferred that cortical networks operate in a regime below criticality. Given the generality of our results across different species, brain areas, and cognitive states, our results suggest self-organization to this reverberating regime as a general organization principle for cortical network dynamics.

The Reverberating Regime Combines Features of AI and Critical State
At first sight,ˆ= m 0.98 of the reverberating regime may suggest that the collective spiking dynamics is very close to critical. Indeed, physiologically a Δ ≈ m 1.6% difference to criticality ( = m 1) is small in terms of the effective synaptic strength. However, this apparently small difference in single unit properties has a large impact on the collective dynamical fingerprint and makes AI, reverberating, and critical states clearly distinct: for example, consider the sensitivity to a small input, i.e., the susceptibility χ = ∂ ∂ = − R h / m 1 1 . The susceptibility diverges at criticality, making critical networks overly sensitive to input. In contrast, states with ≈ m 0.98 assure sensitivity without instability. Because this has a strong impact on network dynamics and putative network function, finely distinguishing between dynamical states is both important and feasible even if the corresponding differences in effective synaptic strength (m) appear small.
We cannot ultimately rule out that cortical networks selforganize as close as possible towards criticality, the platonic ideal being impossible to achieve for example due to finite-size, external input, and refractory periods. Therefore, the reverberating regime might conform with quasi-criticality (Williams-García et al. 2014) or neutral theory (Martinello et al. 2017). However, we deem this unlikely for two reasons. First, in simulations of finite-size networks with external input, we could easily distinguish the reverberating regime from states with = m 0.9999 , which are more than one order of magnitude closer to criticality than any experiment we analyzed. Second, operating in a reverberating regime, which is between AI and critical, may combine the computational advantages of both states : the reverberating regime enables rapid changes of computational properties by small parameter changes, keeps a sufficient safety-margin from instability to make seizures sufficiently unlikely (Priesemann et al. 2014), balances competing requirements (e.g., sensitivity and specificity (Gollo 2017)), and may carry short-term memory and allow to integrate information over limited, tunable timescales (Wang 2002;Boedecker et al. 2012). For these reasons, we consider the reverberating regime to be the explicit target state of self-organization. This is in contrast to the view of "as close to critical as possible," which still holds criticality as the ideal target.

More Complex Network Models
Cortical dynamics is clearly more complicated than a simple branching model. For example, heterogeneity of single-neuron morphology and dynamics, and non-trivial network topology likely impact population dynamics. However, we showed that statistics of cortical network activity are well approximated by a branching model. Therefore, we interpret branching models as a statistical approximation of spike propagation, which can capture a fair extent of the complexity of cortical dynamics. By using branching models, we draw on the powerful advantage of analytical tractability, which allowed for basic insight into dynamics and stability of cortical networks.
In contrast to the branching model, doubly stochastic processes (i.e., spikes drawn from an inhomogeneous Poisson distribution) failed to reproduce many statistical features (Fig. S10). We conjecture that the key difference is that doubly stochastic processes propagate the underlying firing rate instead of the actual spike count. Thus, propagation of the actual number of spikes (as e.g., in branching or Hawkes processes; Kossio et al. 2018), not some underlying firing rate, seems to be integral to capture the statistics of cortical spiking dynamics.
Our statistical model stands in contrast to generative models, which generate spiking dynamics by physiologically inspired mechanisms. One particularly prominent example are networks with balanced excitation and inhibition Sompolinsky 1996, 1997;Brunel 2000), which became a standard model of neuronal networks (Hansel and van Vreeswijk 2012). A balance of excitation and inhibition is supported by experimental evidence (Okun and Lampl 2008). Our statistical model reproduces statistical properties of such networks if one assumes that the excitatory and inhibitory contributions can be described by an effective excitation. In turn, the results obtained from the well-understood estimator can guide the refinement of generative models. For example, we suggest that network models need to be extended beyond the asynchronous-irregular state (Brunel 2000) to incorporate the network reverberations observed in vivo. Possible candidate mechanisms are increased coupling strength or inhomogeneous connectivity. Both have already been shown to induce rate fluctuations with timescales of several hundred milliseconds (Litwin-kumar and Doiron 2012; Ostojic 2014; Kadmon and Sompolinsky 2015).
Because of the assumption of uncorrelated, Poisson-like network firing, models that study single neurons typically assume that synaptic currents are normally distributed. Our results suggest that one should rather use input with reverberating properties with timescales of a few hundred milliseconds to reflect input from cortical neurons in vivo. This could potentially change our understanding of single-neuron dynamics, for example of their input-output properties.

Deducing Network Properties from the Tractable Model
Using our analytically tractable model, we could predict and validate network properties, such as avalanche size and duration, ISI, or activity distributions. Given the experimental agreement with these predictions, we deduced further properties, which are impossible or difficult to assess experimentally and gave insight into more complex questions about network responses: How do perturbations propagate within the network, and How susceptible is the network to external stimulation?
One particular question we could address is the following: Which fraction of network activity is attributed to external or recurrent, internal activation? We inferred that about 98% of the activity is generated by recurrent excitation, and only about 2% originates from input or spontaneous threshold crossing. This result may depend systematically on the brain area and cognitive state investigated: For layer 4 of primary visual cortex in awake mice, Reinhold et al. (2015) concluded that the fraction of recurrent cortical excitation is about 72%, and cortical activity dies out with a timescale of about 12 ms after thalamic silencing. Their numbers agree perfectly well with our phenomenological model: a timescale of τ = 12 ms implies that the fraction of recurrent cortical excitation is , just as found experimentally. Under anesthesia, in contrast, they report timescales of several hundred milliseconds, in agreement with our results. These differences show that the fraction of external activation may strongly depend on cortical area, layer, and cognitive state. The novel estimator can in future contribute to a deeper insight into these differences, because it allows for a straight-forward assessment of afferent versus recurrent activation, simply from evaluating spontaneous spiking activity, without the requirement of thalamic or cortical silencing.