Grid-graph modeling of emergent neuromorphic dynamics and heterosynaptic plasticity in memristive nanonetworks

Self-assembled memristive nanonetworks composed of many interacting nano objects have been recently exploited for neuromorphic-type data processing and for the implementation of unconventional computing paradigms, such as reservoir computing. In these networks, information processing and computing tasks are performed by exploiting the emergent network behaviour without the need of fine tuning its components. Here, we propose grid-graph modelling of the emergent behaviour of memristive nanonetworks, where the memristive behaviour is decoupled from the particular and detailed behaviour of each network element. In this model, the memristive behavior of each edge is regulated by an analytical potentiation-depression rate balance equation deduced from physical arguments. By comparing modelling and experimental results obtained on nanonetworks based on Ag NWs, the model is shown to be able to emulate the main features of the emergent memristive behaviour and spatio-temporal dynamics of the nanonetwork, including short-term plasticity, paired-pulse facilitation and heterosynaptic plasticity. These results show that the model represents a versatile platform for exploring the implementation of unconventional computing paradigms in nanonetworks.


Introduction
In a period where Moore's law is threatened by physical limitations, and classical von Neumann architecture determines a bottleneck for computing performances, the scientific community is moving towards new paradigms: one of the most promising is bioinspired and neuromorphic computing with the aim of emulating the highly-efficient human brain functionalities [1]. Memristive devices represent one of the most promising electronic devices to mimic functionalities of the biological synapse, a fundamental element in neuromorphic computing development [2,3]. Nowadays, besides conventional cross-bar array memristive architectures [4], self-organized memristive networks are attracting growing interest. These kinds of self-organized systems including NW networks, percolating nanostructured networks and atomic switch networks, are characterized by a structural topology that is more similar to biological neural networks and have been shown to exhibit structural plasticity, homo-and hetero-synaptic plasticity, paired-pulse facilitation (PPF) and collective neural-like dynamics [5][6][7][8][9][10][11][12][13][14][15]. These dynamics have been exploited for the implementation of unconventional computing paradigms such as reservoir computing [16][17][18][19][20][21]. While top-down realized crossbar arrays several models have been developed, the development of models able to describe the emergent behaviour of selforganized memristive networks at the macroscale still represent a challenge. In this context, Manning et al [22] tried to reproduce the exact nanowire topology with its junction distribution by exploiting SEM image processing. Each junction is represented as an edge in an abstract graph conceptualization, where each edge dynamically evolves according to a power law as a function of the current. Besides the one-to-one mapping of all the nanowire junctions, generating issues as soon as the system scales up, the model equation represents a semi-empiric formula, not deduced from physical first principles. Hochstetter et al [6] and Zhu et al [12] instead proposed, in their works, the generation of a random distribution of nanowires based on their mean length and variance, and the association of a threshold model equation to model each edge dynamics. This modelling approach, however, lacks an analytical solution due to the piecewise definition of its associate differential equation and considers one-to-one mapping of all the nanowire junctions, hindering the modelling of large scale NW networks. Also, all these models assume the same behaviour of all the NW junctions in contraposition with experimental results that show a distribution of resistance states and memristive response of network junctions due to the stochasticity of the mechanical contact.
Here, the emergent memristive behaviour of nanonetworks is modelled by means of a memristive gridgraph, where the dynamic behaviour of each edge is regulated by a physical-based potentiation-depression rate balance equation. The model is shown to be able to emulate the main features of the emergent behavior of Ag NW networks including spatio-temporal dynamics with short-term synaptic plasticity, PPF and heterosynaptic plasticity. As a consequence, the grid-graph model can be exploited for exploring neuromorphic-type data processing and for implementation of unconventional computing paradigms in nanonetworks.

