Emergent dynamics of neuromorphic nanowire networks

Neuromorphic networks are formed by random self-assembly of silver nanowires. Silver nanowires are coated with a polymer layer after synthesis in which junctions between two nanowires act as resistive switches, often compared with neurosynapses. We analyze the role of single junction switching in the dynamical properties of the neuromorphic network. Network transitions to a high-conductance state under the application of a voltage bias higher than a threshold value. The stability and permanence of this state is studied by shifting the voltage bias in order to activate or deactivate the network. A model of the electrical network with atomic switches reproduces the relation between individual nanowire junctions switching events with current pathway formation or destruction. This relation is further manifested in changes in 1/f power-law scaling of the spectral distribution of current. The current fluctuations involved in this scaling shift are considered to arise from an essential equilibrium between formation, stochastic-mediated breakdown of individual nanowire-nanowire junctions and the onset of different current pathways that optimize power dissipation. This emergent dynamics shown by polymer-coated Ag nanowire networks places this system in the class of optimal transport networks, from which new fundamental parallels with neural dynamics and natural computing problem-solving can be drawn.

The human brain is a product of evolution, tuned and reshaped by an ever-changing environment. The brain's neuronal system is able to achieve the ability to recognize, conceptualize and memorize objects in the physical world. Using environmental information we establish logical associations that ultimately allows us not only to survive, but also to solve highly complex problems 1 .However, in an increasingly connected and interactive world, the volume of information to process has exponentially increased, and in order to extract and synthesize meaningful information, computerized approaches, such as machine learning and its various incarnations have gained tremendous popularity 2 .
Typically, Artificial Neural Networks (ANNs) attain this goal by a very delicate and case-selective combination of learning strategies 3 . Data containing complex or contextual associations between objects normally requires an heuristic sampling which limits their ability to synthesize information. Conventional CMOS architectures also restrains the amount of data that is efficiently processed with ANNs due to power consumption bottlenecks.
Interest in the creation of synthetic neurons that could increase the processing abilities of ANNs has increased considerably with the discovery of nanomaterials with memristive properties 4 . A memristive device is a non-linear two-terminal device in which the resistance shows resilience to change (i.e. memory), manifested in hysteretic behavior when the energy change is reversed or reduced, also termed as resistive switching. The memristor thus has two important neurosynapse-like properties, plasticity and retention. Traditional integrate-and-fire models, that emulate the electrical behavior of neurons using passive circuit elements, can be simulated exclusively with these elements [5][6][7] . Memristive devices have been successfully embedded into various CMOS architectures, enabling the realization of synthetic neural networks(SNN). SNNs imitate the topology of an ANN in a physical layout, typically stacking memristive terminals in cross-bar configurations 8,9 . Using voltage pulses to configure the internal state, or weight, of individual memristors; memorization, learning and classification abilities have been achieved [10][11][12][13] . However promising, this approach remains reliant upon CMOS technology

Structural properties of pVp-Ag nanowire networks
An optical micrograph image of a PVP-coated Ag nanowire network has been acquired after drop-cast deposition on a SiO 2 substrate. As is evident in Fig. 1a, nanowires are dispersed with random orientations and their spatial distribution exhibits areas of denser connectivity interspersed with sparser regions. Higher-resolution scanning electron microscopy (SEM) (Fig. 1b) provides a more detailed view of the nanowire network, revealing a high-degree of interconnectivity between adjacent nanowires, including multiple overlaps. Small amounts of particle-like or irregular Ag impurities not removed by centrifugation are observed as well. The average diameter and length of nanowires was measured to be 360 ± 110 nm and 14 ± 5 μm, respectively. The network density is thus estimated to be 0.76 nanowires/μm 2 , which is controlled by the concentration of nanowires in solution prior to drop-casting onto the surface. Upon further inspection with high-resolution transmission electron microscopy (HR-TEM), a PVP layer enveloping the surface of a Ag nanowire is observed (Fig. 1c), which is primarily responsible for the preferential formation of nanowires in the final product 34 . The average thickness of the PVP layer measured by TEM was found to be 1.2 ± 0.5 nm. Owing to the insulating nature of PVP, the structure of the junction formed between two PVP-Ag nanowires is considered to be a resistive metal-insulator-metal junction, as sketched in Figs. 1d,1e. Resistive switching is observed when an electric field arises at the junction between two PVP-Ag nanowires as formation of a conductive filament through the PVP interlayer is promoted 37,39 . The highly-resistive nature of these junctions thus governs the electrical properties of the network prior to electroformation of conductive filaments across nanowire-nanowire junctions (referred to as junctions for the rest of the text).

