## Abstract

Current efforts to achieve neuromorphic computation are focused on highly organized architectures, such as integrated circuits and regular arrays of memristors, which lack the complex interconnectivity of the brain and so are unable to exhibit brain-like dynamics. New architectures are required, both to emulate the complexity of the brain and to achieve critical dynamics and consequent maximal computational performance. We show here that electrical signals from self-organized networks of nanoparticles exhibit brain-like spatiotemporal correlations and criticality when fabricated at a percolating phase transition. Specifically, the sizes and durations of avalanches of switching events are power law distributed, and the power law exponents satisfy rigorous criteria for criticality. These signals are therefore qualitatively and quantitatively similar to those measured in the cortex. Our self-organized networks provide a low-cost platform for computational approaches that rely on spatiotemporal correlations, such as reservoir computing, and are an important step toward creating neuromorphic device architectures.

## INTRODUCTION

There is currently huge interest in building neuromorphic computational systems that operate on brain-like principles to perform certain tasks (e.g., classification and pattern recognition) much more efficiently than on traditional computers. Neuromorphic processing has been implemented in a number of hardware platforms, using components ranging from traditional silicon-based technology (*1*, *2*), crossbar memristor arrays (*3*–*6*), and novel nanoscale device elements that use atomic-scale dynamics to mimic the behavior of synapses and neurons (*7*–*9*). All of these approaches rely on highly ordered architectures, but neuromorphic computational paradigms such as reservoir computing (RC) (*10*–*12*) require complex networks with strong inherent spatiotemporal correlations (*13*, *14*). Maximal computational performance (*15*–*17*) in these systems requires network architectures that are not only complex but also critical (*18*).

At a critical point (*17*–*20*), the underlying network is scale free (*21*, *22*), and the dynamics are characterized by scale-invariant avalanches, as observed in many systems ranging from earthquakes to biological extinctions and magnetization dynamics (*23*). Put simply, a critical system is poised such that any event is likely to trigger subsequent events, leading to a cascade (*15*). Avalanches of events are an important feature of both in vivo and in vitro recordings of neuronal signals (*24*–*27*), and substantial evidence has accumulated that the brain itself operates at a self-organized critical point (*20*), which optimizes information processing, memory, and information transfer (*15*–*19*). Therefore, in designing brain-inspired computational systems, it is natural to look for systems that might exhibit similar critical behavior and, especially, systems that would allow electronic signal processing.

### The percolation transition

Percolating systems (*28*) are obvious candidates for investigation since the concept of percolation (Fig. 1) is central to discussions of criticality and avalanches (*17*). The two-dimensional percolating system of interest here comprises conducting nanoparticles deposited on an insulating substrate (*29*) (Fig. 1A). During deposition (see Materials and Methods), particles come into contact and form groups that are separated from each other by tunnel gaps. The groups increase in size and complexity as more particles are deposited until the onset of conduction across the system (Fig. 1, B and D) when the fraction of the surface that is covered with conducting elements is *p*_{c} ~68% (*28*–*31*). This onset marks the percolation threshold, which is a critical point between insulating and conducting phases (*28*) at which avalanches might be expected to propagate (*17*) (see next paragraph). We hypothesized that when poised at this phase transition, the percolating assembly of nanoparticles might be promising for neuromorphic computing since (i) the correlation length diverges (*28*) and so long-range correlations are expected; (ii) it has previously been established that synapse-like atomic-scale switching processes occur in the tunnel gaps (*29*); (iii) devices are stable over periods of months (*32*) as required for real-world applications; and (iv) the complex networks are self-organized by simple, cost-effective processes (*32*).

Near the percolation threshold, tunneling across gaps between particles contributes to the conductivity of the network (see Fig. 1 and fig. S1), and it is necessary to go beyond standard approaches to continuum percolation that consider only ohmic connections between well-connected particles (*28*). Percolating-tunneling systems have, however, received relatively little theoretical attention [see (*31*, *33*) and references therein), and quantities such as the relevant percolation critical exponents are yet to be calculated. Experimentally, it is observed that external electrical stimuli cause conductive filaments to be continually formed and broken within the tunnel gaps (Fig. 1A, insets) (*29*), leading to a complex and dynamical network comprising conducting pathways, switches, and tunnel junctions (further illustrated in fig. S2).

