Subsampling scaling: a theory about inference from partly observed systems

TL;DR

This paper develops a theoretical framework to correct biases caused by spatial subsampling in observations of complex systems, enabling accurate inference of their true properties.

Contribution

The authors derive a subsampling scaling theory applicable to various observables, allowing correction of biases and analysis of system criticality from partly observed data.

Findings

Power-law scaling observed only in mature neural networks

Subsampling scaling distinguishes critical from subcritical systems

Framework applicable to diverse systems like epidemics and neural networks

Abstract

In real-world applications, observations are often constrained to a small fraction of a system. Such spatial subsampling can be caused by the inaccessibility or the sheer size of the system, and cannot be overcome by longer sampling. Spatial subsampling can strongly bias inferences about a system's aggregated properties. To overcome the bias, we derive analytically a subsampling scaling framework that is applicable to different observables, including distributions of neuronal avalanches, of number of people infected during an epidemic outbreak, and of node degrees. We demonstrate how to infer the correct distributions of the underlying full system, how to apply it to distinguish critical from subcritical systems, and how to disentangle subsampling and finite size effects. Lastly, we apply subsampling scaling to neuronal avalanche models and to recordings from developing neural networks. We show that only mature, but not young networks follow power-law scaling, indicating self-organization to criticality during development.

Equations

$$P_{\mathrm{sub}}(s)=P(X_{\mathrm{sub}}=s)=\sum _{k=0}^{\infty }P(s+k)\left(\genfrac{}{}{0ex}{}{s+k}{s}\right)p^s(1-p{)}^k.$$

$${\lambda }_{\mathrm{sub}}=\ln \left(\frac{e^{\lambda }+p-1}{p}\right)\iff \lambda =\ln ((e^{{\lambda }_{\mathrm{sub}}}-1)p+1).$$

$$C_{\lambda }/C_{{\lambda }_{\mathrm{sub}}}=1-e^{-\lambda }+p{e}^{-\lambda }=\frac{e^{-{\lambda }_{\mathrm{sub}}}+p-p{e}^{-{\lambda }_{\mathrm{sub}}}}{p}.$$

$$P(s)$$

$$=\left(1-e^{-\lambda }+p{e}^{-\lambda }\right)P_{\mathrm{sub}}\left(\frac{\lambda }{\ln \left(\frac{e^{\lambda }+p-1}{p}\right)}s\right)=G(p,\lambda )P_{\mathrm{sub}}(sF(p,\lambda )).$$

$$P(s)=p^a{P}_{\mathrm{sub}}(p^{b}s),\mathrm{for}\mathrm{any}a,b\in \text{mathds}R\mathrm{with}a-b\gamma =1-\gamma .$$

$$P(s)\approx p{P}_{\mathrm{sub}}(p\cdot s).$$

$$N{M}^{\beta -1}{P}_{\mathrm{sub}}(sN{M}^{\nu -1};M,N)=g(s).$$

$$\text{if }z(x,y,t)\ge 4:$$

$$z(x,y,t+1)=z(x,y,t)-4$$

$$z(x\pm 1,y\pm 1,t+1)=z(x\pm 1,y\pm 1,t)+1$$

$$P_{\mathrm{sub}}(s)=\sum _{k=s}^{\infty }P(k)\left(\genfrac{}{}{0ex}{}{k}{s}\right)p^s(1-p{)}^{k-s}.$$

$$G_{\mathrm{sub}}(z;p)$$

$$G_{\mathrm{sub}}(z;p)=G(1-p(1-z))$$

$$E(X_{\mathrm{sub}})=G_{\mathrm{sub}}^{\prime }(1^{-};p)=G^{\prime }(1-p(1-1^{-}))=p{G}^{\prime }(1^{-})=pE(X)$$

$$P(X=s)=\left(\genfrac{}{}{0ex}{}{s+r-1}{s}\right)p_{\mathrm{NB}}^r(1-p_{\mathrm{NB}}{)}^s,$$

$$m=E(X)=r\frac{(1-p_{\mathrm{NB}})}{p_{\mathrm{NB}}},$$

$$G(z)={\left(\frac{p_{\mathrm{NB}}z-z+1}{p_{\mathrm{NB}}}\right)}^{-r}={\left(1+\frac{m}{r}(1-z)\right)}^{-r}.$$

$$G_{\mathrm{sub}}(z;p)={\left(1+\frac{pm}{r}(1-z)\right)}^{-r},$$

$$\frac{1-p_{\mathrm{NB}}^{\prime }}{p_{\mathrm{NB}}^{\prime }}=p\frac{1-p_{\mathrm{NB}}}{p_{\mathrm{NB}}}.$$

$$P_{\exp }(s)=(1-e^{-\lambda })\cdot e^{-\lambda s},\mathrm{with}s\in \text{mathds}{N}_0.$$

$$\frac{e^{-{\lambda }_{\mathrm{sub}}}}{1-e^{-{\lambda }_{\mathrm{sub}}}}=p\frac{e^{-\lambda }}{1-e^{-\lambda }}\iff e^{{\lambda }_{\mathrm{sub}}}-1=\frac{e^{\lambda }-1}{p}.$$

$${\lambda }_{\mathrm{sub}}=\ln \left(\frac{e^{\lambda }+p-1}{p}\right).$$

$$P_{\mathrm{sub}}(s;\gamma ,p)=C_{\gamma }\sum _{n=0}^{\infty }(s+n{)}^{-\gamma }{p}^s(1-p{)}^n\left(\genfrac{}{}{0ex}{}{s+n}{n}\right).$$

$$P_{\mathrm{sub}}(s;\gamma ,p)\stackrel{s\to \infty }{⟶}{c}_{\gamma }(p)s^{-\gamma }.$$

$$\frac{\partial }{\partial p}{P}_{\mathrm{sub}}(s;\gamma ,p)\stackrel{s\to \infty }{⟶}{s}^{-\gamma }\frac{\partial }{\partial p}{c}_{\gamma }(p).$$

$$\frac{\partial }{\partial p}{P}_{\mathrm{sub}}(s;\gamma ,p)$$

$$C_{\gamma }s\sum _{n=0}^{\infty }(s+n{)}^{-\gamma }{p}^{s-1}(1-p{)}^n\left(\genfrac{}{}{0ex}{}{s+n}{n}\right)$$

$$C_{\gamma }\sum _{n=0}^{\infty }(s+n{)}^{-\gamma }{p}^{s}n(1-p{)}^{n-1}\left(\genfrac{}{}{0ex}{}{s+n}{n}\right)$$

$$s^{-\gamma }\frac{\partial }{\partial p}{c}_{\gamma }(p)\approx \frac{\partial }{\partial p}{P}_{\mathrm{sub}}(s;\gamma ,p)\approx \frac{s}{p}{c}_{\gamma }(p)s^{-\gamma }-\frac{s+1}{p}{c}_{\gamma }(p)(s+1{)}^{-\gamma }.$$

$$\frac{\partial }{\partial p}{c}_{\gamma }(p)=\frac{c_{\gamma }(p)}{p}\underset{s\to \infty }{\lim }\left[s-(s+1){\left(1+\frac{1}{s}\right)}^{-\gamma }\right]$$

$$\underset{x\to 0}{\lim }\frac{(1+x{)}^{\mu }-1}{\mu x}=1,$$

Full text

Subsampling scaling: a theory about inference from partly observed systems

A. Levina1,3∗ and V. Priesemann2,3∗

1Institute of Science and Technology Austria, Am Campus 1, 3400 Klosterneuburg

2Max Planck Institute for Dynamics and Self-Organization and

3Bernstein Center for Computational Neuroscience, Am Fassberg 17, 37077 Göttingen, Germany

∗anna/viola@nld.ds.mpg.de; the authors contributed equally

Abstract

