Avalanches and edge-of-chaos learning in neuromorphic nanowire networks

The brain’s efficient information processing is enabled by the interplay between its neuro-synaptic elements and complex network structure. This work reports on the neuromorphic dynamics of nanowire networks (NWNs), a unique brain-inspired system with synapse-like memristive junctions embedded within a recurrent neural network-like structure. Simulation and experiment elucidate how collective memristive switching gives rise to long-range transport pathways, drastically altering the network’s global state via a discontinuous phase transition. The spatio-temporal properties of switching dynamics are found to be consistent with avalanches displaying power-law size and life-time distributions, with exponents obeying the crackling noise relationship, thus satisfying criteria for criticality, as observed in cortical neuronal cultures. Furthermore, NWNs adaptively respond to time varying stimuli, exhibiting diverse dynamics tunable from order to chaos. Dynamical states at the edge-of-chaos are found to optimise information processing for increasingly complex learning tasks. Overall, these results reveal a rich repertoire of emergent, collective neural-like dynamics in NWNs, thus demonstrating the potential for a neuromorphic advantage in information processing.

T he age of big data has driven a renaissance in Artificial Intelligence (AI). Indeed, the ability of AI to find patterns in big data could arguably be described as super-human; the human brain is simply not designed for brute-force iterative optimisation at scale. Rather, the brain excels at processing information that is sparse, complex and changing dynamically in time. The increasing prevalence of streaming data requires a shift in neuro-inspired information processing paradigms, beyond static Artificial Neural Network (ANN) models used in AI. The brain's unique capacity for adaptive, real-time learning is enabled by the complex interplay between its neuro-synaptic non-linear elements and a recurrent network topology 1 , with information processing manifested through emergent collective dynamics 2 .
Physically implemented neuro-inspired information processing has been demonstrated by nano-electronic device components that integrate memory and computation, enabling a shift away from the traditional von Neumann architecture. Resistive switching memories 3 are an important class of such devices known as memristors 4 . Their electrical response depends not only on applied stimulus, but also on memory to past signals 5 , which can mimic the short-term plasticity and long-term potentiation of neural synapses 6 . Utilising memristive switches as artificial synapses in crossbar architectures 7 has shown great promise in realising physical implementations of popular ANNs such as feed-forward 8 and convolutional neural networks 9 .
Self-assembled networks of metallic nanowires with memristive cross-point junctions [10][11][12][13] are an approach to move from rigid crossbar array architectures towards a more neural-like structure in hardware. Nanowire networks (NWNs) exhibit a small-world topology 14 , thought to be integral to the brain's own efficient processing ability 15 . Electrochemical metallic filament growth through electrically insulating, ionically conducting coatings (e.g. metal-oxides 16 , Ag 2 S 10 or PVP 11 ), enables memristive switching at the metal-insulator-metal junctions 11,17 . The interplay between memristive switching and a recurrent network topology promotes emergence of collective, adaptive dynamics, such as formation of conducting transport pathways that can be dynamically tuned 18,19 , giving networks a biologically plausible structural plasticity 20 . These neuromorphic properties equip NWNs with unique learning potential, with applications ranging from shortest-path optimisation 21 to associative memory 22 . Additionally, the neuromorphic dynamics of NWNs may be exploited for processing dynamic data in a reservoir computing framework 23 , as shown in both experimental 24,25 and simulation [26][27][28] studies.
It has been widely postulated that optimal information processing in non-linear dynamical systems may be achieved close to a phase transition, in a state known as criticality 29 . Distinct phase transitions are associated with criticality, notably 'avalanche criticality' 30 and 'edge-of-chaos criticality' 31 . In avalanche criticality, a system lies at the critical point of an activity-propagation transition where perturbations to the system may trigger cascades over a range of sizes and duration, characterised by scale-free power-law distributions. In sub-critical systems, activity can only propagate locally. Super-critical systems exhibit characteristically large avalanches that span the system. Scale-invariant avalanches, concomitant with avalanche criticality, have been observed in neuronal cultures 29,32,33 and neuromorphic systems comprised of percolating nanoparticles 34 . In edge-of-chaos criticality, dynamical states lie between order and disorder and the system retains infinite memory to perturbations. Edge-of-chaos dynamics have been observed in cortical networks 35 and appear to optimise computational performance in recurrent neural networks 36 , echo state networks 37 and random boolean networks 38 .
The observation of 1/f power spectra in neuromorphic NWNs has led to the suggestion that they may be poised at criticality 12,39 . While necessary, 1/f spectra are insufficient for criticality as 1/f noise can be produced by a diverse range of processes without spatio-temporal correlations, including uncoupled networks of isolated memristors 40 . Here, we present evidence for avalanche criticality in memristive NWNs and show that a critical-like state occurs near a first-order (discontinuous) phase transition. We also present the first evidence for edge-of-chaos criticality in neuromorphic NWNs and demonstrate that information processing is optimised at the edge-of-chaos for computationally complex tasks. Our results reveal new insights into neuroinspired learning, suggesting that in addition to non-linear neuromemristive junctions, the adaptive collective dynamics facilitated by neuromorphic network structure is essential for emergent brain-like functionality.