### Behavior of critical networks

The importance of criticality for neuromorphic computing was discussed in previous work on silver nanowire networks (*13*), but criticality in self-organized nanodevices has never been investigated systematically. Here, we investigate criticality in percolating nanoparticle networks using analysis techniques similar to those previously used on recordings of neuronal signals (*24*–*27*).

Figure 1 and fig. S2 show how avalanches are generated in the dynamical percolating network. When an atomic switch is activated in one tunnel gap, consequent changes in the electric field (and current) in other tunnel gaps cause further events at other locations in the network (*34*). That is, each successive switching event leads to others, creating an avalanche. As illustrated in Fig. 1D and fig. S2, the surface coverage is important in controlling the network connectivity: If too many (or too few) events are triggered, then the avalanche will accelerate uncontrollably (or die quickly)—critical dynamics are characterized by power law distributions in which there is no single characteristic length or time scale (*17*) and in which there is a small but important chance that very large avalanches are observed [see Fig. 1, fig. S2, and (*15*)]. While the propagation of avalanches is commonly described as a percolation phenomenon (*17*), avalanches have not previously been reported in percolating networks.

## RESULTS

### Avalanches of electrical signals and correlations

The simple contact geometry (Fig. 1A) allows us to sensitively observe electrical signals from the entire network, and Fig. 2A shows complex patterns of switching events [observed as discrete changes in the device conductance (*G*)], which are qualitatively similar to those observed due to neuronal avalanches in cortical tissue (*24*–*27*). Figure 2A highlights the self-similar nature of the avalanches, with characteristic bursts of activity on multiple time scales. Quantitative analysis of the avalanches is discussed below. Figure 2 (B and E) shows that the measured changes in conductance (Δ*G*) range over about three orders of magnitude, with different signal sizes reflecting the locations of the switches on different pathways through the highly branched network. The distribution of measured Δ*G* values is heavy tailed (Fig. 2, B and E), consistent with spatial self-similarity in percolating structures (*35*).

Histograms of interevent intervals (IEIs) (Fig. 2, C and F) show clear power law behavior, implying scale-free dynamics. The corresponding autocorrelation functions (ACFs) are shown in Fig. 2 (D and G). Again, power law behavior is observed. The data shown in Fig. 2 (D and G) are obtained with sampling rates that are different by a factor of a thousand, and yet the power law exponents (β) obtained are almost identical (β ~0.2). Together, these data cover more than five orders of magnitude in time and provide strong evidence for temporal self-similarity (*36*). However, so-called renewal processes, where successive IEIs are uncorrelated, can also exhibit power law ACFs (*37*), and so we show in Fig. 2 (D and G) the ACFs obtained when the sequence of IEIs is randomly shuffled to destroy any correlations. The experimental data have a much smaller slope than the shuffled data, indicating long-range temporal correlations (*36*). We note that in the neuroscience literature, it is often difficult to directly measure ACFs, and so detrended fluctuation analysis is commonly used to quantify correlations [see (*36*, *38*) and citations therein]. For comparison, we use the relationship β = 2 − 2*H* (*38*) to estimate the Hurst exponent (*H*). We find *H* ~0.9, consistent with the strong correlations described below. This is a first piece of evidence for highly correlated activity in our networks, but the more detailed analysis below is required to unambiguously demonstrate avalanches and criticality.

### Demonstration of criticality