In real-world applications, observations are often constrained to a small fraction of a system. Such spatial subsampling can be caused by the inaccessibility or the sheer size of the system, and cannot be overcome by longer sampling. Spatial subsampling can strongly bias inferences about a system’s aggregated properties. To overcome the bias, we derive analytically a subsampling scaling framework that is applicable to different observables, including distributions of neuronal avalanches, of number of people infected during an epidemic outbreak, and of node degrees. We demonstrate how to infer the correct distributions of the underlying full system, how to apply it to distinguish critical from subcritical systems, and how to disentangle subsampling and finite size effects. Lastly, we apply subsampling scaling to neuronal avalanche models and to recordings from developing neural networks. We show that only mature, but not young networks follow power-law scaling, indicating self-organization to criticality during development.

1 Introduction

Inferring global properties of a system from observations is a challenge, even if one can observe the whole system. The same task becomes even more challenging if one can only sample a small number of units at a time (spatial subsampling). For example, when recording spiking activity from a brain area with current technology, only a very small fraction of all neurons can be accessed with millisecond precision. To still infer global properties, it is necessary to extrapolate from this small sampled fraction to the full system.

Spatial subsampling affects inferences not only in neuroscience, but in many different systems: In disease outbreaks, typically a fraction of cases remains unreported, hindering a correct inference about the true disease impact [1, 2]. Likewise, in gene regulatory networks, typically a fraction of connections remains unknown. Similarly, when evaluating social networks, the data sets are often so large that because of computational constraints only a subset is stored and analyzed. Obviously, subsampling does not affect our inferences about properties of a single observed unit, such as the firing rate of a neuron. However, we are often confronted with strong biases when assessing aggregated properties, such as distributions of node degrees, or the number of events in a time window  [3, 4, 5, 6]. Concrete examples are distributions of the number of diseased people in an outbreak, the size of an avalanche in critical systems, the number of synchronously active neurons, or the number of connections of a node. Despite the clear difference between these observables, the mathematical structure of the subsampling problem is the same. Hence our novel inference approach applies to all of them.

Examples of subsampling biases, some of them dramatic, have already been demonstrated in numerical studies. For example, subsampling of avalanches in a critical model can make a simple monotonic distribution appear multi-modal [4]. In general, subsampling has been shown to affect avalanche distributions in various ways, which can make a critical system appear sub- or supercritical [5, 7, 8, 9, 10, 11], and sampling from a locally connected network can make the network appear “small-world” [6]. For the topology of networks, it has been derived that, contrary to common intuition, a subsample from a scale-free network is not itself scale-free [3]. Importantly, these biases are not due to limited statistics (which could be overcome by collecting more data, e.g. acquiring longer recordings, or more independent subsamples of a system), but genuinely originates from observing a small fraction of the system, and then making inferences including unobserved parts. Although subsampling effects are known, in the literature there is so far no general analytical understanding of how to overcome them. For subsampling effects on degree distributions, Stumpf and colleagues provided a first analytical framework, stating the problem of subsampling bias [3].

In this paper, we show how to overcome subsampling effects. To this end we develop a mathematical theory that allows to understand and revert them in a general manner. We validate the analytical approach using various simulated models, and finally apply it to infer distributions of avalanches in developing neural networks that are heavily subsampled due to experimental constraints. Finally, we show that finite-size and subsampling effects clearly differ, and derived a combined subsampling-finite-size scaling relation. Together, our results introduce a novel approach to study under-observed systems.

2 Results

2.1 Mathematical subsampling

To derive how spatial subsampling affects a probability distribution of observables, we define a minimal model of “mathematical subsampling”.

We first introduce the variables with the example of avalanches, which are defined as cascades of activity propagating on a network [12, 13], and then present the mathematical definition. The main object of interest is a “cluster”, e.g. an avalanche. The cluster size $s$ is the total number of events or spikes. In general, the cluster size is described by a discrete, non-negative random variable $X$. Let $X$ be distributed according to a probability distribution $P(X=s)=P(s)$. For subsampling, we assume for each cluster that each of its events is independently observed with probability $p$ (or missed with probability $1-p$). Then $X_{\mathrm{sub}}$ is a random variable denoting the number of observed events of a cluster, and $X-X_{\mathrm{sub}}$ the number of missed events. For neural avalanches, this subsampling is approximated by sampling a random fraction of all neurons. Then $X_{\mathrm{sub}}$ represents the number of all events generated by the observed neurons within one avalanche on the full system. Note, that this definition translates one cluster in the full system to exactly one cluster under subsampling (potentially of size zero; this definition does not require explicit binning, see Sec. 2.4 and Sec. 4.2.4). We call the probability distribution of $X_{\mathrm{sub}}$ “subsampled distribution” $P_{\mathrm{sub}}(s)$. An analogous treatment can be applied to e.g. graphs. There a “cluster” represents the set of (directed) connections of a specific node, and thus $X$ is the degree of that node. Under subsampling, i.e. considering a random subnetwork, only connections between observed nodes are taken into account, resulting in the subsampled degree $X_{\mathrm{sub}}$.

As each event is observed independently, the probability of $X_{\mathrm{sub}}=s$ is the sum over probabilities of observing clusters of $X=s+k$ events, where $k$ denotes the missed events and $s$ the sampled ones (binomial sampling):

[TABLE]

This equation holds for any discrete $P(s)$ defined on $\mathbb{N}_{0}$, the set of non-negative integers. To infer $P(s)$ from $P_{\mathrm{sub}}(s)$, we develop in the following a novel “subsampling scaling” that allows to parcel out the changes in $P(s)$ originating from spatial subsampling. A correct scaling ansatz collapses the $P_{\mathrm{sub}}(s)$ for any sampling probability $p$.

In the following, we focus on subsampling from two specific families of distributions that are of particular importance in the context of neuroscience, namely exponential distributions $P(s)=C_{\lambda}e^{-\lambda s}$ with $\lambda>0$, and power laws $P(s)=C_{\gamma}s^{-\gamma}$ with $\gamma>1$. These two families are known to show different behaviors under subsampling [3]:

For exponential distributions, $P(s)$ and $P_{\mathrm{sub}}(s)$ belong to the same class of distributions, only their parameters change under subsampling. Notably, this result generalizes to positive and negative binomial distributions, which include Poisson distributions. 2. 2.

Power-laws or scale-free distributions, despite their name, are not invariant under subsampling. Namely, if $P(s)$ follows a power-law distribution, then $P_{\mathrm{sub}}(s)$ is not a power law but only approaching it in the limit of large cluster size ($s\rightarrow\infty$).

In more detail, for exponential distributions, $P(s)=C_{\lambda}e^{-\lambda s}$, $s\in\mathbb{N}_{0}$, subsampling with probability $p$ results in an exponential distribution with decay parameter $\lambda_{\mathrm{sub}}$ that can be expressed as a function of $\lambda$ and $p$ (for the full analytical derivation see Sec. 5.1.1):

[TABLE]

Likewise, changes in the normalizing constant $C_{\lambda}=1-e^{-\lambda}$ of $P(s)$ are given by:

[TABLE]

These two relations allow to derive explicitly a subsampling scaling for exponentials, i.e. the relation between $P(s)$ and $P_{\mathrm{sub}}(s)$:

[TABLE]

Thus given an exponential distribution $P(s)$ of the full system, all distributions under subsampling can be derived. Vice versa, given the observed subsampled distribution $P_{\mathrm{sub}}(s)$, the full distribution can be analytically derived if the sampling probability $p$ is known. Therefore, for exponentials, the scaling ansatz above allows to collapse all distributions obtained under subsampling with any $p$ (Fig. 1 A,B).

The presented formalism is analogous to the one proposed by Stumpf et al. [3]. They studied which distributions changed and which preserved their classes under subsampling. In the following we extend that study, and then develop a formalism that allows to extrapolate the original distribution from the subsampling, also in the case where an exact solution is not possible.

For power-law distributions of $X$, $X_{\mathrm{sub}}$ is not power-law distributed, but only approaches a power law in the tail ($s\to\infty$). An approximate scaling relation, however, collapses the tails of distributions as follows (mathematical derivation in Sec. 5.1.2). For $s\to\infty$, a power law $P(s)=C_{\gamma}s^{-\gamma}$ and the distributions obtained under subsampling can be collapsed by:

[TABLE]

For any $a,b$ satisfying the relation above, this scaling collapses the tails of power-law distributions. The “heads”, however, deviate from the power law and hence cannot be collapsed (see deviations at small $s$, Fig. 1 D). These deviations decrease with increasing $p$, and with $\gamma\rightarrow 1^{+}$ [3], (5.1.3). We call these deviations “hairs” because they “grow” on the “heads” of the distribution as opposed to the tails of the distribution. In fact, the hairs allow to infer the system size from knowing the number of sampled units if the full systems exhibits a power-law distribution (5.2).

In real world systems and in simulations, distributions often deviate from pure exponentials or pure power laws [14, 15]. We here treat the case that is typical for finite size critical systems, namely a power law that transits smoothly to an exponential around $s=s^{\mathrm{cutoff}}$ (e.g. Fig. 2 A). Under subsampling, $s_{\mathrm{sub}}^{\mathrm{cutoff}}$ depends linearly on the sampling probability: $s_{\mathrm{sub}}^{\mathrm{cutoff}}=p\cdot s^{\mathrm{cutoff}}$. Hence, the only solution to the power-law scaling relation (Eq. 5) that collapses (to the best possible degree), both, the power-law part of distributions and the onsets of the cutoff is the one with $a=b=1$:

[TABLE]

As this scaling is linear in $p$, we call it p-scaling. A priori, p-scaling is different from the scaling for exponentials (Eq. 4). However, p-scaling is a limit case of the scaling for exponentials under the assumption that $\lambda\ll p$: Taylor expansion around $\lambda=0$ results in the scaling relation $P(s)\approx pP_{\mathrm{sub}}(p\cdot s)$, i.e. the same as derived in Eq. 6. Indeed, for exponentials with $\lambda=0.001$ p-scaling results in a nearly perfect collapse for all $p>0.01$, however $p\leq 0.01$ violates the $\lambda\ll p$ requirement and the collapse breaks down (Fig. 1 B, inset). Thus p-scaling collapses power laws with exponential tail if $\lambda$ is small, and also much smaller than the sampling probability. This condition is typically met in critical, but not in subcritical systems (Sec. 5.3).

2.2 Subsampling in critical models

Experimental conditions typically differ from the idealized, mathematical formulation of subsampling derived above: Distributions do not follow perfect power laws or exponentials, and sampling is not necessarily binomial, but restricted to a fixed set of units. To mimic experimental conditions, we simulated avalanche generating models with fixed sampling in a critical state, because at criticality, subsampling effects are expected to be particularly strong: In critical systems, avalanches or clusters of activated units can span the entire system and thus under subsampling they cannot be fully assessed.

We simulated critical models with different exponents of $P(s)$ to assess the generality of our analytically derived results. The first model is the widely used branching model (BM) [8, 16, 17, 18, 19], and the second model is the Bak-Tang-Wiesenfeld model (BTW) [13], both studied in two variants. Both models display avalanches of activity after one random unit (neuron) has been activated externally (drive). In the BM, activity propagates stochastically, i.e. an active neuron activates any of the other neurons with a probability $p_{\mathrm{act}}$. Here $p_{\mathrm{act}}$ is the control parameter, and the model is critical in the infinite size limit if one spike on average triggers one spike in its postsynaptic pool (see Methods). We simulated the BM on a fully connected network and on a sparsely connected network. The avalanche size distributions of both BM variants have an exponent $\approx 1.5$ [17], and for both variants, subsampling results are very similar (Sec. 5.4). Hence in the main text we show results for the fully connected BM, while the results for the sparsely connected BM are displayed, together with results of a third model, the non-conservative model from Eurich, Herrmann & Ernst (EHE-model) [20], in Sec. 5.4. As expected, distributions of all critical models collapse under p-scaling.

In the BTW, activity propagates deterministically via nearest neighbors connections. Propagation rules reflect a typical neural non-leaky “integrate-and-fire” mechanism: Every neuron sums (integrates) its past input until reaching a threshold, then becomes active itself and is reset. The BTW was implemented classically with nearest neighbor connections on a 2D grid of size $M=L\times L$ either with open (BTW), or with circular (BTWC) boundary conditions. For the BTW/BTWC the exponent of $P(s)$ depends on the system size, and for the size used here ($M=2^{14}$) it takes the known value of $\approx 1.1$ [21]. Thus the slope is flatter than 1.29, which is expected for the infinite size BTW  [21, 22].

For subsampling, $N$ units were pre-chosen randomly. This subsampling scheme is well approximated by binomial subsampling with $p=N/M$ in the BM, because the BM runs on a network with full or annealed connections, and hence units are homogeneously connected. In the BTW/BTWC, subsampling violates the binomial sampling assumption, because of the models’ deterministic, local dynamics.

For all models, the avalanche distributions under full sampling transit from an initial power law to an exponential at a cutoff $s^{\mathrm{cutoff}}\approx M$ due to finite size effects (Fig. 2 A). For small $s$, the hairs appear in the BM, originating from subsampling power laws (Fig. 2 B, see Fig. 5 A for a flattened version). These hairs are almost absent in the BTW/BTWC, because the power-law slope is close to unity (5.1.3). The tails, even those of the BTWC, which have an unusual transition at the cutoff, collapse well. The BTWC is an exception in that it has unusual finite size effects, translating to the characteristic tails of $P(s)$. In fact, here the tails collapse better when applying fixed instead of binomial subsampling (Fixed subsampling refers to pre-choosing a fixed set of units to sample from; this may violate mean-field assumptions). This is because loosely speaking, binomial subsampling acts as a low pass filter on $P(s)$, smearing out the peaks, while fixed subsampling conserves the shape of the tails better here, owing to the compactness of the avalanches specifically in the 2D, locally connected BTWC. Overall, despite the models’ violation of mean-field assumptions, the analytically motivated p-scaling ansatz allows to infer $P(s)$ from subsampling, including the detailed shapes of the tail.

2.3 Distinguishing critical from subcritical systems

Distinguishing between critical and subcritical systems under subsampling is particularly important when testing the popular hypothesis that the brain shows signatures of “critical dynamics”. Criticality is a dynamical state that maximizes information processing capacity in models, and therefore is a favorable candidate for brain functioning [18, 23, 24, 25]. Typically, testing for criticality in experiments is done by assessing whether the “neural avalanche” distributions follow power laws [12]. Here, subsampling plays a major role, because at criticality avalanches can propagate over the entire network of thousands or millions of neurons, while millisecond precise sampling is currently constrained to about 100 neurons. Numerical studies of subsampling reported contradictory results [4, 8, 5, 10, 9]. Therefore, we revisit subsampling with our analytically derived scaling, and compare scaling for critical and subcritical states.

In contrast to critical systems, subcritical ones lack large avalanches, and the cutoff of the avalanche size distribution is independent of the system size (if $M$ is sufficiently large; Fig. 9). As a consequence, the distributions obtained under subsampling do not collapse under p-scaling (Fig. 2 C). In fact, there exists no scaling that can collapse all subsampled distributions (for any $p$) simultaneously, as outlined below, and thereby p-scaling can be used to distinguish critical from non-critical systems.

The violation of p-scaling in subcritical systems arises from the incompatible requirement for scaling at the same time the power-law part, the exponential tail, and the cutoff onset $s_{\mathrm{sub}}^{\mathrm{cutoff}}$ . On the one hand, the exponential tail becomes increasingly steeper with distance from criticality (larger $\lambda$), so that the relation $\lambda\ll p$ required for p-scaling (Eq. 6) does not hold anymore for small $p$ (Sec. 5.3). Thus, a collapse of the tails would require the scaling ansatz for exponentials (Eq. 4). On the other hand, slightly subcritical models still exhibit power-law behavior up to a cutoff $s^{\mathrm{cutoff}}:=c$ that is typically much smaller than the system size ($c\ll M$). To properly scale this part of the distribution, p-scaling is required. Likewise, the onset of the cutoff scales under subsampling with $p$: $s_{\mathrm{sub}}^{\mathrm{cutoff}}=c\cdot p$, requiring a scaling of the $s$-axis in the same manner as in the p-scaling. Thus, because the exponential decay requires different scaling than the power law and $s_{\mathrm{sub}}^{\mathrm{cutoff}}$, no scaling ansatz can collapse the entire distributions from “head to tail”.

