Self-organized nanoscale networks: are neuromorphic properties conserved in realistic device geometries?

Self-organised nanoscale networks are currently under investigation because of their potential to be used as novel neuromorphic computing systems. In these systems, electrical input and output signals will necessarily couple to the recurrent electrical signals within the network that provide brain-like functionality. This raises important questions as to whether practical electrode configurations and network geometries might influence the brain-like dynamics. We use the concept of criticality (which is itself a key charactistic of brain-like processing) to quantify the neuromorphic potential of the devices, and find that in most cases criticality, and therefore optimal information processing capability, is maintained. In particular we find that devices with multiple electrodes remain critical despite the concentration of current near the electrodes. We find that broad network activity is maintained because current still flows through the entire network. We also develop a formalism to allow a detailed analysis of the number of dominant paths through the network. For rectangular systems we show that the number of pathways decreases as the system size increases, which consequently causes a reduction in network activity.


Introduction
Complex self-assembled networks of nanomaterials have recently emerged as promising systems for neuromorphic computing [1]. The general concept is that since these networks exhibit structures or behaviours that are similar to those of the brain, they may enable novel methods of brain-like information processing, with potentially attractive features such as ultra-low power consumption. Many different types of neuromorphic networks have been investigated, taking advantage of building blocks that include a variety of types of nanowires and nanoparticles [2][3][4][5][6][7][8][9][10][11]. Recent advances include demonstrations of criticality [4,12] and edge of chaos learning [13]. Furthermore, various computational tasks were successfully demonstrated using approaches [14][15][16] based on reservoir computing (RC) [17], where the non-linear dynamics and memory capacity of the network are exploited for information processing.
One of the attractions of self-assembled networks of nanomaterials is that they allow all-electrical implementations of RC, which are promising for real-world applications [18]. In these implementations electrical contacts (electrodes) are required in order for signals to be input to and read out from the network, and this leads to an important and thus far overlooked issue: different electrode configurations can modify the network dynamics. Figure 1 provides a simple example in which the choice of electrode configuration (red: inputs; pink: outputs) changes the current pathways through a network. Since the components in neuromorphic networks are memristors or switches, and the behaviour of each component depends on the local currents and voltages, differences in the current pathways can lead to differences in the network dynamics. Consequently, the Left: two-electrode system (2ES); middle: four-electrode system (4ES); right: alternating-electrode system (AES). Colours are the same as those used below in simulations of the current flow in the actual percolating network of nanoparticles. Red and pink squares denote input (applied voltage) and output (current to ground) electrodes respectively. High, medium and low currents in each path are indicated schematically with red, purple/dark blue, light blue bonds respectively. electrode configuration has the potential to impact RC performance. In a worst-case scenario, the electrode configuration could cause only a small part of the network to be active (e.g. due to concentration of the current through a small number of paths) potentially leading to poor RC performance. Note that this issue does not arise in RC implementations such as echo-state networks [19] or optical systems [20,21], as in those cases the readout and the network dynamics are decoupled.
Here we investigate the effects of the electrode configuration on the dynamics of networks of nanoscale electronic components, using detailed and realistic simulations. We consider the impact of the configuration of inputs and outputs as well as the size and aspect ratio of the network. We develop a method to quantify the number and length of pathways through the complex network and show that for a constant voltage, the number of pathways decreases with increasing system size, reducing network activity. Despite this, the network dynamics are usually not strongly affected by the electrode configuration and network geometry, and critical dynamics are observed [4]. Criticality is important because it is associated with optimal information processing capability [22], and is a key feature of the biological brain [23][24][25]. Hence we quantify the neuromorphic potential of the networks using the strict mathematical definition of criticality.
This work focuses on practical, experimentally realizable, configurations of outputs, that are suitable for implementation of RC algorithms in real-world devices. In contrast, in the RC literature it is usual to assume that output signals can be extracted from all nodes in the network [9,18,19,26]. While this is possible in simulations, access to the internal nodes is not practically feasible in a physical device.
We expect that the effects of electrode configuration and network geometry will be similar for all types of nanoscale electronic networks but focus here on percolating networks of nanoparticles (PNNs) [4,27], which are illustrated schematically in figure 2. PNNs are particularly attractive for neuromorphic computing because straightforward experimental methods exist for preparing them at a critical point, using simple costeffective deposition techniques [4,27]. The deposition of nanoparticles is stopped close to the percolation threshold (onset of conduction across the network), which ensures that the networks contain many tunneling gaps (i.e. voids between well-connected groups of nanoparticles, figure 2(b)). Tunnel gaps emulate some of the integrate-and-fire characteristics of biological neurons via atomic-scale switching processes [28]. Formation and destruction of atomic-scale filaments occurs in response to applied input signals [27], as illustrated in the inset of figure 2(b). When one switch changes state, consequent changes in the distribution of voltages/currents across the network lead to other switching events, and hence to complex patterns of correlated electrical signals that are distributed across the network [5] (note that the memristive properties of individual tunnel junctions [12] are difficult to probe directly due to this collective behaviour). The resulting switching dynamics are consistent with avalanches observed in cortical tissue [29] and meet strict statistical criteria for criticality [4].

