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

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.

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

$$\begin{equation}{\Delta}{d}_{i}=\begin{cases}{r}_{\mathrm{d}}({E}_{i}-{E}_{\mathrm{T}}),\quad \hfill & \mathrm{i}\mathrm{f}\enspace {E}_{i}\geqslant {E}_{\mathrm{T}}\hfill \\ 0,\quad \hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{cases}\end{equation} \tag{ 1 }$$

Once a filament has formed (with an initial width w₀), its destruction can be explained by electromigration effects [35] where the current (I_j ) through the filament decreases its width (w_j ) according to

$$\begin{equation}{\Delta}{w}_{j}=\begin{cases}{r}_{\mathrm{w}}({I}_{j}-{I}_{\mathrm{T}}),\quad \hfill & \mathrm{i}\mathrm{f}\enspace {I}_{j}\geqslant {I}_{\mathrm{T}}\hfill \\ 0,\quad \hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{cases}\end{equation} \tag{ 2 }$$

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 Ω⁻¹. 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.

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 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 ($\bar{G}$, where $\bar{\,}$ 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.

Zoom In Zoom Out Reset image size

Figure 4(a) shows the properties of square 2ESs as a function of L (left column, W = L), while figure 4(b) shows results for rectangular 2ESs with different aspect ratios (right column, W = 200). Figures 4(a) and (b) shows $\langle \bar{G}\rangle$ (top row), where ⟨⟩ indicates an average over 24 network realizations. Also shown are the average number of active sites (⟨N_A⟩, middle row) [at which a switching event occurs at least once] and the average density of active sites (⟨ρ_A⟩ = ⟨N_A⟩/(LW), bottom row). Results are shown for different levels of voltage stimulus, in order to demonstrate the effect of the voltage distribution throughout the system.

Zoom In Zoom Out Reset image size

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.

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.

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).

Zoom In Zoom Out Reset image size

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.

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: $\bar{G}$ = 0.193; 4ES: $\bar{G}$ = 0.175 Ω⁻¹; AES: $\bar{G}$ = 0.206 Ω⁻¹). 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)).

Figures 7(a)–(c) shows the distributions of inter-event intervals (IEIs) [times between consecutive switching 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.

Zoom In Zoom Out Reset image size

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

$$\begin{equation}P(S)\sim {S}^{-\tau }\end{equation} \tag{ 3 }$$

and

$$\begin{equation}P(T)\sim {T}^{-\alpha },\end{equation} \tag{ 4 }$$

and the mean sizes of the avalanches are related to their durations according to

$$\begin{equation}\langle S\rangle (T)\sim {T}^{1/\sigma \nu z},\end{equation} \tag{ 5 }$$

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]

$$\begin{equation}1/\sigma \nu z=(\alpha -1)/(\tau -1)\end{equation} \tag{ 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].

Figures 8(a) and (b) show that the sizes and durations of the avalanches are distributed as power laws for each of the 2ES, 4ES, and the AES (red, green and blue respectively). Figure 8(c) is a plot of the mean avalanche size given duration $(\left\langle S\right\rangle (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].

Zoom In Zoom Out Reset image size

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.

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.

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.

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.