In critical systems, the sizes (*S*) and durations (*T*) of avalanches have probability density functions (PDFs) that follow well-defined power law distributions$$encoding>P(S)\sim {S}^{\u2013\tau}$$(1)$$encoding>P(T)\sim {T}^{\u2013\alpha}$$(2)and the mean sizes of the avalanches are related to their durations by$$encoding>\u3008S\u3009(T)\sim {T}^{1/\mathrm{\sigma}\mathrm{\nu}\mathit{z}}$$(3)where 1/σν*z* is a key characteristic exponent (*23*). As explained in (*39*), in critical systems, the estimate of 1/σν*z* obtained from Eq. 3 should be consistent with that estimated from the “crackling relationship”$$encoding>1/\mathrm{\sigma}\mathrm{\nu}\mathit{z}=(\alpha -1)/(\tau -1)$$(4)

In addition, the shapes of the avalanches should also exhibit a characteristic “collapse” onto a single curve (*23*, *25*, *40*), yielding an additional independent estimate of 1/σν*z*. To demonstrate criticality (i.e., to distinguish from other processes, such as random walks, that can also lead to power law distributed avalanches), these three independent estimates of 1/σν*z* must be in agreement (*23*, *27*, *39*).

Avalanches in our percolating networks have sizes and durations (see Materials and Methods) that are distributed in accordance with power laws (Eqs. 1 and 2), as shown in Fig. 3 (A and F) and Fig. 3 (B and G). After shuffling the sequence of IEIs in the original data, correlations are destroyed, and both distributions become exponential (gray curves) (*41*), confirming that the original data are highly correlated. Figure 3 (C and H) shows that the mean avalanche size (〈*S*〉) depends on avalanche duration (*T*) as expected from Eq. 3, with 1/σν*z* ~1.6. Physically, the slope of the 〈*S*〉(*T*) plots is related to the fractal dimension of the avalanches (*23*).

Figure 3 (D, E, I, and J) shows that the avalanches exhibit shape collapse. As is the case for neuronal recordings (*25*, *40*), this requires an enormous number of avalanches and necessarily yields an estimate of 1/σν*z* that has a relatively large uncertainty. The avalanche profiles are consistent with the parabolic scaling functions reported for neuronal avalanches (*25*). Detailed investigation of avalanche shapes may eventually yield further understanding of the underlying processes (*42*, *43*), but such an analysis would require a much larger amount of data so as to improve the signal-to-noise ratio in these plots. Nevertheless, as shown in Fig. 4, the values of 1/σν*z* obtained from shape collapse are consistent with those obtained from the crackling relationship (Eq. 4) and the mean avalanche sizes (Eq. 3). We emphasize that the important feature of Fig. 4 is that in almost every case, the blue, cyan, and green symbols agree within the uncertainties. This is true for a number of devices, for a range of stimuli, and for repeated independent measurements on the same devices. The fact that all three estimates of 1/σν*z* are in agreement is very strong evidence for criticality in our percolating devices (*39*). The agreement of all estimates of 1/σν*z* even while α and τ vary (Fig. 4) is evidence for self-tuned criticality (*18*).

## DISCUSSION

The avalanches of switching events observed here exhibit statistics that are qualitatively and quantitatively similar to those of the avalanches observed in cortical tissue (*25*, *41*) and in very recent work on whole organisms (*26*). Within the brain, criticality is hypothesized to play a role in cognition, adaptive behavior, and disease diagnosis and underpins several therapeutic and diagnostic opportunities (*19*). Criticality in our system originates from long-range spatiotemporal correlations in the network of atomic switches, as required (*18*) for neuromorphic computing strategies such as RC (*10*).

In RC, spatiotemporal correlations are required to allow projection of input signals into higher-dimensional spaces; corresponding output signals can then be combined using linear regression to perform computational tasks such as classification, pattern recognition, and time series prediction (*10*–*12*). The need to build critical systems to provide a platform for RC has previously been highlighted in (*17*, *18*). Our percolating networks meet that need.

Previous demonstrations of RC using nanoscale devices have relied on small numbers of device elements (*11*, *12*) because integration into appropriate physical networks is challenging and/or expensive. Furthermore, it is not clear that complex/critical networks could be achieved using lithographic processing. In contrast, the straightforward deposition methods used here immediately provide the required complex network of switches, are inexpensive, and allow electronic readout of correlated signals from the network as required for integration with other electronic components. We believe our processes are scalable because they are compatible with standard complementary metal-oxide semiconductor (CMOS) processing.