2.4 Cluster definition, binning, and subsampling scaling

The main focus of this paper is to show how distributions of avalanches, node degrees or other “clusters” change under spatial subsampling, and how to infer the distribution of the fully sampled system from the subsampled one. To this end, it is essential that the clusters are extracted unambiguously, i.e. one cluster in the full system translates to exactly one cluster (potentially of size zero) under subsampling. This condition is easily realized for the degree of a node: One simply takes into account only those connections that are realized with other observed nodes. For avalanches, this condition can also be fulfilled easily if the system shows a separation of time scales (STS), i.e. the pauses between subsequent avalanches are much longer than the avalanches themselves (see Sec. 4.2.4). With a STS, temporal binning [12] can be used to unambiguously extract avalanches under subsampling. However, the chosen bin size must neither be too small nor too large: If too small, a single avalanche on the full system can be “cut” into multiple ones when entering, leaving, and re-entering the recording set. This leads to steeper $P_{\mathrm{sub}}(s)$ with smaller bin size (Fig. 3 A). In contrast, if the bin size is too large, subsequent avalanches can be “merged” together. For a range of intermediate bin sizes, however, $P_{\mathrm{sub}}(s)$ is invariant to changes in the bin size. In Fig. 3 A, the invariance holds for all bin sizes $32<\mathrm{bin\>size}<\infty$. The result does not depend on the topology of the network (compare Fig. 3 A for a network with sparse topology and Fig. 3 D for fully connected network). If a system, however, lacks a STS, then $P_{\mathrm{sub}}(s)$ is expected to change for any bin size. This may underlie the frequently observed changes in $P_{\mathrm{sub}}(s)$ in neural recordings [4, 7, 12, 26, 27, 28], as discussed in [8].

To demonstrated the impact of the bin size on p-scaling, we here used the branching model (BM), which has a full STS, i.e. the time between subsequent avalanches is mathematically infinite. When sampling $N=2^{7}$ out of the $M=2^{14}$ units, then $P_{\mathrm{sub}}(s)$ deviates from a power law for small bin sizes and only approaches a power law with the expected slope of $1.5$ for bin sizes larger than 8 steps (Fig. 3 A). The same holds for subsampling of any $N$: With sufficiently large bin sizes, $P_{\mathrm{sub}}(s)$ shows the expected approximate power law (Fig. 3 B). In contrast, for small bin sizes avalanches can be cut, and hence $P_{\mathrm{sub}}(s)$ deviates from a power law (Fig. 3 C). This effect was also observed in [8, 9], where the authors used small bin sizes and hence could not recover power laws in the critical BM under subsampling, despite a STS. Thus in summary, p-scaling only collapses those $P_{\mathrm{sub}}(s)$, where avalanches were extracted unambiguously, i.e. a sufficiently large bin size was used (compare Fig. 3 E and F).

The range of bin sizes for which $P_{\mathrm{sub}}(s)$ is invariant depends on the specific system. For the experiments we analyzed in the following section, we found such an invariance for bin sizes from 0.25 ms to 8 ms if $P_{\mathrm{sub}}(s)$ follows a power law, indicating indeed the presence of a STS (Fig. 4 D). Thus our choice of 1 ms bin size suggests an unambiguous extraction of avalanches, and in this range p-scaling works as predicted theoretically.

2.5 Subsampled neural recordings: Learning more by sampling less

We applied p-scaling to neural recordings of developing networks in vitro to investigate whether their avalanches indicated a critical state. To this end, we evaluated recordings from $N=58$ multi-units (see Methods, [29]). This is only a small fraction of the entire neural network, which comprised $M\approx 50.000$ neurons, thus the avalanche size distribution obtained from the whole analyzed data is already a subsampled distribution $P_{\mathrm{sub}}(s)$. To apply p-scaling, we generated a family of distributions by further subsampling, i.e. evaluating a subset $N^{\prime}<N$ of the recorded units. In critical systems, p-scaling is expected to collapse this family of distributions if avalanches are defined unambiguously, as outlined above (Sec. 2.4).

Interestingly, for early stages of neural development, p-scaling does not collapse $P_{\mathrm{sub}}(s)$, but for the more mature networks we found a clear collapse (Fig. 4; for all experiments see Fig. 13). Thus developing neural networks start off with collective dynamics that is not in a critical state, but with maturation approach criticality [30, 31]. Some of the mature networks show small bumps in $P_{\mathrm{sub}}(s)$ at very large avalanche sizes ($s\approx 5000\Leftrightarrow s/N\approx 60$). These very large avalanches comprise only a tiny fraction of all avalanches (about 2 in 10,000). At first glance, the bumps are reminiscent of supercritical systems. However, supercritical neural models typically show bumps at system or sampling size ($s=N$), not at those very large sizes. We discuss this in more detail in Sec. 5.5, and suggest that the bumps are more likely to originate from neurophysiological finite size effects.

For the full, mature network, our results predict that $P(s)$ would extend not only over three orders of magnitude as here, but over six, because $p\approx 10^{-3}$ (Sec. 5.5). Our analysis of neural recordings illustrates how further spatial subsampling allows to infer properties of the full system, even if only a tiny fraction of its collective dynamics has been observed, simply by sampling even less ($N^{\prime}<N$) of the full system.

2.6 Subsampling versus finite size scaling

In the real world we are often confronted with data affected by both subsampling and finite system size effects, i.e. observations originated from a small part of a large, but not infinite system. Thus we need to deal with a combination of both: subsampling effects as a result of incomplete data acquisition and finite-size effects inherited from the full system. To disentangle influences from system size and system dynamics, finite size scaling (FSS) has been introduced [32, 33]. It allows to infer the behavior of an infinite system from a set of finite systems. At a first glance, finite size and subsampling effects may appear to be very similar. However, if they were, then distributions obtained from sampling $N$ units from any system with $N\leq M$ would be identical, i.e. independent of $M$. This is not the case, as e.g. the distributions for fixed $N=2^{6}$ clearly depend on $M$ (Fig. 5 B). In fact, in both models the tails clearly inherit signatures of the full system size. Moreover, in the BM, subsampling a smaller fraction $p=N/M$ of a system increases the “hairs”, an effect specific to subsampling, not to finite size (see the increasing convexity of the flat section with decreasing $p$ in the BM, Fig. 5 B).

Importantly, as shown above, for critical systems one can always scale out the impact of subsampling, and thereby infer the distribution of the full system, including its size specific cutoff shape (Fig. 5 A). Hence, it is possible to combine FSS and subsampling scaling (detailed derivation are in Sec. 5.6): Consider a critical system, where FSS is given by: $M^{\beta}P(sM^{\nu};M)=g(s),$ here $g(s)$ is a universal scaling function. Then FSS can be combined with subsampling scaling to obtain a universal subsampling-finite-size scaling:

[TABLE]

Using Eq. 7 allows to infer the distribution for arbitrary subsampling ($N$) of any system size ($M$), Fig. 6.

3 Discussion

The present study analytically treats subsampling scaling for power laws (with cutoff), exponential distributions, and negative and positive binomial distributions. For all other distributions, utmost care has to be taken when aiming at inferences about the full system from its subsampling. One potential approach is to identify a scaling ansatz numerically, i.e. minimizing the distance between the different $P_{\mathrm{sub}}(s)$ numerically, in analogy to the approach for avalanche shape collapse [7, 34, 26, 35, 36]. We found that for our network simulations such a numerical approach identified the same scaling parameters as our analytic derivations (Sec. 5.7). However, given the typical noisiness of experimental observations, a purely numerical approach should be taken with a grain of salt, as long as it is not backed up by a circular form analytical solution.