Results
A physically motivated model for neuromorphic structure and function in nanowire networks. PVP-coated Ag nanowires selfassemble to form a highly disordered, complex network topology (experiment: Fig. 1a; simulation: Supplementary Fig. 1). As a neuromorphic device, NWNs are operated by applying an electrical bias between fixed electrode locations across the network 12,22 . To gain deeper insight into the neuromorphic dynamics, a physically motivated computational model of Ag-PVP NWNs was developed. The model is briefly described here (see Methods for further details). b Biased nanowire-nanowire junctions promote Ag + migration and redox-induced Ag filament formation through the electrically insulating, ionically conducting PVP layer. Filament gap distance s varies from s max to 0. Junction resistance is modelled as a constant series resistance ( and time-dependent tunnelling resistance (R t ). c Junction conductance (G jn ¼ R À1 jn ) and s, as a function of state variable Λ(t). As |Λ| increases, the junction transitions from high resistance to tunnelling and then ballistic transport.
Ag|PVP|Ag junctions lie at the intersection between nanowires. These electrical junctions are modelled as voltage-controlled memristors 12,41 . Λ(t) represents the state of a junction's conductive Ag filament in the PVP layer (Fig. 1b, Eq. (5)), and is closely related to the physical filament gap distance (Fig. 1b, c, Eq. (8)) that determines the junction's conductance (G jn ðtÞ ¼ R À1 jn ¼ I jn =V jn ). A junction filament grows when voltage exceeds a threshold: |V jn | > V set . The sign of Λ(t) encodes the junction polarity: a filament approaches Λ ¼ Λ max when V jn > V set and approaches Λ ¼ ÀΛ max when V jn < − V set . When the voltage is sufficiently low (|V jn | < V reset ), the filament decays towards Λ = 0, as a result of thermodynamic instability. The rate of filament growth or decay is determined by the voltage difference relative to V set or V reset . Hence, Λ(t) encodes a junction's memory to past electrical stimuli. See Methods (Eq. (5), Fig. 9) for full mathematical details. Figure 1c shows the non-linear dependence of G jn on |Λ| that produces switch-like junction dynamics. When 0 ≤ |Λ| < Λ crit , the junction is insulating (G jn = G off ≪ G on , where G on = G 0 is the conductance quantum). As |Λ| approaches Λ crit , the junction transitions to a tunnelling regime where conductance exponentially grows as |Λ| increases. After the filament has grown (Λ crit ≤ jΛj < Λ max ), the filament tip width is the limiting factor of conductance and transport is ballistic, with G jn = G on = G 0 .
Next, simulations using this model are presented analysing the network level dynamics of this neuromorphic system.
Adaptive network dynamics near a discontinuous phase transition. NWNs driven by a constant bias V converge to a steadystate conductance (simulation: Fig. 2; experiment: Supplementary  Fig. 2a). The transient dynamics of the network conductance time-series, G nw (t), occurs in steps, as the network shifts between metastable states (plateaus in G nw ). Above a threshold voltage V th , NWNs exhibit system-level switching, with G nw dramatically increasing by three orders of magnitude. Post-activation, the conductance increase slows to a rate comparable with preactivation (on logarithmic scale). Increasing V increases the activation rate and diminishes the step-like increases in conductance. These qualitative features of G nw (t) are independent of the specific network and voltage in simulation and experiment.
Simulations reveal that in the Low Conductance State (LCS: G nw ≲ 10 −3 G 0 for this network), all network junctions exhibit low tunnelling conductance and no junctions are in a ballistic transport regime. The High Conductance State (HCS: G nw ≳ 10 −1 G 0 ) exhibits pathways of high-conductance junctions (G jn ≈ G 0 ) spanning from source to drain electrodes. Figure 2b shows a single pathway, however the HCS may also consist of multiple, parallel high conductance pathways at higher voltages (Supplementary Fig. 3), and for large or dense networks with many equivalently short source-drain paths. Figure 2c shows the steady-state G nw (re-scaled to G Ã 1 ¼ nG nw ðt ¼ 1Þ=G 0 , where n is the graph length of the shortest path between source and drain electrodes) as a function of V * = V/V th for two different network initial states: an inactive (LCS) network, with all junctions initialised to Λ = 0 (black circles); and a network pre-activated (HCS) to G Ã nw ¼ nG nw =G 0 % 4:7 by a 1.8 V = 20V th DC bias for 10 s (red triangles). For V * < nV reset , networks always converge to a stable LCS, irrespective of initial conditions. Above this level, multiple stable G Ã 1 states emerge. Multi-stable conductance states for a wider range of initial conditions are shown in Supplementary Fig. 4. Networks exhibit hysteresis, since initial junction filament states control which stable state is reached. In simulation, for all initial conditions, G Ã 1 is always non-decreasing with V * , but this is not observed experimentally ( Supplementary Fig. 2c) as it is difficult to prepare a network in exactly the same state before each activation (due to long-term memory of junctions). Additionally, persistent fluctuations from junction noise, not included in the model, can facilitate transitions between multi-stable states 12,39 ( Supplementary Fig. 2b).
A discontinuity in the order parameter (steady-steady conductance, G Ã 1 ) between LCS (G Ã 1 $ 0) and HCS (G Ã 1 $ 1) is revealed as the control parameter (V) is varied for any given network and initial conditions. This is found in both simulation (e.g. V * = 0.5 for pre-activated and V * = 1 for inactive networks in Fig. 2c) and experiment ( Supplementary Fig. 2c). These distinctive discontinuities (V * = 0.5 for pre-activated and V * = 1 for inactive), along with the properties of multistability and hysteresis, indicate that the formation/annihilation of the first conducting pathway in the network corresponds to a first-order (discontinuous) non-equilibrium phase transition, marking an abrupt change in the global state of the system. In the model, the location of the discontinuous transition is universal (Supplementary Fig. 5), as it depends only on the length of the shortest path between the electrodes and on the initial state (either LCS or HCS). Networks in a LCS (no conducting pathways) transition to NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-24260-z ARTICLE a HCS (conducting pathways exist) when V * > nV set . Networks in a HCS transition to a LCS when V * < nV reset . In real (i.e. nonideal) NWNs, V set and V reset may vary between junctions due to variations in the thickness of PVP coating around nanowires.
Additional discontinuities in G Ã 1 may be observed at higher integer multiples of thresholds (V = mV set , V * = kV reset ). For example, at V * = 11/9 ≈ 1.2 (corresponding to V = 11V set ) in Fig. 2c. These correspond to the formation/annihilation of additional transport pathways. In the limit of large network size, only the first discontinuity (between G Ã 1 $ 0 and G Ã 1 $ 1) in G Ã 1 vs V remains concomitant with formation of the first conducting pathway ( Supplementary Fig. 5).
These results demonstrate that NWNs can adaptively respond to external driving and can undergo a first-order phase transition between bi-stable states (LCS and HCS). These global network dynamical states arise from the recurrent connectivity between junctions and their switching dynamics, as shown next.
Collective junction switching drives non-local transport. Network activation or de-activation as described above can be understood as a collective effect emerging from recurrent connections between junctions. By Kirchoff's Voltage Law (KVL), the sum of voltages along any path between source and drain is equal to the externally applied voltage, while the sum of voltages around any closed loop in the network is zero. Since the voltage across a given junction V jn (t) controls the future evolution of its filament state Λ(t) and hence conductance G jn (t), KVL couples the dynamics of memristive junctions to the network topology. A simulation at V * = 1.1, with uniform initial conditions (Λ = 0), showing G jn (t) (Fig. 3a), V jn (t) (Fig. 3b) and corresponding visualisations at three time-points ( Fig. 3c-e), reveals the qualitative characteristics of collective switching dynamics in NWNs.
Initially, the disordered network structure ensures the voltage distribution ( Fig. 3b) is inhomogeneous, with junctions near the electrodes typically having the highest V jn . For junctions with |V jn | < V set , filament states remain static, or otherwise decay to |Λ| = 0, if Λ ≠ 0 and |V jn | < V reset . When |V jn | > V set , filaments grow at a constant rate dΛ/dt = |V jn | − V set until G jn (Fig. 3a) rapidly increases from the onset of tunnelling, coinciding with a decrease in V jn (Fig. 3b). G jn plateaus once |V jn | reaches V set . This first occurs at t ≈ 0.8 s for junction #1. Voltage is redistributed locally to neighbouring junctions 42 (junctions near opposite electrodes, e.g. #1 and #9, can be considered connected by an edge, with R jn = 0, at the power source). If the subsequent voltage increase results in |V jn | > V set , the filament growth rate increases for these junctions. For junctions already in a tunnelling regime (G off < G jn < G on = G 0 ), G jn briefly increases until V jn returns to V set , whereas insulating junctions experience a delay before G jn increases as their filament grows. Hence, G nw displays small, steplike increases coinciding with junctions transitioning from insulating (G jn = G off ) to tunnelling (G off < G jn < G on ) and groups of already active junctions adaptively adjust their G jn such that their V jn lies on or below the threshold. These active junctions effectively wait for other junctions to activate before switching on, resulting in a switching synchronisation phenomenon similar to that seen in 1D memristive networks 43,44 .
Through these switching events, the network self-organises to prevent a large voltage drop (greater than V set ) across any given junction, funnelling most of the current down a single (or a few) transport pathways, which grow from the electrodes towards the centre of the network (Fig. 3c, d). For any given pathway of length m, if V * < m V set , then after higher voltage junctions adjust their voltage to V set the pathway ceases to grow. This is the fate of all paths in networks with V * < 1. When, V * > 1, nearby pathways grow in competition, with the shortest pathway forming first.
Once the final junction along a source-drain pathway enters a tunnelling regime (t ≈ 6 s), G jn of each junction on the path adjusts ( Fig. 3e) such that they receive the same voltage (V jn = V/n > V set ), since a larger V jn increases dΛ/dt and hence, G jn , which in turn reduces V jn . At this point, junctions along the path behave collectively like a single memristive junction under constant bias: for each of these junctions, Λ grows linearly in time, leading to exponential growth in G jn and hence, G nw . This manifests itself as a large, rapid increase in G nw (t). This growth ceases when all junctions along the path reach the ballistic transport regime (Λ = Λ crit , G jn = G 0 ), where they stably remain at later times. For higher voltages or denser networks, additional transport pathways may form, resulting in later steps in G nw (cf. Fig. 2a). For sufficiently high voltages (when V jn > V set for all junctions along the path), the plateaus become less pronounced and conductance increases continuously with time, however synchronous switching is still observed (Supplementary Fig. 6).
These results show that the emergence of transport pathways can be attributed to coupling between the complex network topology and memristive junction switching. Cascades of activity are induced as junctions transition into the conducting regime, adaptively redistributing voltage to their neighbours. This collective switching activity near a phase transition is reminiscent of avalanche dynamics, investigated next.
Avalanche switching dynamics. Avalanches with scale-free size and lifetime event statistics-a hallmark of critical dynamicshave been observed in neuronal populations 29,32,33 and other neuromorphic systems 34,45 . For NWNs, the conductance timeseries of each junction (from the model described above) is converted to discrete switching events by introducing a threshold for the conductance change ΔG jn between adjacent time-points. A similar procedure is applied to experimental measurements of G nw . As in neuronal data 32,33 , events are binned discretely in time into frames (width Δt). This allows an 'avalanche' to be defined as a sequence of frames containing events preceded and followed by an empty frame. Avalanche size (S) is defined as the total number of events in the avalanche. Avalanche life-time (T) is defined as the number of frames in the avalanche. See Methods for full details.
By altering the strength of the driving voltage away from the threshold V th , the avalanche distributions begin to deviate from a power-law (Fig. 5). When V * < 1, pathways are unable to form across the network and the switching events result in small-scale avalanches (black points). As V * approaches 1, the distribution elongates and becomes a power-law (red points). Above V * = 1 (when the networks activate), a bi-modal distribution is evident, as avalanches of large characteristic sizes and lifetimes emerge above the power-law tail (cyan and blue points). As network size increases, the probability density of the bump relative to the power-law region grows ( Supplementary Fig. 8d-f). This suggests these anomalously large avalanches are consistent with supercritical states.
Order-chaos transition from polarity-driven switching. The response of a NWN to constant voltage stimulus, examined above, has revealed rich collective dynamics, but ultimately conductance converges to a steady state. When networks are driven by unipolar periodic stimuli, they converge to periodic attractor states. However, in response to a periodic, alternating polarity driving signal, NWNs exhibit a more diverse range of dynamics, from ordered to chaotic, depending on the amplitude A and frequency f of the driving signal. This is quantified by calculating the maximal Lyapunov exponent λ for simulated networks driven by a triangular AC signal (see Methods).   increasing A results in an increase in λ towards zero at A ≈ 3 V. Above this frequency, increasing A can lead to a λ > 0 and the onset of chaotic dynamics in the network. For A ≳ 2.5 V and f ≳ 1.5 Hz, stronger chaotic dynamics are inferred, with λ ≳ 10 s −1 (yellow region). Further increasing f for a given A leads to a decrease of λ. For a given λ, the collective switching states of NWNs was characterised by the ratio of maximum to minimum network conductance, r, averaged over each AC cycle. When 1 ≲ r ≲ 10, the network remains in either LCS (no pathways) or HCS (pathways exist). When r ≳ 10 3 , the network is persistently driven back and forth across the transition between these states. Figure 6b reveals that only networks exhibiting strong switching (r ≳ 10 3 ) can exhibit chaotic dynamics (λ > 0), whereas both edgeof-chaos dynamics (λ ≈ 0) and ordered dynamics (λ < 0) are observed over the range of possible collective switching states.  Supplementary  Fig. 10. In the λ < 0 regime, I − V (Fig. 7a) and G − V (Fig. 7d) cycles exhibit symmetric, repetitive hysteresis. For these states the minimum conductance occurs when V = 0, meaning the network fully activates and de-activates before polarity reverses. In the edge-of-chaos regime (λ ≈ 0), I − V (Fig. 7b) and G − V (Fig. 7e) cycles are still symmetric and repetitive, but now the network does not deactivate completely before polarity reverses. In the chaotic regime (λ > 0), I − V (Fig. 7c) and G − V (Fig. 7f) trajectories diverge over successive cycles (evident from the thickening of blue curves). These chaotic trajectories are not unbound, but are instead confined to a reduced region of phase space (a 'chaotic attractor'). Examples of chaotic trajectories of different Lyapunov exponents are shown in Supplementary  Fig. 11.
Qualitatively, the observed Lyapunov exponents can be understood by considering a small perturbation (δΛ) to the filament state (Λ i ) of the i-th junction. When tunnelling is absent (cf. Fig. 1c), δΛ changes G jn and hence, leaves the voltage distribution of the network unchanged. Thus, as the network continues to evolve, δΛ is remembered by the junction (the junction states remain a fixed distance δΛ apart). When Λ i approaches its boundaries (cf. Fig. 1), either from when the filament decays to zero (Λ i → 0 when |V i | < V reset ) or saturates (jΛ i j ! Λ max when |V| > V set ), δΛ shrinks. When a junction is in a tunnelling regime, however, any perturbation to Λ i is exponentially amplified in terms of conductance. Perturbations that increase a junction's conductance decrease its voltage, slowing filament growth if |V i | > V set , or increasing filament decay if |V i | < V reset . The effect on neighbouring junctions is the opposite. Under slow driving, the network can adapt to such perturbations and retain the size of the perturbation. Under fast driving, however, the network is unable to adapt and perturbations grow during activation and de-activation, leading to separation of nearby network states. The frequency which constitutes fast or slow depends on the amplitude applied (cf. Fig. 6) and the network structure (size and density). The dynamical balance between the mechanisms that enforce order (perturbations shrink) and create chaos (perturbations grow) determines the stability of the global network dynamics. Hence, tuning the driving signal allows control over system dynamics. As shown next, this may be advantageous when utilising neuromorphic NWNs for information processing 36,45 .
Information processing optimised at the edge-of-chaos. The information processing capacity of neuromorphic NWNs in different dynamical regimes was tested with the benchmark reservoir computing 51,52 task of non-linear wave-form transformation previously demonstrated in NWN experiments 24,25 and simulations 26,27 . In this task a triangular wave is input into the network, nanowire voltages are used as readouts and are trained using linear regression to different target wave-forms. Examples of target waveforms obtained are shown in Supplementary Fig. 12 and sample network readouts (nanowire voltages) are shown in Supplementary Fig. 13. Figure 8 shows the network performance in transforming an initial triangular wave of given A and f to different target waveforms as a function of the maximal Lyapunov exponent λ for each  Fig. 6). Over the range of λ, the rank order of relative accuracy between tasks is consistent with the similarity between the input (triangular) and target waveforms, namely sin > square > π/2-shift ≈ 2f. The similarity between waveforms can be understood by considering the Fourier decomposition of each signal. To convert to a sine wave, higher harmonics must be removed from the triangular input signal, whereas conversion to a square wave requires additional odd higher harmonics. For double frequency conversion, odd higher harmonics must be removed, with even harmonics added. For phase shift conversion, the network must produce a lag to the input signal, i.e. coefficients of cosine terms in the Fourier series become coefficients of sine terms. The sine wave transformation accuracy is always ≥ 0.98 before decreasing to ≈ 0.95 above λ ≈ 0. The square wave and π/2shifted wave accuracy increase as λ increases towards zero, peaking when −2 s −1 ≲ λ ≲ 0 s −1 at 0.86 for square wave and 0.69 for π/2-shifted wave. Like the sine wave, both the square wave and π/2-shifted wave accuracy tapers off rapidly when the network is in a chaotic regime (λ > 0), to an approximately constant accuracy of~0.6 for square and 0.3 for π/2-shift.
Notably, for the π/2-shifted wave, the accuracy for chaotic states lies above the minimum accuracy for ordered states. For the 2f target wave, the accuracy is zero for λ ≲ − 0.5 s −1 , before rapidly increasing to peak near the edge-of chaos regime λ ≈ − 0.1 s −1 at an accuracy of 0.79. Like the other wave-forms, accuracy decreases as the network becomes chaotic (λ > 0), but accuracy is higher than for ordered dynamics (λ < 0). For all tasks, strongly chaotic states underperform compared to the edge-of-chaos. This suggests that while more ordered (sine wave target) or slightly chaotic (2f target) dynamical regimes may be optimal depending on the computational complexity of the task, only a relatively   narrow range of network states, those tuned to the 'edge-ofchaos', are robust performers over a diverse set of target waveforms.