Our percolating networks are naturally suited to RC since deposition onto multielectrode arrays will straightforwardly provide the required multiple inputs and outputs. Ultimately, we envisage architectures in which the percolating network is deposited onto chips with predefined CMOS circuitry that processes input and output signals. The question of how best to code input/output signals for particular computational tasks has so far been given little consideration in the RC literature and will require further investigation since combinations of temporal and spatial input/output coding are possible (*12*).

We believe that percolating devices could provide many further opportunities for brain-like computing; for example, three-dimensional percolating systems could provide additional network complexity. While our attempts to tune the operating point of the system by varying the surface coverage were unsuccessful (see Materials and Methods), it may yet be possible to tune the system through the critical point using other external control parameters. In addition, strategies can be envisaged that would incorporate a wide range of alternative switching elements within the networks, such as phase-change materials, oxide-based memristors, and electrochemical switches (*6*–*9*). Last, we highlight the need for theoretical work/simulations to support device development. Percolating-tunneling networks have so far been relatively little studied (*31*, *33*), and important quantities such as the degree distribution for network connections (*21*), percolation critical exponents (*28*), and precise universality class (*23*) are yet to be elucidated.

## MATERIALS AND METHODS

### Device fabrication

Our percolating devices were fabricated by simple nanoparticle deposition processes (*32*, *44*). Seven-nanometer Sn nanoparticles were deposited between gold electrodes (spacing, 100 μm) on a silicon nitride surface. Deposition was terminated at the onset of conduction, which corresponds to the percolation threshold (*28*, *44*). The deposition takes place in a controlled environment with a well-defined partial pressure of air and humidity, as described in (*32*). This process leads to controlled coalescence and fabrication of robust structures, which function for many months but which yet allow atomic-scale switching processes to take place unhindered. The individual atomic-scale filaments formed within tunneling gaps exhibit quantized conduction (*29*), but, as shown in fig. S1 (E and F), the complex percolating-tunneling network comprises many parallel and series paths, and so the device conductance is not quantized.

As described in the discussion of Fig. 4, data from all devices presented here are consistent with criticality. To explore noncritical networks, several devices were fabricated with surface coverages that were either below or above the percolation threshold (see Fig. 1, B and D). In all cases, however, it was not possible to observe avalanches of events. When stimulated, low-coverage devices exhibit a number of irreversible drops in conductance and rapidly become open circuit: Atomic switches can be opened, but the prevalence of large gaps in the network makes closing switches difficult. High-coverage devices exhibit no switching events for small stimulus voltages and exhibit irreversible drops in conductance for larger stimulus voltages due to melting of particles when large currents flow.

### Electrical stimulus and measurement

Electrical stimuli were applied to the contact on one side of the percolating device, while the opposite edge of the system was held at ground potential. While a variety of types of stimuli (voltage pulses and ramps) can be applied to the system, constant DC voltages were used in this work because they facilitate observation of ongoing reconfigurations of the states of the switches in the devices and, in particular, the analysis of avalanches (see below). Measurements over long time periods are also necessary to avoid significant cutoffs in the power law distributions (*17*). Here, we presented data from DC stimulus of four devices, but the data presented are consistent with that obtained from DC, pulsed, and ramped voltage stimulus of a further 10 devices.

Our electrical measurements were performed using two distinct sets of measurement electronics to allow measurement of the device conductance on two distinct time scales. The first method relies on a picoammeter and is limited to a relatively slow sampling rate (0.1 s sampling interval). The second method uses a fast digital oscilloscope to allow a much higher sampling rate (200 μs sampling interval for the data presented here). As shown in Figs. 2 to 4, both methods resulted in qualitatively and quantitatively similar data, with similar power law exponents for each of the main quantities of interest. Hence, our results and conclusions are not influenced by the sampling rate.

### Data analysis