electrical Activation of current transport in pVp-Ag nanowire networks
Critical, or threshold, activation is a known feature of metallic nanowires, atomic switches and nanoparticle networks. It results from the formation of low resistance conduction paths between electrodes 31 . In a typical I-V experiment, the threshold voltage depends on many factors, including network density, the nature of the individual junctions and the interelectrode distance. For a density of 0.76 nanowires/μm 2 , threshold voltage for a PVP-Ag nanowire network is relatively insensitive to variations in electrode separation and lies between 1.5-2 V 40 . We have tracked the evolution of electrical current on a time series at fixed electrode separation when a constant voltage bias larger than 2 V was applied to the network. Experimental setup is schematized in Fig. 2a, left. At both ends of the network, two probes are contacted and used as electrodes for biasing the network. Figure 2b shows the time evolution of current measured over the period under which networks undergoes electrical activation. The complete activation series resembles a sigmoidal shape. In the initial moments, current is negligibly small. An exponential increase is then observed, under which the conductance of the network increases more than two orders of magnitude in tenths of seconds. Just after activation occurs, a period of logarithmic increase in conductance is followed, before finally reaching the saturation range of the measurement system.
Detailed segments of the dynamics during activation have been plotted in the insets of Fig. 2b to outline characteristic noise fluctuations prior to and after activation. As can be seen in the insets, the activation curve does not increase monotonically but rather in a step-like manner. Neither the length of the plateaus nor their relative height follows a regular pattern. Irregular plateaus can be also observed just after network transitions to a high-conductance state. We used an atomic switch-like model for the nanowire-nanowire junction 21 to simulate these features. The model accounts for the junction-to-network dynamics by combining electrical network solutions using Kirchhoff equations, with a filament growth model for the junctions which is functionally similar to the classical HP model of the memristor 4 . The filament length continuously changes between two junction states: open, in which filament length is less than the thickness of the PVP layer (high resistance), and closed, for which filament length is equal to the thickness of the PVP layer (low resistance). Junction resistance switches back and forth between these two regimes upon changes in the electrical current going through it. Details of the model can be found in the Supplementary Information. The modelled network that was used to simulate activation series has been depicted in Fig. 2a on a two-dimensional grid, with similar density as the physical system. The graph of this network is also shown on Fig. 2a, right. The nanowire network layout closely recreates the inhomogeneous distribution of nanowires on the substrate produced by our experimental deposition method (Fig. 1a). The simulated current after applying a fixed bias voltage is displayed in Fig. 2c. The shape of the simulated activation curve closely resembles the experimentally measured curve (Fig. 2b). It exhibits a characteristic sigmoidal shape, with exponential activation at low conductance, followed by a critical activation period in which the network very rapidly switches from a low-conductance(LCS) to a high-conductance state(HCS).
In the simulation, the network current spans approximately two orders of magnitude. Due to the reduced size of the simulated network, conductance is higher in the simulated activation curve, as the bias voltage is similar in both cases. Qualitative agreement is nevertheless observed by comparing the overplotted experimental and simulated activation curves in the inset of Fig. 2c. Typically, the activation process in random networks is probed by acquiring I-V curves. In an activation time series, several network parameters such as density, local homogeneity and interprobe distance can influence the time it takes to activate the network whenever voltage is larger than the threshold, as well as the number and shape of plateaus. Furthermore, local junction activation in random discrete areas of the network due to random thermal fluctuations can preclude the formation of certain current pathways at the expense of others, as shown recently 41 . Despite this, the time series profile exhibits a universal sigmoidal shape. This universality was confirmed by activation curves acquired with different Ag-PVP nanowire networks (Fig. S1).
Further insight can be gained about fundamental physical processes occurring during activation by carefully inspecting the current and voltage distribution in the network represented as a graph, at two selected times (t 1 and t 2 from Fig. 2c). Figure 2a(right), 2d and 2e plot a graph representation of the Laplacian matrix of the electrical network. In this representation, nodes (black dots) represent nanowires and edges (red lines) represent individual nanowire-nanowire junctions. A force-directed graph layout was used to better visualize the underlying topological structure of the electrical circuit formed by the nanowire network 45 . In network regions where the local nanowire density is higher, the network is more highly interconnected. As a consequence, these areas can be regarded as a highly parallelized circuit in which current is more efficiently distributed through many individual open junctions. The local voltage variations between individual junctions are much smaller in highly interconnected regions. These regions of higher connectivity are bridged by few individual junctions connected in series, in which local voltage variations are larger (as shown in the equipotential voltage map of Fig. 2d, up). By comparing Figs. 1a and 2a with the graph, it is evident that the topological structure of the network electric circuit is comprised of a combination of parallel and series pathways between both electrodes (nodes labelled '1' and '2'). The green line in the graph of Fig. 2a (right) highlights the edges along the shortest-path between both electrodes. A breadth-first search algorithm, assuming equal weight in all edges, was used to compute this path.
In the initial state of the network, with all junctions opened, the inhomogeneous density of the nanowire deposition produces an irregular voltage distribution. Figure 2d, down and 2e show the distribution of opened junctions at t 1 and t 2 . For clarity, only the edges for which the nanowire-nanowire junctions are switched on (low resistance) are drawn as red lines, leaving the other edges blank. Larger voltage drops appear in the more resistive areas of the network, as evident from Fig. 2d, which implies that these junctions are the first to switch to a closed state. Junctions closing near the electrode and in between areas of higher density will result in discrete jumps in current, as the effective resistance of the network is reduced. The exponential increase in current just before transitioning to a high-conductance regime occurs as junctions in denser areas of the network switch in an avalanche-process, bridging the low resistance pathway. Figure 2e outlines the state of the network just before activation, in which a single pathway of closed junctions forms a connected channel. Subsequent discrete current plateaus observed after activation ( Fig. 2b (inset) and 2c) arise as parallel low resistance channels are formed in the network. When all available pathways are active, the network reaches a final conductance state. Analyzing the current distribution per junction in the modelled network, we found that only 25% of the junctions carry more than 1 µA (Fig. S2).