Discussion
The formation (or destruction) of long-range transport pathways between electrodes is a ubiquitous feature of disordered memristive networks with threshold-driven junction switching 53 , including nanowires 12,18,19 and nanoparticles 54 . Here, it was found that the network undergoes a discontinuous phase transition when the first pathway forms, with a linear relation between threshold voltage and source-drain path length. A similar discontinuous transition coinciding with pathway formation was reported in simulations of adiabatically driven memristive networks with large G on /G off ratios (ratio between maximum and minimum possible G jn ) 53 . Smaller G on /G off ratios were instead found to result in a continuous transition. When G on /G off is large, changes to G jn lead to large changes in the voltage distribution, strongly coupling the dynamics of junctions in the network. This results in a switching synchronisation phenomenon (cf. Fig. 3a), analogous to that observed in 1D memristive networks 43,44 .
Conversely, when G on /G off is small, junctions changing conductance have a smaller effect on their neighbours. In this regime, junctions still collectively switch between metastable states, but do not exhibit switching synchronisation ( Supplementary Fig. 14).
Hence, by tuning G on /G off , which can be achieved experimentally by varying thickness of PVP coating, the collective behaviour of NWNs can be controlled. Finding power-law avalanche distributions over a few orders of magnitude near a first-order (discontinuous) transition is surprising, as they are usually associated with second-order (continuous) transitions 55 . Disorder and hysteretic behaviour are likely key ingredients for scale-free avalanches observed in other systems driven near a discontinuous phase transition, such as magnetic materials 56 and neural networks 57 . These properties are certainly present in NWNs.
Critical avalanches were demonstrated here by tuning voltage to the activation threshold, but the role of self-organisation in achieving and sustaining critical avalanches requires further investigation. Experiments show NWNs continue to exhibit persistent, bi-directional conductance fluctuations (Supplementary Fig. 2b) as a result of junction noise, electromigrationinduced junction breakdown events, and subsequent recurrent feedback by the network 12,58 . These mechanisms may allow NWNs to self-organise to an avalanche critical state, allowing mapping to models such as 'self-organised criticality' (SOC) 30 (self-organisation to a continuous transition), 'self-organised bistability' (SOB) 59 (self-organisation to a discontinuous transition), quasicriticality 60 (departs from criticality with crackling noise equation (Eq. (4)) obeyed) and their non-conservative counterparts 61 . Over long time scales, if anomalously large avalanches coexist with scale-free avalanches, then SOB would be a better description of NWNs than SOC. However, fluctuations could plausibly make the transition (cf. Fig. 2c) continuous resulting in a SOC-like state.
Critical dynamics has previously been observed in selfassembled tin nanoparticle networks (NPNs) 34,62 . There are notable differences between dynamics observed in nanowire and nanoparticle networks. In NWNs, resistive switching is facilitated by filament growth through an insulating layer. In NPNs, the nanoparticles are not coated with an insulating layer and resistive switching is due to tunelling-driven filament growth across nanogaps between nanoparticles. Thus, in the absence of filament growth, nanoparticles in contact are conductive, whereas nanowires in contact are insulating. Consequently, resistive switching dynamics and critical avalanches are only observed when the nanoparticle density is finely tuned to the percolation threshold. NWNs, on the other hand, exhibit resistive switching and avalanches at densities close to and well above the percolation threshold, such as twice the threshold (cf. Supplementary Figs. 7  and 8). Conversely, in NPNs, the voltage does not need to be finely tuned to achieve critical avalanches, provided networks are on the percolation threshold. Breakage of conductive filaments from Joule heating and electromigration effects self-tune NPNs to a critical state. In NWNs, critical avalanches with power-law sizes and life-times are observed when tuning networks close to the threshold voltage. At voltages below the threshold, avalanches do not span the network and exhibit exponentially decaying avalanche distributions (cf. Fig. 5), while at voltages significantly above the threshold, large avalanches of a characteristic size and duration are observed, corresponding to formation and annihilation of non-local conducting pathways. Thus, passivation of nanoscale metallic components affords the advantage of not having to fine-tune density.
The universality of avalanches poses an interesting question for future studies on neuromorphic systems. The experimental studies of avalanches in percolating NPNs presented non-universal avalanche exponents that satisfy criteria for avalanche criticality 34,62 . Notably their model 62 exhibits avalanche exponents (τ = 2.0, α = 2.3, 1/στz = 1.3) very close to NWN values, suggesting they may belong to the same universality class. A capacitive breakdown model of current-controlled NWNs, in a lower current regime than studied here, 63 found universal avalanche exponents consistent with the classic SOC sandpile model 30 . However, in other neuromorphic systems such as adiabatically driven memristive networks 53 and spiking neuromorphic networks 45 , as well as neuronal culture experiments 32 , avalanche statistics match that of a branching process 64 , a member of the universality class of directed percolation (τ = 1.5, α = 2, 1/στz = 2).
While ordered attractor (λ < 0) states in models of networks of voltage-controlled memristors under alternating polarity stimuli have previously been observed 65 , the diverse edge-of-chaos and chaotic dynamics of neuromorphic NWNs have not been previously shown in memristive networks. Unlike certain types of memristors 66 , individual junctions driven by periodic stimuli are incapable of exhibiting chaos, but chaos emerges in these networks due to strong recurrent coupling between components 67 . This requires an alternating polarity pulse where activation and de-activation are both strongly driven: for unipolar periodic pulses only λ ≤ 0 is found. As the 'edge-of-chaos' (λ ≈ 0) can be reached in a strongly driven regime (high f), it does not necessarily coexist with 'avalanche critical' states. Avalanche distributions (P(S) and P(T)) near λ ≈ 0 do not follow power-laws ( Supplementary Fig. 15): power-law fits fail Kolmogorov-Smirnov test unless range of fit is made very small (x max =x min ≲ 2), unlike the DC case at V * = 1. This deviation may be attributed to the fast driving signal which ensures activity is injected into the network at a non-uniform rate while avalanches propagate, breaking time-scale separation between driving and network feedback (avalanches), obfuscating the distinction between consecutive avalanches. This reinforces the point often overlooked in neural network models that despite often coinciding 68 , activity propagation (avalanche) and order-chaos transitions are distinct 69 .
The observation of optimal overall performance on the nonlinear transformation task at the 'edge-of-chaos' corroborates the popular hypothesis of robustness of information processing near phase transitions 31 and is consistent with simulations in other types of recurrent networks [36][37][38] . Despite this, the task dependence of accuracy is striking. For the simplest task (transformation to sine wave) the 'edge-of-chaos' state afforded no computational advantage. On the other hand, the greatest payoff from the 'edge-of-chaos' state was found for the most dissimilar target wave-forms (double frequency, phase shifted). This result corroborates a previous study using a spiking neuromorphic network that found critical dynamics maximises the abstract properties of the system (auto-correlation time, susceptibility and information theoretic measures) and hence, performance in tasks of non-trivial computational complexity, yet for simpler tasks ordered dynamical states (away from criticality) may perform more optimally 45 .
Neuromorphic NWNs may be utilised for a range of information processing tasks. Information may be stored in memristive junction pathways for static tasks such as associative memory 22 and image classification 26 . However, it is the coupling of memristive junctions with recurrent network topology that makes NWNs ideal for temporal information processing when implemented in a reservoir computing framework, such as signal transformation 24,25,27 and non-linear time-series forecasting 26,70 . These applications suggest NWNs are a promising neuromorphic system for adaptive signal processing of streaming data.
The rich dynamical behaviour revealed here may be observed on other network topologies, provided networks are highly recurrent and disordered. Recurrent networks have many short loops 71 , allowing junctions to be strongly coupled (e.g. short loops in Fig. 2b which are coupled by Kirchoff's laws), thus generating the diverse range of time scales for avalanche events and chaotic dynamics to be observable. This recurrent topology is important for information processing in networks, such as recurrent neural networks 52,71 . Disorder ensures spatiotemporal inhomogeneity of the voltage distribution, maximising the degrees of freedom of networks. Understanding how to harness network structure for optimal information processing provides an exciting future challenge for neuromorphic network research.
In conclusion, neuromorphic NWNs respond adaptively and collectively to electrical stimuli and undergo a first-order phase transition in conductance at a threshold voltage determined by the shortest path length between electrodes. Due to recurrent loops in the network, junctions switch collectively as avalanches. These avalanches are consistent with avalanche criticality at the critical voltage. At higher voltages, anomalously large avalanches coexist with avalanches spanning a range of scales. Under alternating polarity stimuli, networks can be tuned between ordered and chaotic dynamical regimes. The edge-of-chaos is the most robust dynamical regime for information processing over a range of task complexities. These results suggest that neuromorphic NWNs can be tuned into regimes with diverse, brain-like collective dynamics 2,35 , which may be leveraged to optimise information processing.

Methods
Simulations. NWNs were modelled ( Supplementary Fig. 1a) by randomly placing lines of varying length (randomly sampled from a gamma distribution with mean 10 μm and standard deviation 1 μm) and angular distribution (sampled uniformly on [0, π]) on an 30 × 30 μm 2 grid (centre locations sampled from a uniform distribution).
The network was converted to a graph representation ( Supplementary Fig. 1b) in which nodes and edges correspond to equipotential nanowires and Ag|PVP|Ag junctions, respectively. The largest connected component was used in the simulations. Unless stated otherwise, a network with 100 nanowires and 261 junctions was used.
A voltage bias is applied between source and drain nodes (chosen at opposite ends of the network) and the model solves Kirchoff's laws to calculate the voltages V jn (t) across each junction at each time point. At each junction (edge), electrochemical metallisation is phenomenologically modelled with a conductive filament parameter (Λ(t)). The filament parameter is restricted to the interval ÀΛ max ≤ ΛðtÞ ≤ Λ max and dynamically evolves according to a threshold-driven bipolar memristive switch model (Fig. 9, Eq. (5)) 41,72 : Junction conductance (cf. Fig. 1b, c) is modelled as a constant residual resistance of the insulating PVP layer (R r ), in parallel with constant filamentary resistance (R f ¼ G À1 0 ( R r ) and Λ-dependent tunnelling resistance (R t ), with tunnelling conductance G t calculated using the low voltage Simmon's formula (Eq. (7)) for a MIM junction 73 .
with effective mass m * = 0.99m e and PVP layer (assumed homogeneous) thickness s max ¼ 5 nm . The potential barrier ϕ = 0.82 eV is the difference between Fermi levels of PVP and Ag. A = (0.41 nm) 2 = 0.17 nm 2 is the area of a face of the silver unit cell. Nanowire resistance is considered negligible compared to junction resistance 11 . G t thus introduces an additional non-linear dependence on V, through the filament growth parameter s = s(Λ(V)), that modulates junction switching due to filament formation (cf. Supplementary Fig. 20). Free parameters are chosen such that activation time is comparable to experiment. Values used are V set = 0.01 V, V reset = V set /2, Λ crit = 0.01 Vs, Λ max ¼ 0:015 Vs, R off ¼ 10 3 ðG 0 Þ À1 ¼ 12:9 MΩ, R f ¼ ðG 0 Þ À1 ¼ 12:9 kΩ and b = 10. For the Lyapunov analysis and non-linear transformation task V set = V reset was used. Simulations use the Euler method with time-step dt = 10 −3 s. The effect of model parameters on results presented here is discussed in Supplementary Information.
Experimental. Physical NWNs were synthesised with the polyol process 74 using 1,2-propyleneglycol (PG) as an oxidising agent for silver nitrate (AgNO 3 ) and were drop-cast onto a glass substrate 12 . NWN size was 500 × 500 μm. Nanowires had mean length 7.0 ± 2.4 μm, diameter 500 ± 100 nm and density ≈ 0.1 nw/μm 2 determined with a high amplification optical microscope. The PVP-coating thickness is 1.2 ± 0.5 nm 12 . Networks were electrically stimulated with pre-  5)). Arrows are proportional to filament growth rate (dΛ/dt) and point in the direction of filament growth.
patterned rectangular gold electrodes of width 500 μm and current was read out using an in-house amplification system 12 at a sampling rate of 1 kHz.
Avalanche analysis. For simulation, an event is defined as when 1 G jn j ΔG jn dt j exceeds a certain threshold (10 −3 s −1 ) before returning below the threshold. A similar procedure was applied to experimental data using the network level (rather than junction level) conductance time-series, with an event defined as when 1 G nw j ΔG nw dt j exceeds 1 s −1 before returning below the threshold, or when ΔG nw , exceeds a threshold 5 × 10 −8 S before returning below the threshold. Changing the event detection threshold leads to negligible change of avalanche statistics.
As in studies of avalanches in neural cultures 32 and nanoparticle networks 34 , the natural choice of frame (Δt) is the average inter-event-interval, 〈IEI〉, over the time-scale that events occur, where IEI is the time between adjacent events. Effect of changing frame on avalanche statistics is shown in Supplementary Fig. 16. Avalanche size and life-time distributions, P(S) and P(T), were binned linearly for simulations (logarithmically for experiments due to poorer resolution of the tail) and fit using a Maximum Likelihood (ML) power-law fit (P(x) ∝ x −σ ). Scaling exponent (σ) were determined to precision 10 −2 , lower (x min ) and upper (x max ) cut-offs 75 were determined to nearest integer. Statistical significance of the ML fit was tested using the Kolmogorov-Smirnov (KS) test 76 . 500 synthetic data sets were generated from the distribution of the ML fit. The KS distance between each and the ML fit, was compared with the KS distance between the empirical (simulation or experiment) data and the ML fit. The p value is the fraction of fits where the KS distance is smaller for the empirical data than the synthetic data. In all cases where power-laws were presented above, the null hypothesis that data followed a doubly truncated power-law was accepted to the chosen statistical significant (p > 0.5). For power-law fits with p > 0.5, no significant change was found to exponents (within uncertainties) for different x min and x max values up to a cut-off determined by system size (Supplementary Fig. 17). Hence, x min and x max are chosen such that log ðx max =x min Þ is maximised for fits with p > 0.5. As average avalanche size (〈S〉(T)) is not in the form of a probability distribution (so ML methods do not apply), it was fit using linear regression on a log-log plot.
In simulation, as networks converge to a steady-state under constant DC bias (cf. Fig. 2), to generate avalanches, either a voltage pulse is applied to a network in a steady state (allowing certain junctions to activate/deactivate, subsequently triggering avalanches of other switching junctions), or junctions are manually perturbed by changing their state (e.g. switching them from high G jn to low G jn or vice versa). Results based on the former method are presented here. Avalanches are calculated from the transients as networks relax to their steady state using the binning method described above. To obtain enough statistics to sample the avalanche statistical distributions an ensemble of 1000 randomly generated networks of physical dimensions 150 × 150 μm 2 (100 × 100 μm 2 for Fig. 5), density 0.10 nw (μm) −2 and with nanowire length 10 ± 1 μm were simulated for 30 s each starting from all filament states from Λ = 0. All nanowires with physical centre location within 2.5% (3.75 μm for this network size) of top of network were selected as source electrodes. All nanowires with physical centre location within 2.5% of bottom of network were selected as drain electrodes. This ensured the whole network is stimulated, reducing finite-size effects.
In experiment, noise and junction breakdown events 12 perturb networks from their steady-state allowing spontaneous initiation of avalanches at constant voltage. NWNs were stimulated with DC biases of gradually increasing voltage (with waiting time of 3 h between measurements) until the switching threshold was determined. The network was then stimulated with a voltage just above threshold for 1000 s, or until current exceeded 8 × 10 −5 A (to prevent network damage), recording data at 30 kHz. This was repeated 3 times (with 3 h wait between measurements) on the same network. Time-series are truncated, first by excluding initial times for which network current I < 1 × 10 −8 A and then between the first and last detected event. Avalanche statistics from all data-sets (at fixed voltage just above threshold) were combined to produce avalanche statistical distributions. Combining data-sets is valid under the assumption that avalanches at fixed voltage on an individual network follow the same statistical distribution.
Unlike in simulations, in experimental data, switching activity and avalanches in sub-threshold networks cannot be identified as they fall below the experimental noise floor. In contrast, junction conductance time-series can be simulated at all voltages. Examples of both experimental and simulated avalanches are shown in Supplementary Fig. 18.