The data analysis methods used in this work are substantially the same as those developed in the neuroscience community to analyze microelectrode array recordings from biological brain tissue (*24*–*27*). Events are defined as changes in the conductance signal that exceed a threshold value. We show in fig. S3 that, as in the neuroscience literature, the choice of threshold in a reasonable range does not significantly affect the avalanche analysis. Avalanches were defined following the procedure described in (*24*, *45*): The data were analyzed by defining time bins with width equal to the mean IEI (*24*), and the duration of the avalanche was defined as the number of consecutive nonempty bins between empty bins. Figure S4 shows the effect of bin size on the avalanche analysis.

### Power law exponents and “fitting” processes

We followed the maximum likelihood (ML) approach of (*46*), as implemented in the Neural Complexity and Criticality MATLAB toolbox by (*45*), to statistically validate criticality. For all data that can be plotted in the form of a PDF (IEIs, *S*, and *T*), the ML estimators were obtained for both power law and exponential distributions. We used the Bayesian information criterion (*47*) to identify which distribution is more likely and find in all cases that it is the power law (*46*), as shown in fig. S5.

We generated 500 random power laws using the calculated ML estimators and iterated the choice of cutoffs (*x*_{min} and *x*_{max}) for the data range, finding the probability of obtaining a Kolmogorov-Smirnov (KS) distance (*46*) at least as great as for that of experimental data. In all cases, we failed to reject the null hypothesis that avalanche distributions are power law distributed (we required *P* > 0.2). Conversely, we rejected the null hypothesis that avalanche distributions are power law distributed for shuffled data (see Fig. 3, A, B, F, and G) (*41*), i.e., when we removed the correlations in the experimental avalanche data, we found that the distributions are exponential: The correlations and avalanches are intimately linked.

The KS test is extremely sensitive to small deviations from a mathematical power law (*45*), and since neuronal avalanche data are commonly observed to have curvature at small *x* values (*48*), a minimum *x* value was introduced (*46*). The deviation from power law behavior in our data is rather small (see Figs. 2 and 3), but we nevertheless followed (*45*, *46*) and allowed the data to be truncated to show definitively that the distributions are power law. The IEI distributions were found to be power law over about two orders of magnitude in time, with both a lower and upper cutoff. The avalanche size and duration distributions are always found to be power law up to the maximum range of the data, i.e., the data are power law over at least three orders of magnitude for the sizes and two for the durations. In all plots shown here, the range in which the power law model is in agreement with the data is indicated by the solid line. The estimated uncertainties in α and τ are ±0.1. Uncertainties estimated from the ML procedure are smaller but do not reflect the scatter in the data.

The ML methods cannot be applied to data which is not in the form of a probability distribution, and so, as in (*23*, *25*, *45*), standard linear regression techniques were used to obtain the measured exponents for 〈*S*〉(*T*) and *A*(*t*).

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/11/eaaw8438/DC1

Fig. S1. Illustration of different types of two-dimensional (2D) percolation model.

Fig. S2. Illustration of 2D percolation with tunneling and atomic switches.

Fig. S3. Choice of threshold.

Fig. S4. Choice of time bin size.

Fig. S5. Comparison of power law and other fits to the avalanche size and duration distributions.

Fig. S6. Scatter plot of α(τ).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

**Acknowledgments: **We thank D. O’Neale, M. Plank, and D. Stouffer for helpful discussions and support with ML estimation and associated methods, and J. Beggs, N. McNaughton, and C. Young for discussions on neuronal avalanches and criticality. **Funding:** This project was financially supported by The MacDiarmid Institute for Advanced Materials and Nanotechnology, the Ministry of Business Innovation and Employment, and the Marsden Fund UOC1604. **Author contributions:** S.A.B. conceived the study and initiated the project. J.B.M. initiated and performed the avalanche and criticality analysis. J.B.M., S.S., S.K.B., S.K.A., and E.G. performed the experiments. J.B.M., S.S., S.K.A., S.K.B., and S.A.B. performed the data analysis. S.A.B. and J.B.M. wrote the manuscript with comments from all authors. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).