Short-term Memory and Stochastic Breakdown
The interplay between adaptability, structural inhomogeneities and current transport across electrical nanowire networks gives rise to non-linear temporal dynamics in the activated network after the voltage is removed, as shown in Fig. 3. Here, the nanowire network was electrically stimulated by a square pulse signal with fixed amplitude (5 V) and duration (10 seconds). This setting was chosen to ensure complete activation of the network in each series. After each activation, the voltage bias was reduced to 10 mV and current was acquired during the following 100 seconds. The network's conductance response post-stimulus (Fig. 3a) reveals a complex deactivation cycle, non-monotonically decreasing and abruptly transitioning between plateaus with different conductances at irregular time intervals. This is in stark contrast to a fixed resistor network, which would maintain the same conductance, but almost immediately transition to proportionally lower currents when reducing voltage. The observed nanowire network behavior is instead consistent with a short-term or retention memory that can be attributed to the collective response of the individual junctions. Current transport is necessary to stabilize each junction 46 , but can then remain in a closed state even when the bias is reduced to a value much lower than the activation voltage (or sub-threshold level). Activation-deactivation cycles were repeatedly applied to the network and different deactivation traces were recorded. A selection of typical traces is shown in Fig. 3b, which shows that in the first 100 seconds after the pulse is applied, the network deactivates and transitions to a low conductance state. Deactivation traces regularly exhibit similar characteristic patterns as in Fig. 3a, with sharp drops in conductance followed by irregular plateaus. In many traces, the network abruptly loses its conductance after voltage removal. In other cases, HCS is retained for long periods of time, up to 60 s (bright green line), while in intermediate cases network loses the HCS in step-wise irregular plateaus at intermediate conductance states. We (2019) 9:14920 | https://doi.org/10.1038/s41598-019-51330-6 www.nature.com/scientificreports www.nature.com/scientificreports/ have analyzed the distribution of times between switching events (drops in conductance) to outline the random nature of these (Fig. S3).
This activation-deactivation measurement sequence was simulated to recreate deactivation traces. Figure 3c shows the simulated network conductance response using the same network parameters and geometry as in Fig. 2a. The shape of the simulated response (blue line) closely resembles the experimentally measured response (Fig. 3a), showing similar features: high-conductance state memory retention and abrupt conductance drops followed by irregular plateaus. In the model, resistance between junctions is controlled by an internal state variable, the length of the Ag filament. This length grows with electrical current, which promotes the development of the junction, and shrinks through an exponential decay term. This term is used to reproduce the inherent instability of atomic switches at room temperature 28,46 . The physical mechanism whereby Ag filaments decay could be related to spontaneous clustering of the Ag filament whenever voltage difference is absent, in order to minimize surface energy of the filaments 47 . Thus, we have included a dissolution threshold in the filament length. If the filament length falls below this dissolution threshold, individual junctions immediately switch back to an open state. When the voltage applied to the activated network drops below a certain level, the junctions will become unstable and susceptible to breakage or shrinking. Although junctions near the electrodes or spanning denser regions of the network show stronger resilience during voltage reduction, they eventually become unstable and break, producing the distinctive sharp transitions in conductance, as highly interconnected areas are successively disconnected within the network. As can be noted by comparing Fig. 3a,c, the agreement between the simulation and the experiment is only qualitative. The typical decay times for the experimental curves are unpredictable and about 10 times longer than the simulated one (Fig. 3b), while the model, even with this dissolution threshold approximation, is deterministic, and as such, the decay time is determined by the decay term used in the simulation, which was 1 s. To further clarify the validity of this approximation, in qualitative terms, we have included in Fig. 3c the decay curve (pink) using a continuous model for the resistance evolution when the filament length shrinks 21 . In this model, there are no dissolution threshold below which junction resistance drops due to filament breakage. The exponential decay curve represents an average decay of the ensemble of individual junctions participating in current transport, which produces a smooth reduction in the equivalent resistance of the network. Evidently, this junction resistance model is not able to capture the underlying transport hierarchy existing within the network pathways, and neglects the underlying potentiation-in-decay mechanism responsible for the dynamical behavior of the network. The resistance model with dissolution threshold (or binary) combined with the observed decay traces (Fig. 3b), on the other hand, suggests that the different deactivation responses can be attributed to two effects: the stochastic nature of the filament-breaking junctions 28,46 and the potentiation of www.nature.com/scientificreports www.nature.com/scientificreports/ the backbone current pathways when very low voltage is present. This combined effect gives rise to the creation of multiple meta-stable states of the network, appearing as irregular plateaus at different conductances with different lifetimes. Thus, this mediation of stochastic thermal junction breakdown not only probes the existence of several conduction channels actively participating in current transport, but also produces a continuous and adaptive readjustment of the preferred conductance pathways when the network loses conductance. The fact that we observe conductance plateaus signals the extraordinary fault-tolerance of the network, as the backbone of the network must carry more current to maintain the same conductance level. Figure 4 shows the nanowire network evolution upon applying consecutive I-V cycles. The first upward ramp (blue line, labelled '1') shows the network transitioning from a disconnected or low conductance regime, which produces no detectable current, up to the exponential activation of the network when the voltage is above activation threshold (cf. Fig. 2). The voltage is immediately ramped down and in the first downward cycle, the network current reduces approximately linearly. The fluctuations observed in the downward trace can indicate a low degree of connectivity in the network. On the next upward cycle (blue line, labelled '2'), sharp activation occurs again but at a lower threshold voltage. The tendency is repeated in the two following upward cycles. As the I-V waveform was 10 Hz, the time between consecutive ramps is less than the typical decay times (short-term memory) observed in Fig. 3, and as memory is retained in successive I-V cycles, upward voltage sweeps strengthen the network connectivity until the final cycle ('6'), in which both upward and downward curves show a completely stable linear profile, indicative of an ohmic system. This state is stable as long as I-V cycles are driven in the same positive voltage range and compliance current.