Results and discussion
The proposed grid-graph model of memristive nanonetworks relies on two main steps: (i) approximation of a high-density memristive nanonetwork as a continuous medium, (ii) parcellation of the 2D continuous domain and subsequent approximation as a regular grid graph. As theoretically proven by Forrò et al [23] and experimentally shown by Sannicolo et al [24], for sufficiently dense conductive NWs, the voltage distribution over the sample is linear, resembling the electrical behaviour of a continuous medium. The approximation turns out to be valid even over few NW lengths scale for normalized nanowire density D > 2D c , where D c = 5.63 is the percolation critical density [25]. Also, previous works [26,27] on NW network conductivity mapping by electrical resistance tomography have outlined a nearly homogenous conductivity distribution across NW network by tailoring drop-casting. Experimental data of Ag NW networks [5] considered here have been measured on a sample (figure 1(a)) with an estimated normalized density D = 39.06, well above the percolation limit. In these devices, the emergent behaviour arises from resistive switching phenomena in NWs and in NW junctions [5]. The typical I-V characteristics of the Ag NW network showing the emergent memristive behaviour is reported in supplementary figure S1 (available online at http://stacks.iop.org/NCE/2/014007/mmedia). When a voltage difference is present in between two Ag NWs, an Ag conductive filament connecting Ag metallic cores can be formed across the polyvinylpyrrolidone insulating shell that surrounds the NW as a consequence of the Ag NW synthesis process. These Ag filaments are volatile, where the driving forces of the spontaneous dissolution are related to the minimization of the interfacial energy and electromotive force effects [28,29]. As a consequence, the volatile behaviour of Ag filaments reflects in short-term emerging dynamics of Ag networks.

Modelling nanonetworks as grid-graphs
The 2D plane representing the continuous material can be parcelled by an N × M mesh ( figure 1(b)). Experimental two-terminal measurements on these devices [5] exhibit a memristive behaviour highlighting shortterm plasticity (STP) property, hence the pixeled network can be represented as a regular grid of memristors that endow short-term memory dynamics ( figure 1(c)). Note that a similar parcelling of a continuous material followed by approximation as a grid of memristors or resistors has been performed in simulations to model a memristive material [30], as well as for modelling the insulating layer of a ReRAM cell [31,32]. Also, a regular square array of memristors have been considered from the theoretical point of view to show the possibility of memristive networks to solve shortest-path optimization problems [33] and to solve mazes [34]. Here, the grid graph model is built by associating a grid-graph node to each pixel of the parcelled domain, where memristive edges let the neighbour nodes communicate. Each grid-graph node was assumed to be equipotential. Moreover, in order to take into account the isotropic nature of homogenous deposited NW networks, diagonal memristive edges are added to the graph with a random orientation between diagonal nodes. The exclusion of an orientation for each possible diagonal is necessary to avoid the generation of a 3D graph, keeping the model as a 2D framework. Note that the proposed model considers the effective behaviour of the network over a large spatial scale, without mapping the exact topology of the network with all its nanowires and junctions, overcoming the impossibility of one-to-one mapping of the exact network topology in a spatial scale of the order of cm 2 through SEM images and to deal with the stochasticity of memristive behaviour of single NW junctions. Thus, the memristive behaviour of each edge of the grid-graph represents an average memristive behaviour resulting from a multitude of NW junctions in the corresponding area of the sample.

Memristive dynamics of grid-graphs
The emergent network dynamics, driven out of equilibrium by input voltage signals, depends on the evolution of each single memristive edge in the model design. Based on a previous work [35], where a model for mimicking STP in a single memristive cell was developed, each edge dynamics can be modelled by an a potentiation-depression rate balance equation: with g the normalized conductance and k p , k d the potentiation and depression coefficient respectively, defined as a function of input voltage only: Equations (2) and (3) derive from the physical description of the ionic diffusion [29,36,37], making the model a physics-based one. Defining the conductance range of each edge as [g min , g max ], the current flowing into it will follow the first Ohm's law: The low computational cost of the model is linked to the analyticity of the solution, which can be casted into a recursive fashion on a discretized time domain, emphasizing the memory property of the memristive edges: In this context, the model consists of 7 parameters (k p 0 , k d 0 , η p , η d , g 0 , g min , g max ) that can be retrieved from interpolation of experimental measurements. A similar approach was exploited in a previous work [35], where parameters of the potentiation-depression rate-balance equation were retrieved from experimental data in case of single NW memristive elements. Note that the same equation was exploited to describe dynamics of both diagonal and non-diagonal edges. The parameter g 0 is the initial conductance value of the memristive element from which the recursive solution (equation (5)) starts. In the original model conception, g 0 is physically associated to the pristine state of the memristive edge before the stimulation. Moving from single memristive cell to a network of them, the natural approach would be to generalize the concept of a single value g 0 to a tensor value g 0 to describe the initial state of each edge. However, usually, the initial conductance distribution across the physical nanowire network is not known. By assuming a uniform and homogeneous network in the pristine state, g 0 was assumed to be identical for each edge. The g 0 value for each edge was defined as the value such that the experimental pristine conductance measurement of the network (i.e. the effective conductance in between two points of the network) across two terminal is reproduced through data interpolation. This approach is further corroborated by the fact that the effect of initial state of the network vanishes after few iterations of the model solution. This simplified approach may also lead to a value of g 0 falling outside the range [g min , g max ], which is the result of decoupling the meaning of g 0 from a pure physical quantity.
While equation (5) describes each memristor evolution according to the input voltage on each single edge, the modified voltage node analysis [38] (MVNA) algorithm is exploited to compute, at each time step, the voltage distribution across the network by solving the Kirchhoff current law for each node of the network. The MVNA is a generalization of the voltage node analysis (VNA), offering the additional possibility to deal with voltage sources and branch currents. The voltage distribution then serves for the successive update of the network edges conductance. Remarkably, the MVNA does not deal with transients, but just looks at the circuit as a network of resistances (memristors frozen in a particular state) in a certain time instant, providing the regime solution of the circuit for each timestamp. The computational algorithm for modelling the memristive dynamics of the memristive network is conceptualized in figure 2.

Modelling emergent behavior of nanonetworks
Next, it is shown that the model is able to capture the main experimental features of the emergent behaviour of memristive neuromorphic nanowire networks including short-term homo-synaptic plasticity, PPF and heterosynaptic plasticity [5]. Homo-synaptic plasticity is an input-specific phenomenon, generating a change in the synapse strength (namely, the effective electrical conductance in the modelled device) in between pre-and post-synaptic neurons which are involved directly in the stimulation. When two or more identical electrical stimuli are close enough in time (temporally correlated), this kind of devices exhibits PPF [28,39]: the synaptic weight potentiation benefits from the facilitation of previous stimulus, progressively leading to a higher conductance state over stimulation pulses. Hetero-synaptic plasticity, instead, is an input-unspecific phenomenon where synaptic changes are observed to be induced also at synapses that were not directly stimulated [40]. For this purpose, the model graph has been defined as a 21 × 21 node structure with random diagonals added to a regular grid of memristors, whose dimension guarantee spatial positioning of electrodes according to fabricated devices reported in a previous work [5], by keeping relative distances among them. In the following, the position of electrodes reflects the geometry of the experimentally measured devices. Additionally, the extra constraint of no common edges between electrodes has been considered for the grid graph dimension definition. Figure 3 depicts the input-specific experimental data and the model reproduction of the network behaviour when stimulated with an input voltage pulse (4 V, 10 s) and following relaxation where the conductance state was probed with a reading voltage of 5 mV. As can be observed in figure 3(a), the grid-graph model is in good agreement with experimental results, both during potentiation and relaxation. Even though the good agreement between experimental data and modelling, it is necessary to point out that the model is not able to reproduce discrete steps in network conductance related to stochastic switching events in NW junctions. Such discrete steps, that can occur during both potentiation and relaxation, are hardly reproducible during the experiments, and they were reported also by other authors [9,22]. The recursive solution of the state equation in equation (5) does not take into account of such discrete changes of device conductance, since the model aims to reproduce the emergent collective behavior of the system (including potentiation/relaxation temporal dynamics and heterosynaptic phenomena). Comprising the stochastic behavior of single NW junctions would lead to a more accurate simulation, but probably at the cost of a remarkably increased computational time and power. Indeed, a key aspect of our results is that the equations exploited in our model can show the emergent behavior and dynamics of the network with a simple personal computer in the timescale of few minutes/hours. Indeed, the model represents a useful tool to investigate the conductance distribution evolution across the network and its voltage distribution ( figure 3(b)). The emergent behaviour of the collective memristors evolution results in the formation of a conductive path connecting source and ground electrodes. This is in agreement with direct experimental observation of a conductive path connecting stimulated areas by means of lock-in thermography, as reported by Li et al [41].

Effect of stimulation amplitude
The model allows also to investigate the effect of the voltage stimulation amplitude on the conductive path morphology. Figure 4(a) shows modelling and experimental results of the effective conductance potentiation and subsequent spontaneous relaxation by increasing the voltage stimulation amplitude in the range 1-4 V, respectively. As can be observed, in all cases the modelled effective conductance well interpolates experimental data. Also, the (b) panels report the corresponding obtained conductive path morphology after 10 s of stimulation (at the end of the stimulation, before relaxation). As can be observed, a higher voltage amplitude stimulation results in a strengthening of the conductive path. This is in agreement with experimental observations by Manning et al [22] that observed by means of passive voltage contrast SEM images an increase of the conductive pathway size by increasing the stimulation strength. The spatiotemporal growth and dissolution of the of the conductive path during stimulation with 1-4 V pulses and following relaxation can be visualized in supplementary movies 1-4, respectively.

Paired-pulse stimulation
Additionally, the grid-graph model was exploited for simulating the gradual increase of the effective network conductance upon stimulation with temporally correlated pulses to emulate PPF. Figure 5(a) shows both the input stream and the associated measured and modelled current upon network stimulation with a train of 20 pulses with a 10 V amplitude and a 1 kHz frequency. Thanks to the memory property, the potentiation of  In panel (a), the dashed line divides the stimulation region (in grey) and the following relaxation. In panel (b), red intensity of edges is proportional to the edge conductance, blue intensity of each node is proportional to the node voltage, black nodes represent input pads while arrows indicate the current direction. The left electrode was biased while the right one was kept as ground. a given pulse facilitates the potentiation of successive pulses, leading the current to increase more and more up to a saturation region. This gradual increase of the effective NW network conductance results from the competition in between conductive pathway formation during pulses and spontaneous relaxation in the waiting time in between pulses. As can be observed, the model clearly follows the same experimental trend. It is worth noticing that the model works well also on the millisecond timescale, beside the second timescale shown previously. Again, investigation on the conductive path formation can be performed, here as a function of the number of pulses sent to the network ( figure 5(b)), emphasizing the widening of the conductive region as the number of pulses increases.

Heterosynaptic plasticity
The spatial information provided by the model can be exploited to model hetero-synaptic plasticity, where a direct stimulation of the network across two terminal is responsible for changes in the effective conductance (b) Spatio-temporal evolution of the network conductance distribution as a function of the pulse number. Red intensity of edges is proportional to the edge conductance, blue intensity of each node is proportional to the node voltage, black nodes represent input pads while arrows indicate the current direction. The left electrode was biased while the right one was kept as ground.
across other pathways in between electrodes. Figure 6(a) shows the network schematization where five electrodes have been considered, in order to model the experimental results reported in reference [5]. An example of modelling hetero-synaptic plasticity is reported by direct stimulation of electrodes I and II. Here, by interpolating the input-specific curve associated to the potentiation of terminals I-II ( figure 6(b)), the conductive path generation after the end of stimulation can be observed between the involved electrodes (figure 6(c)). In order to assess the quantitative analysis of the model, correlation matrix of resistance variation after the stimulus has been evaluated both for experiment and the model (figures 6(d) and (e), respectively). Note that the resistance variation was calculated as ΔR = −(R post_stimulus − R pre_stimulus ) and the variation of ΔR has to be intended as a decrease in resistance of the synaptic pathway after stimulation. Besides the good reproduction of resistance variation just after the stimulus, the dynamics predicted by the model is in good agreement with the experimental measurement, as can be observed from figures 6(f) and (g). Comparisons of simulated with experimental data of heterosynaptic plasticity effects by direct stimulation of synapses I-III, I-IV and I-V are reported as supplementary figures S2-S4, respectively. It is worth noticing also that the here reported grid-graph dynamic model enables one to model variations of resistances in between electrodes quantitatively, while the previously reported static model [5] allowed only for a qualitative analysis.
Short-term dynamics and heterosynaptic plasticity effects of NW networks endow the system the capability of processing spatio-temporal inputs in multiterminal configuration for the implementation of reservoir Fitting of input-specific response across electrodes I-II due to an input pulse of 8 V amplitude and 1 s duration, followed by a long relaxation region at low input voltage. Inset shows details on the potentiation region. (c) Conductive path formation across stimulated electrodes. Red intensity of edges is proportional to the edge conductance, blue intensity of each node is proportional to the node voltage while arrows indicate the current direction and black nodes represent input pads. Pad I was biased while pad II was kept as ground. (d) Experimental and (e) modelled correlation matrix of resistance variation (ΔR) before and after stimulation of pads I-II (e). (f) Experimental and (g) modelled evolution of resistance variation (ΔR) before and after stimulation of pads I-II. In panels (d)-(g), directly stimulated synapses are highlighted in blue. computing, where the reservoir state is represented by the network conductivity map [16]. In this framework, the model that is able to reproduce both temporal dynamics and heterosynaptic behavior of the reservoir state for arbitrary input voltages (including pulses) can be exploited for simulating different strategies of computing implementation to solve both static and time-dependent tasks. Furthermore, despite non-idealities, it can even be beneficial for reservoir computing in order to enhance the extraction of relevant features, the effect of stochastic discrete steps-observed in experimental data, but not reproduced by the model-has to be evaluated.

Conclusions
In conclusion, we have reported grid-graph modelling of the emergent behaviour of memristive nanonetworks, where the dynamical evolution of the grid-graph edges is regulated by a physical-based potentiation-depression rate balance equation. In good accordance with experimental data obtained in multiterminal Ag NW networks, the grid-graph model is shown to emulate the emergent spatio-temporal dynamics of memristive nanonetworks, including short-term synaptic plasticity, PPF and heterosynaptic plasticity. These results show that the model can be exploited as a versatile tool for the implementation of unconventional computing paradigms in nanonetworks towards the realization of self-assembled neuromorphic systems able to process data similarly to our brain.