Simulations of percolation with tunneling and switching
It is well-established that the experimental nanoparticle networks can be well-described by simulational models of 2D continuum percolation [27,28,[30][31][32]. [For completeness we include a brief comparison with other types of percolating system in supplementary section I (https://stacks.iop.org/NCE/2/024009/mmedia).] The deposited particles are represented by overlapping disks, which form groups as they land on each other, as in the experiments. Below the percolation threshold (surface coverage p < p c ∼ 0.68), the overlapping disks do not form a continuous path and so conduction is possible only due to tunneling through the gaps between groups of particles/disks (figure 2). The conductance of each gap is G i = A exp(−δL i ), where A and δ are constants and L i is the gap size in units of the particle diameter (for convenience the particle diameter is set to 1 and A = 1Ω −1 ; it has been found that δ = 100 gives good agreement with experiments) [27,32].

Key parameters
We present results for p = 0.65 as this coverage is close enough to p c to allow observation of critical dynamics but far enough away to ensure that unrealistic network configurations (e.g. those with very large groups and hence a small numbers of gaps between the electrodes) can be avoided. Note that we have simulated systems with a range of surface coverages and find similar results throughout the range 0.64 < p < p c [5]. We focus on square systems of size 200 × 200 particle diameters, as this is large enough to produce critical dynamics without being computationally exhaustive. Supplementary section II presents similar results for 400 × 400 systems with a greater number of electrodes.

Model of switching behaviour
Upon applying an external voltage stimulus, atomic scale filaments can form within the tunnel gaps. Atomic filament formation is described by a deterministic model [28] in which the formation is driven by the electric field (E i ) in the ith gap via the processes of electric field induced evaporation and electric field induced surface diffusion [33,34]. The size of each tunnel gap (d i ) changes according to Once a filament has formed (with an initial width w 0 ), its destruction can be explained by electromigration effects [35] where the current (I j ) through the filament decreases its width (w j ) according to E T (I T ) represents the electric field (current) threshold, above which d i increases (w i decreases) at a rate of r d (r w ). E T = 10 V (per unit length, measured in particle diameters) and I T = 0.01 A are chosen to best reproduce experimental results [27,28]. After the formation of a filament, the conductance is set to 10 Ω −1 . Note that the precise conductance values could be scaled to match the experimental values more closely but do not affect the overall results and so are chosen to maintain consistency with previous work [32]. This model is described in detail in [28] and exhibits criticality and avalanches of signals that are similar to both the experimental data and to neuronal signals in the cortex [4,28].
We emphasise that this model of the switching behaviour is similar to general models of memristive switching [36]. In fact the equivalence between the models can be made exact for appropriate choices of the resistances in the 'on' and 'off' states of each switch. Hence, taken together, our models of the network structure and switching dynamics capture the essential features of a wide range of self-organised nanoscale networks.

Two-electrode systems-basic properties
Experimental studies of PNNs [4,5] have focused on simple two-electrode systems (2ESs), with the input and output electrodes on the left and right hand edges, respectively (figure 2(a)). Therefore we begin by The input voltage is applied to the electrode on the left side (3V) while the electrode on the right side is grounded. Inset 1 (left): initial increase in conductance upon applying a voltage. As filaments are formed in multiple gaps throughout the network, G increases. Inset 2 (right): patterns of switching behaviour. After ∼10 4 time steps the conductance fluctuates around an average value (G, red horizontal line) as the system is in a meta-stable state in which a dynamical equilibrium is maintained (where filaments are continuously formed and broken). (b) Illustration of the groups of particles, which are separated from each other by tunnel gaps. The groups forming the input electrode (left) and output electrode (right) are shown in black with the other colours representing different particle groups. (c) Active sites (tunnel gaps that switched at least once) over 5 × 10 5 time steps are shown in red while the shades of blue represent different particle groups. (d) The distribution of potentials at each group. Input and output electrodes are shown in black with the colours representing the potential gradient from high (red) to low (blue). (e) Current map. Each filled circle/black dot represents the centre of a group of particles with lines denoting a connection between groups via a tunnel gap. Red and pink circles represent the groups connected to the input and output electrodes, respectively. The line colours indicate the magnitude of the current through the network, red (high) to blue (low). Note that the grey lines represent tunnel gaps carrying no current. (d)-(e) Are representative snapshots of the network at time step 5 × 10 5 .
illustrating the basic properties of those devices. Figure 3(a) shows the conductance as a function of time for a long simulation with a DC applied voltage. The increases (decreases) in conductance each correspond to switching events in which a filament forms (or ruptures). These results are very similar [28] to those observed experimentally. The left inset shows the initial increase in conductance before the system settles into a dynamical equilibrium where formation and destruction of filaments are persistent and conductance fluctuates about some constant mean value (G, where indicates an average over time). The right inset shows the patterns of switching behaviour, which include avalanches of activity as previously discussed in [28]. Figure 3(b) shows a map of the percolating network, with the different groups of particles shown using different colours. The active switching sites (tunnel gaps) are shown in red in figure 3(c). Figures 3(d) and (e) show snapshots of the voltage gradient and current paths though the network respectively. Note that the 2ES is used as an example here and systems with different sizes, aspect ratios, and electrode configurations are discussed below.

Results
As noted above, previous experimental studies of PNNs focus on simple 2ESs (figure 2(a)) but devices with multiple electrodes (figure 2(b)) are required for applications like RC. In these devices, the positions of the electrodes around its periphery define the size and shape of the active part of the network. We first simulate the effect of the size of the system on the switching activity and conductance (section 3.1) and the number and length of the current pathways through the network (section 3.2). We then consider the effect of shape in 2ESs with different aspect ratios, L/W. Note that here we use 'length' (L) for the size of the system in the horizontal direction (i.e. electrode separation) and 'width' (W) for the size of the system in the vertical direction. We specifically consider the effect of increasing L when W is fixed (as in figure 2(b) where W L for opposite electrodes; in our experimental multicontact devices W L ∼ 1 14 [37]). Finally, we investigate the effect of electrode configuration in systems with multiple electrodes (section 4). In square systems (figure 4(a)), for a 3 V stimulus, the conductance decreases with increasing system size. This is because at large L current pathways between electrodes contain more tunnel gaps in series and so the average voltage across each gap is lower. Consequently, the decrease in ρ A shows that there is reduced switching activity due to a lower proportion of tunneling gaps maintaining a significant electric field across them (i.e. E > E T , above the threshold for switching in equation (1)). For a 50 V stimulus, the conductance stays nearly constant with increasing system size because at this stimulus level the proportion of tunnel gaps that exhibit switching activity remains constant. Hence, ρ A is nearly constant and N A increases in accordance with the increase in area.

Size and aspect ratio-effect on switching activity and conductance
In rectangular systems ( figure 4(b)), the conductance decreases with increasing L, for all voltage stimuli. This is because, as with the square systems, the number of tunnel gaps in series increases with L. The average voltage across each gap then decreases, resulting in a dramatic reduction in switching activity: N A decreases to zero for L > 500 with a 2 V voltage stimulus. For a large stimulus (50 V), the average voltage across each gap increases causing N A to increase for all L.

Size and aspect ratio-effect on current pathways
The size and aspect ratio of PNNs also affect the number (N P ) and length (L P ) of the dominant current pathways that span the input and output electrodes. Figure 5(a) shows N P and L P (determined via the algorithm described in supplementary section III) as a function of L for square systems (left column). N P increases with L initially (from L = 100 to 200) but decreases for L > 200. This is because for L 200 the 3 V stimulus is sufficient to maintain E > E T across many tunneling gaps. However, as L increases, L P increases, as shown in figure 5(c). This is because the increasing L increases the number of tunnel gaps required for a path to span the network, reducing the average voltage across each gap. As already mentioned in section 3.1, some gap voltages consequently fall below the threshold for switching so that they remain open, decreasing N P . This results in relatively low currents, even along the highest conductance paths. The increase in L P and the overall decrease in N P both act to decrease the mean conductance, as shown in the upper panel of figure 4(a).
The number and length of the pathways through rectangular systems evolve in a similar way, as shown in figures 5(b) and (d). Since W is fixed, as L is increased the number of parallel paths decreases even more quickly than in square systems; L P increases roughly in proportion to L and N P decreases to zero for large L.

Systems with multiple electrodes
For RC we require multi-electrode systems (MESs) with multiple inputs and outputs [14,15,26,38], as shown in figure 2(b). Data for MESs with many electrodes is presented in supplementary section II, but in this section we illustrate the key points by considering simpler systems with just four electrodes.
We compare three cases for the same 200 × 200 network: a 2ES, a four-electrode system (4ES) and an alternating-electrode system (AES). These configurations are the same as those illustrated schematically in figure 1 but the network being modelled is a PNN, rather than a resistor array. As described above, the standard 2ES (see current map in figure 6(g)) has input and output electrodes along the left and right edges of the system, respectively. For the 4ES (figure 6(h)), input (output) electrodes are connected to two different groups on the left (right) edge. For the AES (figure 6(i)), the input electrodes are connected to groups on the left and right edges, while the output electrodes are connected to groups on the top and bottom edges. Note that in simulations the electrode positions for the 4ES and AES are selected to ensure they are well-connected to the network, i.e. the electrode has a high conductance connection to a group of particles. In experimental devices the electrodes are much larger (as is the system size: L 10 4 particle diameters), which ensures that they are well-connected to at least one group. Also note that the configuration of the electrodes in experiments can be varied in many ways, but these arrangements are a reasonable starting point for investigation of the impact of electrode placement on the network dynamics.
In section 4.1 we illustrate the effects of the electrode configuration on the basic device properties (current flow and switching dynamics). In section 4.2 we investigate whether the critical avalanche dynamics observed in 2ES [4] are still observed in 4ES and AES geometries. Figure 6 shows the conductance variation with time (top panels) for a 2ES (figure 6(a)) compared to a 4ES ( figure 6(b)) and an AES (figure 6(c)). The switching dynamics shown are qualitatively similar for all three systems, with similar average conductances (2ES: G = 0.193; 4ES: G = 0.175 Ω −1 ; AES: G = 0.206 Ω −1 ). The maps in figures 6(d)-(f) show that the switching activity is qualitatively similar for all three systems but is concentrated near the electrodes (the groups of particles coloured in red) in the AES because the tunnel gaps on shorter paths are more likely to exhibit switching behaviour. This is because there are fewer tunnel gaps in series, reducing the path resistance and increasing local potentials. Supplementary section IV shows an example of a 400 × 400 AES which more clearly demonstrates the concentration of the current. In both cases, despite the concentration of current near the electrodes, currents flow throughout the entire network (figures 6(g)-(i)). L P . The current threshold is 1% (see supplementary section III for further details). The averages N P and L P were calculated for 24 realizations of the network. Each panel shows the mean value (green) as well as the maximum (blue) and minimum (brown) values. Note that these values are calculated from one representative instant in time after the system has reached dynamical equilibrium (as in figure 3) but the trends remain the same across multiple instants in time.

Figures 7(a)-(c) shows the distributions of inter-event intervals (IEIs) [times between consecutive switch-
ing events] are power law distributed with similar slopes for all of the 2ES, 4ES, and AES configurations. This is indicative of scale-free temporal dynamics [5]. Additionally, the autocorrelation functions (ACFs, figures 7(d)-(f)) for all three systems show strong correlations between events. Hence the 2ES, 4ES, and AES all exhibit complex, correlated switching dynamics with power-law distributed IEIs. These are necessary but not sufficient criteria for criticality and so we explore criticality more rigorously in section 4.2.

Critical avalanches
Critical avalanche dynamics have been observed both in cultured cortical tissue [39] and in experimental 2ES PNNs [4]. As critical dynamics are thought to maximize several key information processing metrics [22], it is crucial to determine whether criticality in PNNs is conserved under different electrode geometries. In accordance with previous analyses of critical systems [40], the size (S) and duration (T) of each avalanche of signals is defined by counting the total number of events in the avalanche and the number of time bins over which the avalanche propagates. For critical systems, S and T have power-law distributions and and the mean sizes of the avalanches are related to their durations according to where 1/σνz is a key characteristic exponent. The critical exponents are discussed in detail in [41] and the analysis methods in [4].
The key test for criticality is that the exponent obtained from equation (5) should be consistent with two other estimates of 1/σνz. The first is obtained from the 'crackling relationship' [42] 1/σνz = (α − 1)/(τ − 1) (6) and the second from avalanche shape 'collapse' i.e. scaled avalanche profiles should collapse onto a single curve, yielding an independent estimate of 1/σνz [41].  Figure 8(c) is a plot of the mean avalanche size given duration ( S (T)) for each system and figures 8(d)-(f) demonstrates shape collapse. Table 1 shows the critical exponents obtained from these distributions (following the detailed methodology of [4,39]). Importantly, table 1 shows that in each case the three estimates for the critical parameter 1/σνz are in agreement within statistical uncertainty, thereby demonstrating that rigorous criteria for criticality are satisfied for all three systems [4,39,41].
It is interesting to note that some previous studies focused on the response of percolating networks to a stepwise increase (decrease) in the input voltage [43,44]. In that case the networks exhibit rather different avalanches, in which the active switches all turn on (off) in succession, i.e. the conductance transitions occur in a single direction, causing the mean conductance to increase (decrease). The critical size exponent is then τ = 1.5 [44], which is much smaller than that shown in table 1. This difference arises from the fact that the avalanches considered here occur in a balanced state [28] (for a constant DC applied voltage) where there are both increases and decreases in conductance around some mean value, as shown in figures 3 and 6.

Discussion
In the preceding sections criticality was observed for a range of network geometries and electrode configurations. In this section we discuss cases in which critical dynamics are not observed and practical factors which can make the demonstration of criticality difficult in large systems. Finally, we comment on possible choices of the type of electrode.

Non-critical cases
Simulations in which criticality is not observed are reported in detail in supplementary section V. In these cases regular spiking behaviour occurs, where the same switch turns on and off repeatedly. [Similar spiking behaviour was also reported in recent experiments [37].] The relevant quantities (S and T) then no longer follow power law distributions and the criteria for criticality are not met. Spiking dominates the pattern of switching activity in these electrode configurations because a small number of switching sites are dominant, as shown in figure S8b. This problem can be exacerbated in larger systems (see next section) where there are a large number of candidate switching sites, each of which may exhibit spiking behaviour at different times. While in principle the regular spiking can be removed from the data by filtering or thresholding, this is not straightforward in networks with short path lengths where large amplitude spiking behaviour can be observed, or when different sites cause spiking at different times.

Practical factors in larger systems
Supplementary section II presents simulated data for 400 × 400 square systems. While these larger systems are critical, there are two contraposed factors that make performing the criticality analysis challenging. Firstly, larger systems have more switching sites, which all have the potential to contribute to the observed switching activity, and at high applied voltages could lead to very high switching rates. The criticality analysis [4,39], however, requires identification of avalanches of switching activity that are clearly separated from each other. This means that the applied voltage cannot be too high. Secondly, as discussed in sections 3.1 and 3.2, since the applied voltage is distributed across more tunnel gaps (in series with each other), a higher applied voltage is required to achieve electric fields larger than the threshold E T (see equation (1)) in a significant number of tunnel gaps across the network.
These contradictory requirements lead to a narrow range of voltages in which the criticality analysis can be performed. At lower voltages too few events are observed to perform the analysis, and at higher voltages there are too many. This constraint is important in the analysis of simulations, where a relatively high event rate is desirable to reduce computation time. Similar issues can be avoided in experiments performed at high sampling speeds-e.g. MHz sampling rates [37] ensure that it is straight forward to demonstrate criticality.  Table 1. Demonstration of criticality for 200 × 200 systems with different device geometries. The critical exponent 1/σνz is obtained from the crackling relationship (α − 1)/(τ − 1), mean avalanche size given duration S (T ), and avalanche shape collapse. The agreement of these three independent estimates of 1/σνz is a rigorous requirement for criticality. For more criticality analysis details see [4,39,41].

Types of electrode
Finally we comment on the differences in the types of electrodes that must be used for measurement of nanoscale electronic networks and in neuroscience. An important difference is that the propagation of electrical signals through nanoscale networks can be modified by the number and placement of electrodes whereas in neuroscience it is assumed that the electrodes are passive measuring devices and have no impact on the behaviour of a network of neurons. If high impedance voltage probes were connected to nanoscale electronic networks, i.e. they recorded voltages and not currents, they would not drain current and would have less influence on the network dynamics. They would therefore perform a similar role to those used in neuroscience. However, for all electronic systems, there must be at least one electrode through which current flows to ground and the placement of the grounded electrode or electrodes will always modify the current flow, and thereby influence the dynamics, to some extent. The remarkable feature of these results is that, despite differences in the local currents (see figure 6), the critical network dynamics are generally not affected by the electrode placement.

Conclusions
We have demonstrated the impact of device geometry and electrode configuration on the correlated switching behaviour of nanoscale electronic networks through the example of PNNs. In particular, we have shown that increasing the system size and aspect ratio results in a distribution of the applied voltage across a greater number of junctions, reducing the number of 'on' switches and decreasing the system conductance. The finite width of systems with larger aspect ratios restricts the number of conducting pathways, exacerbating this effect. We have also compared a two-electrode configuration to a multi-electrode configuration and found that the switching dynamics and correlations are qualitatively similar, and that despite the concentration of currents near the electrodes, current still flows through the entire system. We have demonstrated that the addition of multiple electrodes and the increase of system size (which are essential for RC applications) do not alter the critical dynamics previously demonstrated in two-electrode systems. The multi-electrode configurations investigated here provide an intuitive basis for a RC [17] approach where the electrodes act as the required input and output nodes. Importantly, since criticality is observed in all geometries investigated, the computational benefits of criticality are retained which is necessary for efficient RC [22][23][24][25]45]. These results therefore support further experimental efforts directed at the physical implementation of RC approaches using nanoelectronic systems.