Transport Optimization and 1/f-Power Law Dynamics
The time series profiles in Figs 2 and 3 reveal underlying features of the switching across the network when voltage is applied continuously or in discrete pulses, while the IV cycle in Fig. 4 contains information about the evolution of the network toward a preferred or definitive high conductance state. To gain further insight into the interplay between the switching dynamics and the characteristic activation cycle, we devised a new measurement scheme to combine the time series and I-V acquisitions. In each voltage step of the I-V cycle, current is acquired for a given amount of time (20 s), the time series for each step are analyzed independently instead of averaged. The current time series acquired during two consecutive voltage steps during the upward ramp is shown in Fig. 5a. The network was in a voltage close to the critical activation. The current increases approximately linearly in this range as voltage is ramped up. The current fluctuations are in the order of tenths of μA. After acquisition, we reconstruct the I-V curve by plotting every voltage step during the upward (and downward) cycles against the span of current measured during that fixed amount of time. Thus, current time series becomes a vertical line when plotted as current versus voltage, as shown in Fig. 5b. This is in contrast to the conventional averaging method that results in a single point for every voltage (black dots in Fig. 5b). This analysis method allows us to observe not only the current vs. voltage distribution of the I-V cycle, but also the relative stability of the network www.nature.com/scientificreports www.nature.com/scientificreports/ against voltage variations (as accounted for by longer (shorter) lines when network current fluctuates more (less) for a given voltage). The resulting curve is displayed in Fig. 5c, in which the lines corresponding to the upward (downward) cycle are colored blue (orange). The activation cycle is similar to that seen in the I-V curves in Fig. 4, but with noticeable differences. First, the onset of a high conductance regime occurs at a slightly larger voltage. This implies that the activation dynamics of the network is altered when following this measurement scheme. The duration of every voltage step (20 s) is longer than the typical time to deactivate a network (Fig. 3b), and the height of each voltage step (14 mV) is small compared to the threshold voltage (~2 V). The increase in fluctuations just before network activation (Fig. 5a) suggests that there is an increased competition between thermal junction breakdown and junction formation. Additionally, the network topological layout and nanowire density inhomogeneities can also affect the activation sequence.
We computed the power spectral density (PSD) of every current time series. Similarly to performing a Short-Time Fourier Transform (STFT), time series are subdivided in fixed intervals (one second each in this case). Spectral analysis is computed on each of these separately allowing to disentangle different dynamical www.nature.com/scientificreports www.nature.com/scientificreports/ regimes occurring during the I-V cycle. Figure 5d shows a PSD computed from the time segment shown by a red box in Fig. 5a. Changes in the characteristic exponent of the resulting power-law dependence reflect dynamical changes in the network. The exponent, also termed β 48 is the slope of the line of best fit to the PSD in a log-log plot. This non-zero slope arises as a consequence of correlations in the current fluctuations of different magnitude caused by the switching dynamics 18 . Figure 5e shows how β changes during the complete I-V cycle, as a function of the average conductance during each time frame.
Two different network conductance states are apparent upon inspection: the low conductance state regime (LCS, blue dots), which spans a range 10 −10 -10 −6 S, and the high conductance state (HCS, orange dots), above 10 −5 S. The LCS and HCS regimes coincide with voltages corresponding to the upward and downward ramps, respectively. In the LCS regime, all β values are scattered between 0 and −2, and exhibit a steep linear increase up to ~10 −9 S, followed by an approximately constant but disperse distribution of values around −1.5 until the activation voltage, at which conductance is 10 −6 S. The HCS regime shows a much different distribution, with β values between −0.5 and −2.2 centered in a narrow conductance band near 2-3 × 10 −5 S. The magnitude of the power-law exponent β is regularly used in the analysis of different dynamical processes. A magnitude of the exponent close to 1 is often referred to as "1/f noise", and has been found to appear in many dynamical processes. Despite its ubiquity, the origin of such behavior is still a matter of debate 48 , but may be related to the system in a state of self-organized criticality 49 , in which the interaction between low-frequency and high-frequency processes is hierarchically organized in a scale-free, coherent manner. This self-organized state has been observed in neural dynamics 43,50 . Processes with β exponent close to 2 are characteristic of Brownian systems 51 .
When the network is fully inactive (<10 −9 S conductance), it can be considered as a system with very large equivalent resistance, as current falls in the noise level of the amplificator. The approximately flat PSD is thus a consequence of thermal (white) noise induced by the sub-nA current through the amplification system, and any physical process occurring in the network are masked by a larger noise amplitude. As voltage increases, junctions start to close, diminishing the equivalent resistance in discrete steps (Fig. 2c), and different regions of the network becoming connected. In Fig. 5e, the linear dependence in the β values at |β| > 0.5 reflects the degree to which the current fluctuations produced by discrete junction switching overcomes the amplificator's noise. Beyond the nA regime, the signal-to-noise ratio is sufficiently high that the switching dynamics of the network is the dominating factor in the entire spectra (Fig. 5d). As discussed before, in this pre-activation regime, the avalanche-like junction opening that produces sharp and strong transitions in the I-V cycle is counterbalanced by thermal desorption of single junctions due to large acquisition times for every point, and β values are distributed around -1.5. Just before activation (rightmost 10 blue circles in Fig. 5e), the much larger and irregular current fluctuations (rightmost blue line in Fig. 5c) are distinguished by β values close to −2, after which a sharp transition to the HCS occurs (red points in Fig. 5e). Figure 5f shows the corresponding probability distribution of |β|values when the network is in the HCS. The β distribution is centered closed to -1(average value of −1.2.) showing an asymmetric tail at |β| > 1.5. As is evident from comparison of the corresponding data (red dots) in the downward ramp of Figs. 5c,e, the network HCS is very stable against current fluctuations. This resilience against even the strongest current fluctuations (longest vertical lines in the downward ramp data in Fig. 5c) can be attributed to its ability to adapt and redistribute signal power to lower-frequency collective modes. We ascribe the tail of the β distribution observed in Fig. 5f to these unstable states, manifested as sharp, irregular and random fluctuations. In this particular example, the network is relatively stable during the whole downward cycle, and we estimate from the cumulative distribution in Fig. 5f that the probability to enter this unstable regime is 3% (shaded box in Fig. 5f). By repeating the experiment for different networks, the transition to a definitive or stable high-conductance state is also observed, but this state is sometimes more unstable so the whole distribution is shifted towards larger β values, as can be seen in the Supplementary Information (Fig. S4).