Our analytical derivations assumed annealed sampling, which in simulations was well approximated by pre-choosing a random subset of neurons or nodes for sampling. Any sampling from randomly connected networks is expected to lead to the same approximation. However, in networks with e.g. local connectivity, numerical results depend strongly on the choice of sampled units [4]. For example, for windowed subsampling (i.e. sampling a local set of units) a number of studies reported strong deviations from the expected power laws in critical systems or scale free networks [4, 5, 6]. In contrast, random subsampling, as assumed here for our analytical derivations, only leads to minor deviations from power laws (hairs). Thus to diminish corruption of results by subsampling, future experimental studies on criticality should aim at implementing random instead of the traditional windowed sampling, e.g. by designing novel electrode arrays with pseudo-random placement of electrodes on the entire area of the network. In this case, we predict deviations from power laws to be minor, i.e. limited to the “hairs” and the cutoff.

We present here first steps towards a full understanding of subsampling. With our analytical, mean-field-like approach to subsampling we treat two classes of distributions and explore corresponding simulations. In future, extending the presented approach to a window-like sampling, more general forms of correlated sampling, and to further classes of distributions will certainly be of additional importance to achieve unbiased inferences from experiments and real-world observations.

4 Methods

4.1 Analytical derivations

The analytical derivations are detailed in the Supplementary Information.

4.2 Simulations

We simulated two models, the Bak-Tang-Wiesenfeld Model (BTW) with open and with circular (i.e. periodic) (BTWC) boundary conditions, and the branching model (BM) with full and with annealed sparse connectivity.

4.2.1 Bak-Tang-Wiesenfeld Model

The Bak-Tang-Wiesenfeld Model (BTW) [13], was realized on a 2D grid of $L\times L=M$ units, each unit connected to its four nearest neighbors. Units at the boundaries or edges of the grid have either 3 or 2 neighbors, respectively (open boundary condition). Alternatively, the boundaries are closed circularly, resulting in a torus (circular or periodic boundary condition, BTWC). Regarding the activity, a unit at the position $(x,y)$ carries a potential $z(x,y,t)$ at time $t$, ($z,t\in\mathds{N}_{0}$). If $z$ crosses the threshold of 4 at time $t$, its potential is redistributed or “topples” to its nearest neighbors:

[TABLE]

$z(x\pm 1,y\pm 1)$ refers to the 4 nearest neighbors of $z(x,y)$. The BTW/BTWC is in an absorbing (quiescent) state if $z(x,y)<4$, for all $(x,y)$. From this state, an “avalanche” is initiated by setting a random unit $z(x,y)$ above threshold: $z(x,y,t+1)=z(x,y,t)+4$. The activated unit topples as described above and thereby can make neighboring units cross threshold. These in turn topple, and this toppling cascade propagates as an avalanche over the grid until the model reaches an absorbing state. The size $s$ of an avalanche is the total number of topplings. Note that the BTW/BTWC are initialized arbitrarily, but then run for sufficient time to reach a stationary state. Especially in models with large $M$ this can take millions of time steps.

The BTW and the BTWC differ in the way how dissipation removes potential from the system. Whereas in the BTW potential dissipates via the open boundaries, in the BTWC an active unit is reset without activating its neighbors with a tiny probability, $p_{\mathrm{dis}}=10^{-5}$. For BTW an additional dissipation in a form of small $p_{\mathrm{dis}}$ can be added to make the model subcritical.

4.2.2 Branching model

The branching model (BM) corresponds to realizing a classical branching process on a network of units [8, 17, 18]. In the BM, an avalanche is initiated by activating one unit. This unit activates each of the $k$ units it is connected to with probability $p_{\mathrm{act}}$ at the next time step. These activated units, in turn, can activate units following the same principle. This cascade of activations forms an avalanche which ends when by chance no unit is activated by the previously active set of units. The control parameter of the BM is $\sigma=p_{\mathrm{act}}\cdot k$. For $\sigma=1$, the model is critical in the infinite size limit. We implemented the model with full connectivity ($k=M$) and with sparse, annealed connectivity ($k=4$). The BM can be mathematically rigorously associated with activity propagation in an integrate and fire network [37, 38].

For implementation of the BM with full connectivity ($k=M=2^{14}$), note that the default pseudo-random number generator (PRNG) of Matlab(R) (R2015b) can generate avalanche distributions that show strong noise-like deviations from the expected power-law distribution. These deviations cannot be overcome by increasing the number of avalanches, but by specifying a different PRNG. We used the “Multiplicative Lagged Fibonacci” PRNG for the results here, because it is fairly fast.

4.2.3 Subcritical models

To make the models subcritical, in the BM $\sigma$ was set to $\sigma=0.9$, and in the BTW/BTWC the dissipation probability $p_{\mathrm{dis}}$ was set to $p_{\mathrm{dis}}=0.1$, which effectively corresponds to $\sigma=0.9$, because 90% of the events are transmitted, while 10% are dissipated.

4.2.4 Avalanche extraction in the models

The size $s$ of an avalanche is defined as the total number of spikes from the seed spike until no more units are active. Under subsampling, this translates to the total number of spikes that occur on the pre-chosen set of sampled units (fixed subsampling). In principle, the avalanches could also have been extracted using the common binning approach [12], as all the models were simulated with a separation of time scales (STS), i.e. the time between subsequent avalanches is by definition much longer than the longest-lasting avalanche. Hence applying any bin size that is longer than the longest avalanche, but shorter than the pauses between avalanches would yield the same results for any subsampling.

4.3 Neural recordings

4.3.1 Data acquisition and analysis

The spike recordings were obtained by Wagenaar et al. [29] from an in vitro culture of $M\approx 50,000$ cortical neurons. Details on the preparation, maintenance and recording setting can be found in the original publication. In brief, cultures were prepared from embryonic E18 rat cortical tissue. Recording duration of each data set was at least 30 min. The recording system comprised an $8\times 8$ array of 59 titanium nitride electrodes with 30 $\mu$m diameter and 200 $\mu$m inter-electrode spacing, manufactured by Multichannel Systems (Reutlingen, Germany). As described in the original publication, spikes were detected online using a threshold based detector as upward or downward excursions beyond 4.5 times the estimated RMS noise [39]. Spike waveforms were stored, and used to remove duplicate detections of multiphasic spikes. Spike sorting was not employed, and thus spike data represent multi-unit activity.

For the spiking data, avalanches were extracted using the classical binning approach as detailed in [12, 8]. In brief, temporal binning is applied to the combined spiking activity of all channels. Empty bins by definition separate one avalanche from the next one. The avalanche size $s$ is defined as the total number of spikes in an avalanche. The bin size applied here was 1 ms, because this reflects the typical minimal time delay between a spike of a driving neuron and that evoked in a monosynaptically connected receiving neuron, and because 1 ms is in the middle of the range of bin sizes that did not change the avalanche distribution $P_{\mathrm{sub}}(s)$ (Sec. 2.4, Fig. 4 D).

Application of p-scaling by definition requires that one avalanche in the full system translates to one avalanche (potentially of size zero) under subsampling, i.e. an avalanche must not be “cut” into more than one, e.g. when leaving and re-entering the recording set. This can be achieved in experiments that have a separation of time scales by applying a sufficiently large bin size, because this allows for an unambiguous avalanche extraction [8] (Sec. 2.4). Indeed, the in vitro recordings we analyze here appear to show a separation of time scales: We found that varying the applied bin size around 1 ms hardly changed $P_{\mathrm{sub}}(s)$ (Fig. 4 D). In contrast, using too small bin sizes would have led to “cutting” avalanches, which impedes the observation of power laws, and consequently prevents the collapse (illustrated for the BM, Fig. 3).

4.3.2 Data availability and selection criteria

We evaluated 10 recordings for each day, because then the naïve probability of finding the expected behavior consistently in all of them by chance is at most $p=(1/2)^{10}<0.001$. The experimental data was made available online by the Potter group [29]. In detail, we downloaded from the “dense” condition the in vitro preparations 2-1, 2-3, 2-4, 2-5, 2-6, 6-1, 6-2, 6-3, 8-1, 8-2, 8-3, and for each preparation one recording per week (typically days 7, 14, 21, 28, 34/35, but for some experiments one or two days earlier), except for experiment 6-2 where we only got the first three weeks, and 6-3 where we got the last two weeks. We analyzed and included into the manuscript all recordings that we downloaded.

Acknowledgments and Funding

The authors thank Georg Martius and Michael Wibral for inspiring discussions and support. We thank Jens Wilting, Johannes Zierenberg, and João Pinheiro Neto for careful proof reading. AL received funding from the People Program (Marie Curie Actions) of the European Union’s Seventh Framework Program (FP7/2007-2013) under REA grant agreement no. [291734]. VP received financial support from the Max Planck Society. AL and VP received financial support from the German Ministry for Education and Research (BMBF) via the Bernstein Center for Computational Neuroscience (BCCN) Göttingen under Grant No. 01GQ1005B.

Author contributions

AL and VP contributed equally to the work. All results were developed jointly.

5 Supplementary Information for “Subsampling scaling: a theory about inference from partly observed systems”

In the following, we derive in detail the novel subsampling scaling. We first introduce the definition and basic results of subsampling in analogy to Stumpf et al. [3], who treated subsampling of graphs. We then extend the aforementioned study as follows:

• •

First, we focus on an analytical inference of the distribution of the full system from the subsampling, a topic that was not touched by Stumpf et al. To this end we derive (a) the exact subsampling scaling for negative binomials and exponentials, and (b) the approximate scaling for power-law distributions.

• •

Second, we explicitly show how to derive the system size from subsampling induced deviations from power laws (“hairs”).

• •

Third, we treat the relation between subsampling and finite size effects.

• •

Last, we apply our analytically derived subsampling scaling ansatz to infer the probability distribution of avalanche sizes in developing neural networks.

5.1 Mathematical subsampling

Let $X$ be a discrete, non-negative random variable with probability distribution $P(X=s)=P(s)$, with $s\in\mathbb{N}_{0}$, then $G(z)=\sum_{s=0}^{\infty}z^{s}P(s)$ is the corresponding probability generating function (PGF). For distributions such as power laws, where $s=0$ is not supported, $s$ is constrained to $s\in\mathbb{N}$, and the probability distribution needs to be normalized accordingly. $X$ represents the size of a set of “events” that comprise a “cluster”, e.g. the number of spikes in an avalanche, or the degree of a node. For subsampling, we assume that each of the events in the cluster is sampled independently with probability $p$, resulting in a random variable for the observed cluster size, $X_{\mathrm{sub}}$ [3]. Thus the probability $P_{\mathrm{sub}}(X_{\mathrm{sub}}=s)$ to observe a cluster of size $s$ is derived using a binomial distribution:

[TABLE]

The PGF $G_{\mathrm{sub}}(z;p)$ for $X_{\mathrm{sub}}$ with given $p$ is thus:

[TABLE]

Thus the PGF of $X$ and $X_{\mathrm{sub}}$ show a direct relation [3]:

[TABLE]

As a consequence, the expected values of $X$ and $X_{\mathrm{sub}}$ also are closely related: Using the expression for the expected value of $X$, $E(X)=G^{\prime}(1^{-})$

[TABLE]

These relations hold for any $P(s)$, however, only for specific $P(s)$, namely positive and negative binomials, the full and subsampled system’s $P(s)$ follow the same family of distributions [3], e.g. if $P(s)$ is a binomial distribution, then $P_{\mathrm{sub}}(s)$ also is a binomial, but with different parameters.

5.1.1 Subsampling of negative binomial and exponential distributions

Assuming that $X$ follows a negative binomial distribution $X\sim{\mathrm{NB}}(r,p_{{\mathrm{NB}}})$,

[TABLE]

then the expectation of $X$ is given by

[TABLE]

and the PGF is

[TABLE]

Using Eq. 8 then returns the PGF under subsampling,

[TABLE]

which corresponds to the negative binomial distribution with the same $r$, but different $p_{{\mathrm{NB}}}$, selected such that

[TABLE]

In the special case of $r=1$, the negative binomial is a geometric distribution with probability parameter $p_{{\mathrm{NB}}}$, and the discrete exponential distribution is a particular parametrization of the geometric distribution $1-p_{{\mathrm{NB}}}=e^{-\lambda}$.

[TABLE]

Using equation 11, the relation between $\lambda$ and $\lambda_{\mathrm{sub}}$ is:

[TABLE]

Solving this equation with respect to $\lambda_{\mathrm{sub}}$ we obtain:

[TABLE]

5.1.2 Subsampling of power-law distributions

To derive an approximate scaling for power-law distributions under subsampling, we expand on the work by Stumpf et al. [3]. Consider mathematical subsampling as defined in the main text, and a power-law distribution $P(s)=C_{\gamma}s^{-\gamma}$ with exponent $\gamma>1$, and normalization $C_{\mathrm{\gamma}}=1/\zeta(\gamma)$, where $\zeta(\gamma)$ is the Riemann zeta function. Then $P_{\mathrm{sub}}(s;\gamma,p)$ is a binomial subsampling with sampling probability $p$:

[TABLE]

Building on the work by Stumpf et al., we assume that for the tail (i.e. for $s\rightarrow\infty$) the subsampled distribution is approaching an appropriately scaled power-law with slope $\gamma_{\mathrm{sub}}=\gamma$, i.e. we assume

[TABLE]

$c_{\gamma}(p)$ is the subsampling-dependent normalization constant. To derive how $c_{\gamma}(p)$ depends on $p$, we need to assume that

[TABLE]

This is a strong assumption, because typically an exchange of differentiation and limit is only possible in case of uniform convergence of the derivatives [40], which is not the case here. However, all the functions we consider are monotonous in all parameters and numerical results support the assumption above.

In the following we assume that $s$ is large enough so that Eq. 14 can be taken as an identity. Then

[TABLE]

The first term can be approximated as:

[TABLE]

The second term, after introducing $k=n-1$, reduces to:

[TABLE]

The “$\approx$” is inherited from Eq. 14, which is only exact for $s\rightarrow\infty$. Combining the two terms, we obtain for Eq. 15:

[TABLE]

From this, $\frac{\partial}{\partial p}c_{\gamma}(p)$ can be expressed as:

[TABLE]

For solving the limit, we use the known identity

[TABLE]

which can be restated by replacing $x$ by $1/s$ and $\mu$ by $-\gamma$:

[TABLE]

Thus

[TABLE]

This differential equation is solved by:

[TABLE]

For $p=1$ we know that $C^{*}=C_{\gamma}$, because sampling all units does not change the distribution. The final expression for $c_{\gamma}(p)$ is thus

[TABLE]

With this we can derive scaling parameters $a,b$ that collapse the distribution’s tails, i.e. $p^{a}P_{\mathrm{sub}}(p^{b}s)=P(s)$ for large $s$. Using Eq. 20:

[TABLE]

Thus for any $a$ and $b$, such that $a-b\gamma=1-\gamma$, the scaling ansatz leads to a collapse. One of the members of this scaling family is $b=0$, $a=1-\gamma$, which scales the $y$-axis only. As shown in Fig. 7, this scaling collapses the tails of distributions perfectly. For small $s$, however, there are systematic deviations under subsampling, which increase with smaller $p$. We call them “hairs”, because they grow on the head of the distribution, as opposed to the tails.

A different member of the scaling family is $a=b=1$. This scaling is especially attractive, because it does not require information about the exponent $\gamma$ of the power law (see Fig. 1). As this scaling is linear in $p$, we call it p-scaling.

5.1.3 Power-law exponent close to unity

Here we show why the “hairs” become smaller, i.e. converge to zero, in the limit of the power-law exponent $\gamma\to 1$. It is in agreement with results of Stumpf et al. [3], stating that “hairs” are growing with increase of the exponent. Mathematically, the exponent of the power-law distribution cannot be exactly equal to one or smaller, because in this case the distributions cannot be normalized. Thus without loss of generality we consider truncated power laws: $P(s)=C\cdot s^{-1}$ for $s\leq s_{\mathrm{max}}$ and $P(s)=0$ for $s>s_{\mathrm{max}}$. The normalizing constant $C$ depends on $s_{\mathrm{max}}$. In this case the subsampled distribution $P_{\mathrm{sub}}(s)$, with sampling probability $p$ can be written explicitly