Discussion
It has been shown in previous sections that PVP-Ag nanowire networks evolve towards a definitive state of conductance after becoming fully activated, either by means of repetitive I-V cycles or applying super-threshold voltage over longer times. In the activated network, one or many pathways transport electrical currents. The modelling results (Fig. 2a,d,e) suggest that the pathway that opens first is the shortest-distance topological path between the electrodes. However, the model assumes a homogeneous initial distribution of junction resistance, which may not be the case in the real system, which has a larger interelectrode distance, as well as non-ideal contact geometry between nanowires and between nanowires and electrodes. Additionally, the self-assembled network may exhibit larger homogeneity variations that are not accounted for in the model. This implies that filament formation at the level of individual junctions might not always lead to an efficient activation sequence. For example if one junction in a less optimized pathway change its resistance state first, it may promote the activation of junctions in its immediate vicinity, thus giving rise to a winner-takes-all phenomenon, as recently observed by Manning and colleagues 41 in a smaller PVP-Ag network. In any case, once a pathway is active, the network will keep further opening parallel pathways in order to optimize the conductance, or else, minimize power consumption, in order to comply with the laws of conservation of electrical current in a resistor network Previous studies 36 have shown the stability of active networks to last extremely long times as long as the energy pumped from a voltage or current source is high enough. Thus, once in the ohmic regime, a network can remain in this state unless, by means of a higher voltage or current injection, other physical phenomena such as joule heating or electro-migration breaks the individual junctions, thereby resetting them 36,52 . In our study, we have inspected the dynamics of the active network not only near the threshold voltage, but also in the sub-threshold regime. In this regime, the stability of the networks is still very high, but as the current through individual junctions decreases, the balance between junction formation by means of current injection and junction dissolution triggered by stochastic fluctuations is reduced. Decreasing the voltage in a linear sweep while acquiring time www.nature.com/scientificreports www.nature.com/scientificreports/ series facilitates the inspection of the dynamical events, signaling network instability by single junction breakdowns. This is assessed by analysis of the PSD power-law exponent (β) distributions in Figs. 5d,e. The β distribution in Fig. 5e shows that most of the magnitudes are centered around 1, with a small tail of values rounding up to a magnitude of 2. Typically, a 1/f power-law spectrum is indicative of a system in which scale-free dynamics is at play. To better understand how the network changes its dynamical behavior, we dissected the downward voltage ramp of the I-V cycle of Fig. 5c and examined the conductance time series at a voltage in which this dynamical balance was evident. Figure 6a shows that the network maintains a very stable (small fluctuations) conductance near 2.6 × 10 −5 S up to approximately 15 s, after which it suddenly becomes unstable, exhibiting sharp and irregular conductance drops interspersed with short plateaus of unequal duration and magnitude, as well as large spikes. This unstable regime last for 3 s, after which the network automatically recovers the previous stable state, only this time with a slightly higher conductance. The level of fluctuations in this new stable regime appears to be slightly lower. Figure 6a also shows β values calculated at every one second interval in the conductance time series. This analysis reveals that when the network is in a dynamical regime that is stable against fluctuations, β values are closer to 1, whereas they approach 2 when the network succumbs to instability. Figure 6b shows PSDs calculated for a stable and unstable region in the conductance time series in Fig. 6a, indicated by the black triangle and pink diamond, respectively. The PSDs clearly distinguish the different dynamical regimes: the PSD is steeper in the unstable regime due to the stronger low-frequency fluctuations. Upon further inspection of Fig. 6a, with reference to the full distribution of sliced time series in Fig. 5c, it is evident that the onset of the unstable regime is not triggered by the negative voltage sweep, nor it is related to changes in the network appearing immediately after one voltage step is performed in the downward ramp, but instead appears randomly and unpredictably, and with an irregular duration. The unstable regime ends when the network recovers a stable conductance state, as signaled by reduced noise fluctuations. Therefore, the network not only adapt its conductive pathways to the nanowire junction topology whenever energy is supplied either from current or voltage sources, but can also adapt dynamically to random changes in the network, reconfiguring the underlying connectivity to find a new optimal pathway.
Few nanowire network models have addressed stochasticity in the sub-threshold regime 21 . Our modelling of the activation and deactivation cycles (cf. Figs 2 and 3) provides new insights into the process leading to these changes in network. When the applied voltage is much lower than the threshold voltage, the transport hierarchy is revealed through the different irregular decay curves (Fig. 3b), attributed to stochastic junction breakdown. In reducing the voltage from threshold to sub-threshold level in small steps and long waiting times, we were able to track the transition from the stable state down to the regime in which random fluctuations have a larger influence on the dynamics. Whenever β is close to 1, individual junctions within clusters of high network connectivity are continuously switching off and on, with three different dynamical process (potentiation, inhibition and random breakdown) interacting and evolving in a self-organized way such that the connectivity backbone of the network is not compromised. Therefore, a junction that breaks in a cluster of high connectivity is propagated as a small change in conductance, and, depending on the local topology of these clusters, induces the breakdown of neighboring junctions, as the current in them is instantaneously reduced. But at the same time, when these secondary junctions are broken, the backbone of the network needs to carry more current to maintain conductance, so if a region is locally connectivity-depleted, another region will have to sustain a larger current, and thus, will potentiate the switching-on of junctions over that area. However, the more junctions are switched on, the higher the probability of breaking junctions in these newly connected areas. In this way, self-organized feedback produces the scale-free power-law distribution, as fast individual junction switching leads to changes in conductance for clusters of junctions over a slower time scale, and changes in cluster conductance produce changes in network pathway connectivity on an even slower time scale. This hierarchy is responsible for the characteristic 1/f spectrum in the PSD. From this interpretation, we can infer that the steepening of the power-law spectrum associated with abrupt changes in conductance (cf. Fig. 6a) occurs when a critical junction (one bridging areas of larger connectivity) breaks down, and splits a network pathway, reducing the total network conductance (by as much as half in the ideal case of two equivalent parallel pathways). Then, self-organized switching is disrupted by larger junction connectivity-depletion, and the connectivity of the network is spontaneously remapped until the broken pathway, or another parallel pathway, is switched back on. This remapping (or reconfiguring) process produces more low-frequency power in the PSD, as clusters of junctions or secondary pathways collectively redistribute connections, which occurs on longer time scales (lower frequencies) than the higher frequency switching fluctuations of individual junctions. When the network transitions back on to a highly conductive state, the conductance is slightly different, implying this complex dynamical process produces long-lasting changes in the backbone connectivity of the network. We can draw a parallel with long-term retention memory in the brain, which requires not only synaptic potentiation in neural connectivity but also persistent structural reconfigurations in connectivity, referred to as plasticity 53 .
Further insight into the relation between dynamics and connectivity can be gleaned by encompassing the neuromorphic nanowire network within a broader class of systems, optimal transport networks. In general, these are networks that can be represented by a weighted graph (directed or undirected), with edges representing the conductance (or resistance) between the different parts of the network, and whose solution can be obtained by solving a generalized version of Kirchhoff 's equations, that is, flow and mass conservation. Examples of such systems vary from artificial systems, such as power or water supply networks, to natural systems, e.g. vascular networks in plants or animals 54 . The optimality condition, though not always solvable, is computed by considering that transport is optimized whenever the connectivity map (weight between edges) minimizes power consumption. For example, for a fixed resistor network, current is hierarchically and homogeneously distributed in proportion to the equivalent resistance of the different parallel plus series pathways. In contrast, in vascular systems, which adapt flow to pressure differences 54 , dissipation is minimized by concentrating flow on the edges with larger conductance, thus transport is regulated by the formation of tree-like structures 55,56 . Neuromorphic networks clearly show this adaptive behavior, with the system evolving to create few, large connectivity pathways rather than many different pathways, as our activation curve and other recent research shows 41 . However, once the network is connected, parallel pathways, most likely in the form of closed loops or circuits that assist or reinforce the main pathway also appear. This further reduces global power dissipation. Interestingly, research on optimal transport structures suggests that the combination of a fixed unique backbone for transport with different recursively nested loops or cycles is the preferred connectivity map for adaptive networks, and appears once variations in the transport dynamics are introduced in the form of defects or fluctuating currents in the junctions, sources or sinks [55][56][57][58] .
Finally, it is worth considering whether the network dynamics observed in our neuromorphic nanowire system are characteristic of biological neural systems. In particular, we always observe a continuously fluctuating state, as if the system needs to be constantly oscillating, searching the phase space for a solution that guarantees long-term stability. It has been suggested that this edge-of-chaos state could be not only an optimal state of the brain, but a necessary condition for consciousness to exist 43,59 . However, although neural dynamics research has tended to focus on wave-like collective oscillations (e.g. in EEG signal acquisitions), other time series analyses continue to be developed 50,60,61 . Power-law PSD spectra, such as that observed in our nanowire networks, are indeed also observed in neural dynamics (e.g. neuronal avalanche size distribution), although whether these spectra can be interpreted as evidence for the brain operating at or near criticality remains a contentious issue 62 . Changes in power-law slopes, as observed in our nanowire network system, have also been observed in brain measurement data, where a subject performed simple tasks 60 . In our neuromorphic nanowire network, large changes in the power-law dynamics arise as random fluctuations break critical junctions. Introducing a control mechanism to disrupt selected regions of the network could force this remapping mechanism, that ultimately will allow us to harness the natural nonlinear feedback dynamics, thus navigating the large solution space of the network to solving different optimization or recognition tasks.

conclusions
We have analyzed the emergent dynamical properties of a neuromorphic system formed by polymer coated silver nanowires that self-assemble into a network with resistive switch electrical junctions. Our electrical measurements, performed with a large interprobe distance across the network, reproduced two important current-voltage properties observed in similar neuromorphic systems: critical activation at a threshold voltage; and evolution towards an optimal conductance state under successive cycles. We devised a new measurement scheme that unveiled even richer and more complex network dynamics arising from multi-scale interactions of the switch junctions, including: (i) a collective memory response in the sub-threshold voltage regime; (ii) distinct low and high network conductance states distinguished by different power-law fluctuation scalings; and (iii) reconfiguration dynamics, similar to synaptic plasticity in the brain, enabling resilience and adaptation. In particular, we found that the emergent temporal correlations arising in the active network under changing voltages exhibit a 1/f power-law spectrum, which is a consequence of fundamental changes in the connectivity map of the network and is primarily mediated by stochastic fluctuations at the individual junction scale that lead to hierarchical collective dynamics, from the scale of junction clusters to the whole-of-network scale. We envisage that by accessing the network's states through its spectral signatures, we may not only gain a better understanding of its near-criticality dynamics, but may also discover new strategies for controlling training and learning for reservoir computing applications.