[TABLE]

We are interested in the behavior of the “hairs” and thus consider small $s$. In this case, we can approximate $P_{\mathrm{sub}}(s;p)$ by the infinite sum, and make use of the geometric series

[TABLE]

to obtain, with a variable exchange $m=l-s$,

[TABLE]

Thus we showed that in the limit $\gamma\to 1$ subsampling of the power law converges to the original power law.

5.2 Inferring the system size from the subsampled distribution

The deviations from power laws (i.e. the hairs), which emerge under subsampling, allow to infer the system size $M$ from the subsampled distribution $P_{\mathrm{sub}}(s)$ alone, given that $P(s)$ follows a power law. This is because the hairs are a function of the sampling probability $p=N/M$. The hairs are most pronounced for $P_{\mathrm{sub}}(s=1)$ (except for $P_{\mathrm{sub}}(s=0)$, which may remain unknown in experiments). Therefore, the inference of system size in experiments is most accurate if it is based on $P_{\mathrm{sub}}(s=1)$. We explore this in the following derivations. Derivations based on other (small) $s$ can be performed analogously.

Quantitatively, using the explicit relation for subsampling of power laws (Eq. 13) with $l=n+1$ results in:

[TABLE]

where $\zeta$ is the Riemann zeta function, and $Li_{\gamma}(z)=\sum_{k=1}^{\infty}z^{k}/k^{\gamma}$ is the polylogarithm function. This relation is exact if $P(s)$ is a true power law. For application to the real data obtained from subsampled observation the following algorithm allows to infer $p$:

Check whether the experimentally obtained empirical distribution $P_{\mathrm{emp}}(s)$ is likely to originate from a system that under full sampling shows a power-law distribution. If not, the method cannot be applied. 2. 2.

Estimate the power-law slope $\gamma$ of the power-law tail of the distribution to obtain $\hat{\gamma}$. 3. 3.

Solve the following equation for $p$:

[TABLE]

This will return the sampling probability $p$. From this, the system size can be inferred if $N$ is known. This approach is also applicable approximately if the full system does not display a pure power law, but a power law with cutoff at large $s$. Then the power-law slope $\gamma$ has to be inferred on an appropriate interval between the hairs and the cutoff.

We applied this method numerically to the data generated by the critical branching model of size $M=1024$, subsampled to $N=2^{0},2^{1},\ldots 2^{9}$ units based on $10^{7}$ avalanches in the full system. Indeed, the full system size could be inferred by $\hat{M}=\hat{p}N$ with high precision (Fig 8): The maximal deviations were smaller than 6%.

5.3 Subcritical systems

As outlined in the main text, avalanche distributions collapse under p-scaling for critical systems, but not for subcritical systems. The main reason is that for subcritical systems the exponential tail is too steep, i.e. the requirement $\lambda\ll p$ is violated. We in the following derive an approximate relation between $\lambda$ and the distance to criticality ($\varepsilon=\sigma_{\mathrm{crit}}-\sigma=1-\sigma$). We show that for more subcritical systems, $\lambda$ becomes increasingly larger (see also Fig. 9 (left)). To approximately derive the relation between $\lambda$ and $\varepsilon$ (or $\sigma$), we used the branching process [17], because it allows to easily control the distance to criticality by changing the branching ratio $\sigma$, and because it is independent of finite size effects. This is a reasonable assumption, because in subcritical systems $P(s)$ is not affected by changing the system size for any $M>M_{0}$. Only for very small systems sizes there are finite size effects (Fig. 9 (right)).

To derive heuristically the slope of the exponential tail $\lambda$ as a function of the control parameter $\sigma$, consider a branching process with branching ratio $\sigma<1$, and assume that an avalanche starts with 1 neuron firing. Then on expectation in the second time step there are $\sigma$ neurons firing, in the third time step $\sigma^{2}$, and so forth. Thus we obtain an expression for the average avalanche size $\langle s\rangle$:

[TABLE]

In the subcritical regime, the distribution of the avalanche sizes is dominated by the exponential cutoff. We consider that they are well approximated by the power law with slope $\gamma$ and an exponential cutoff parametrized by $\lambda$

[TABLE]

The mean value of this distribution is given by:

[TABLE]

Where $Li_{\gamma}(z)=\sum_{k=1}^{\infty}z^{k}/k^{\gamma}$ is again the polylogarithm function. The relation between $\lambda$ and $\sigma$ is thus

[TABLE]

and hence $\lambda$ approaches zero when approaching the critical point ($\sigma\rightarrow 1$). As $\lambda$ decays slowly as a function of $\sigma$, except in the very close vicinity of the critical point, the requirement for p-scaling, $\lambda\ll p$, is only satisfied in the close vicinity of the critical point. Else p-scaling does not apply.

To compare our analytical with numerical results, we used the same data as in Fig. 9. We first estimated $\gamma\approx 1.3$ from the distributions, and with this solved equation 22 numerically. The analytical results closely fitted the slopes $\lambda$ of the exponentials from the simulations (Fig 10).

5.4 Subsampling scaling for the EHE-model and the sparsely connected branching model

In this section, we investigate whether subsampling scaling also applies to other models than the ones treated in the main manuscript. In particular, we treat here first the Eurich, Herrmann & Ernst (EHE) model [20], a classical extension of the BTW model to neural networks, and then a realization of the BM with sparse connectivity ($k=4$, see Methods). The details of the EHE model can be found in [20, 37]. This model produces power-law distributions of the avalanche sizes with slope $\approx 1.5$ that indicates that it belongs to the same university class as the branching model. However, activity transmission is not stochastic as in BM, but deterministic as in BTW. Another peculiarity of the model lays in its dissipative nature: for the finite system sizes $M$ each spike leads to dissipation of $\Delta\approx 1/\sqrt{M}$ [20]. Thus only in the limit $M\to\infty$ the model is both truly critical and conservative. We simulated the EHE model with both, the classical fully connected graph topology and also with random connectivity probability $p_{\mathrm{conn}}$. For both models, the dynamics is as follows: Each neuron $i$ is a non-leaky integrator, and its membrane potential is denoted by $h_{i}\in[0,1)$. When $h_{i}$ crosses the threshold $\theta=1$ the neuron fires and is reset $h_{i}\mapsto h_{i}-1$ and all its postsynaptic connections receive an input of strength $\alpha/M$, where $\alpha$ is the control parameter in the model, the strength of interaction. For the fully connected network, it is known that $\alpha\approx 1-M^{-0.5}$ leads to an approximate power-law distribution of the avalanche sizes. For not-fully connected networks the connection probability $p_{\mathrm{conn}}$ needs to be included, and thus the condition to achieve approximate power-law distributions generalizes to $\alpha\cdot p_{\mathrm{conn}}\approx 1-M^{-0.5}$.

As the avalanche size distribution in the EHE model can be directly mapped to the branching model [37], subsampling scaling is expected to behave the same as in the BM, producing “hairs” but resulting in a good collapse. We tested this for the model of $M=1024$ neurons with $p_{\mathrm{conn}}=0.1$ and obtained, as expected, a collapse under subsampling scaling (Fig. 11).

The distributions of the fully and the sparsely connected BM are very similar (Fig. 12). The only difference is a slightly more pronounced lack of small avalanches in the fully sampled sparse BM (Fig. 12 C), which translates to somewhat less pronounced “hairs”, in particular under “mild” subsampling ($N\geq 2^{10}$).

5.5 Detailed discussion of the experimental results

Figure 13 displays $P_{\mathrm{sub}}(s)$ for all recordings of developing neural cultures we evaluated (see Methods). As discussed in the main text (Fig. 4), with maturation the $P_{\mathrm{sub}}(s)$ approached power-law scaling, which for the fully sampled culture is expected to extend over almost six orders of magnitude. In addition to the power laws, about half of the mature cultures also showed a bump in $P_{\mathrm{sub}}(s)$ at very large sizes ($s\approx 5000$). These very large avalanches comprise only a tiny fraction of all avalanches ($\approx 0.02\%$). Such bumps are a priori not expected for critical systems. The collapse of the bumps itself is a manifestation of the activity spread during the large avalanches that hit the sampled set proportionally to the number of sampled units. In the following we discuss first whether the distributions with the bumps are expected to collapse under p-scaling, and then the potential origin of the bumps.

Regarding the questions whether the distributions observed here are expected to collapse, the answer is straight forward: The avalanches in the tail make only a tiny fraction of all observed avalanches (about 2 in 10,000), while the other 99.98% avalanches follow a power law for about 3 orders of magnitude. (It is the log-log scale together with the logarithmic binning that might make the bumps appear more prominent than they are.) With only 0.02% of avalanches not following a power law, a decent collapse is to be expected.

Regarding the origin of the bumps, there are a number of potential explanations, which we outline in the following. At first glance, the bumps are reminiscent of supercritical systems (hypothesis 3), however, they do not occur at $N$, the number of sampled units, as expected for supercritical systems. More likely, they may represent finite size effects (hypothesis 1), or alternatively transient switches to a bursty state (hypothesis 2). All three hypotheses are detailed in the following:

Biological finite size effects in a critical system.

Assume the neural cultures were precisely at a critical point. Thus the distribution of the avalanche sizes would be a perfect power law without cutoff. However, in biological systems the avalanche size cannot go to infinity, because biological mechanisms (e.g. depletion of synaptic resources, shortage of Ca2+ or homeostatic mechanisms) limit their maximal size. All these avalanches that are larger than some $s_{\mathrm{trans}}$ (in our data $s_{\mathrm{trans}}\approx 3000$) are thus expected to be distributed around a characteristic, biologically determined size, which here is about $s\approx 5000$. The probability $p_{\mathrm{bump}}$ to observe an avalanche larger than some maximal size $s_{\mathrm{trans}}$ is given by the Hurwitz zeta function. Indeed, in agreement with this hypothesis, the number of the avalanches observed in the “bump” agrees with the probability $p_{\mathrm{bump}}$ for perfectly critical system. Thus the data support our hypothesis that the bump represents the collection of all avalanches that would, in an ideal system, be larger than 3000. (Note that all avalanche sizes $s$ given here are the sizes observed under subsampling). Thus biological finite size effects are a probable origin for the bumps. 2. 2.

Criticality alternates with a state that gives rise to large avalanches

The in vitro neural networks we analyzed could in principle alternate between different states. While in one state, which comprises about 99.98% of the avalanches, the system is critical, in the other state it displays unusually large avalanches that run multiple times over the entire system and give rise to population bursts, i.e. they manifest as the observed bumps. The precise fraction of “burst avalanche” can depend on the properties of each individual culture (some showing none at all), and it could be pure coincidence that the fraction of burst avalanches is in agreement with the fraction expected for the avalanche tail (see hypothesis 1). 3. 3.

A novel form of slight supercriticality in a finite system.

While it is straight forward to identify “subcriticality” (no avalanches covering the full system size, no power-law behavior of distributions, but a prominent exponential tail), it is much trickier to identify “supercriticality” in neural systems by pure observation, potentially because supercriticality in the thermodynamic limit implies a non-zero fraction of infinite avalanches, but in finite systems it depends on the type of system how these infinite avalanches manifest. For supercritical systems in neuroscience, the bump in $P_{\mathrm{sub}}(s)$ occurs typically at $N$, i.e. the system size or the number of sampled neurons [12, 41]. However, here in all experiments where the bump is observed, it is around 80 times $N$ (i.e. $s\approx 5000$ from sampling up to 60 electrodes). Thus here the bumps do not indicate supercritical behavior resembling that of previous studies. However, it could indicate a novel form of supercriticality on a finite system.

How to distinguish between these potential causes of the bump appearance remains an open question for further experimental investigations (e.g. changing the network size; making the network on purpose supercritical). In the experiments evaluated here, the presence of the data collapse in the more mature networks predicts a power-law distribution for $P(s)$ of the full neural system that spans approximately 6 orders of magnitude. However, whether such power-laws scaling is sufficient to infer criticality, is still under debate.

5.6 Combining subsampling scaling and finite-size scaling

As demonstrated in section 2.6, there is a fundamental difference between subsampling scaling that deals with partial observations of a system, and finite-size scaling (FSS) that extrapolates from models of finite size to infinite size systems. Here we show how to combine both scaling ansätze to obtain a universal scaling.

The finite-size scaling ansatz for a critical system is formulated as:

[TABLE]

where $g(s)$ is a scaling function. The formulation for the subsampling scaling in a system with $M$ units is:

[TABLE]

This can be re-written as:

[TABLE]

Our goal is to combine the finite size scaling and the subsampling scaling relations (Eqs. 23 and 24) to factorize out the dependence of $g_{\mathrm{sub}}$ on $M$, and hence be able to collapse subsampled distributions from different system sizes $M$. To find the appropriate scaling, we need to identify the exponents $\delta$ and $\kappa$ such that

[TABLE]

To this end, recall that $P_{\mathrm{sub}}(s,M;M)=P(s,M)$. In the following we first use subsampling scaling to express $P_{\mathrm{sub}}(s,M;M)$ in terms of $P(s,M)$ using Eq. 24:

[TABLE]

Thus the solutions for the exponents is given by $\delta=\beta-1$ and $\kappa=\nu-1$, and hence the general subsampling-finite-size scaling is given by:

[TABLE]

We tested this relation numerically for the case of the branching model (BM). For this model, FSS is given by $\beta=1.5$, $\nu=1$ and thus the subsampling-finite-size scaling is given by:

[TABLE]

Indeed, with this scaling we obtained as expected a good collapse for combining different sampling sizes $N$ and system sizes $M$ (Fig. 6).

5.7 Numerical estimation of optimal scaling

Throughout the manuscript we used an analytical approach to determine the optimal scaling for the subsampled distributions. In this section we confirm that our analytical results coincide with a direct numerical estimation of the scaling constants $a,b$. To achieve an optimal scaling collapse, we numerically estimated the parameters $a$ and $b$ that minimize the distance $d(a,b)$ between the rescaled distributions $P_{\mathrm{sub}}(s;a,b,N)=N^{a}P_{\mathrm{sub}}(N^{b}s)$ under subsampling, and the rescaled distribution under full sampling, $P(s;a,b,N=M)$. In more detail, we first estimated for each $N\in[8,16,32,\dots,M]$ the distance $d(a,b;N)$ as follows:

[TABLE]

The mean $\left\langle\cdot\right\rangle_{s}$ was taken over all $s\in[1,10\cdot N]$, and $s^{\prime}$ are the support points in $P(s^{\prime};a,b,M)$ corresponding to those in $P_{\mathrm{sub}}(s;a,b,N)$, i.e. $s^{\prime}/M^{b}=s/N^{b}$ is fulfilled. As $s^{\prime}$ may take non-integer values, the values $P_{\mathrm{sub}}(s^{\prime};a,b,N)$ are obtained, if necessary, by linear interpolation between the nearest integers. Then the weighted average over all $d(a,b,N)$ is taken to obtain $d(a,b)=\left\langle d(a,b,N)\right\rangle_{N}$. The parameter combination $(a^{*},b^{*})$ that minimizes $d(a,b)$ provides numerically the optimal collapse. We scanned $a$ and $b$ in steps of $0.01$ and found for both, the BM and the BTWC an optimal collapse at $a^{*}=b^{*}=1$, as predicted analytically (Fig. 14). For the BTW, the optimal collapse was at $a^{*}=1.01,b^{*}=1.02$, but the value of $d(a,b)$ in the point analytically obtained $a=b=1$ deviated only by $3\%$ from the absolute minimum (Fig. 14). Thus overall, our numerical results match very well the theoretical prediction.