Neural Functional Connectivity Reconstruction with Second‐Order Memristor Network

The advances of neural recording techniques have fostered rapid growth of the number of simultaneously recorded neurons, opening up new possibilities to investigate the interactions and dynamics inside neural circuitry. The high recording channel counts, however, pose significant challenges for data analysis because the required time and computational resources grow superlinearly with the data volume. Herein, the feasibility of real‐time reconstruction of neural functional connectivity using a second‐order memristor network is analyzed. Spike‐timing‐dependent plasticity, natively implemented by the internal dynamics of the memristor device, leads to the successful discovery of temporal correlations between pre‐ and postsynaptic spikes of the simulated neural circuits in an unsupervised fashion. The proposed system demonstrates high classification accuracy under a wide range of parameter settings considering indirect connections, synaptic weights, transmission delays, connection density, and so on, and enables the capturing of dynamic connectivity evolutions. The influence of device nonideal factors on detection accuracy is systematically evaluated, and the system shows robustness to initial weight randomness, and cycle‐to‐cycle and device‐to‐device variations. The proposed method allows direct mapping of neural connectivity onto the artificial memristor network and can lead to efficient front‐end data analysis of high‐density neural recording systems and potentially directly coupled bioartificial networks.

timing-based learning rules. On the contrary, "traditional" memristors lack such internal timing mechanism, and would require external circuitry to convert relative timing information in the input to overlapping pulse widths or heights. [20] By using the STDP learning rule, applications such as coincidence detection [21] and temporal pattern recognition [22] have been recently experimentally demonstrated in a memristor system.
In this work, we analyze the ability of memristor networks to process neural recording data from a simulated large neural circuit. We show that a second-order memristor array has the potential to natively reconstruct the connectivity pattern from the recorded spike trains in an unsupervised manner, leading to inhomogeneous conductance distributions in the artificial network that match the structure of the biological networks. The high scalability, [14] parallel real-time processing, [13] and low power consumption [17] provide opportunities to integrate the memristor network with high-density neural recording systems for efficient front-end data analysis. [9] A large neural circuit consists of complex and possibly recurrent connections, and spurious results may arise due to the influence of indirect connections. Different transmission delays and neural jitters may also negatively affect the reconstruction accuracy, as STDP is very sensitive to timing information. To understand the impact of these factors on the system performance, we tested and verified the accuracy by comparing the device conductance map to the true biological connection map under a wide range of network settings, using data simulated in a randomly connected network of leaky integrate-and-fire (LIF) neurons. In addition, the memristor system allows the real-time capturing of dynamic connectivity evolutions in the neural work, in contrast to conventional techniques such as CC and TE which assume static network connectivity. Figure 1a shows the design of the memristor-based system for uncovering the connection pattern embedded in the neural network. The system allows direct coupling between the memristor synaptic network and the biological synaptic network, where spikes recorded from the simulated neural network are fed to the memristor network for analysis in an unsupervised fashion. As the synaptic connections are unknown, every recording site in the biological network is connected to both a row electrode and a column electrode of the memristor network, as shown in Figure 1a. The memristor conductance updates upon stimulation by the spikes, based purely on internal physics. [18] Specifically, the device we use here is a tantalum oxide-based second-order memristor. The first state variable (conductance) is determined by both the external stimuli and the status of a second Cross-correlogram computed from three pairs of neurons with excitatory, inhibitory, and no connection, respectively. c) Ground truth of the simulated two-layer network of LIF neurons. d) Ground truth synaptic weight map. e) Ground truth connection map. Connections induced by the excitatory and inhibitory connections are marked with magenta and cyan outline, respectively. f ) Conductance map of the memristor devices, measured after 600 s of recording duration. g) Conductance evolution from all devices in the memristor network over 600 s recording duration (%6000 spikes). state variable (internal temperature). When stimulated, the local temperature increases and then decays spontaneously after the removal of the stimulation as shown in Figure 2b. If a second spike arrives close to the first, heat effect from the first spike has not been completely dissipated and the oxygen vacancy drift and diffusion process will be accelerated by the higher temperature, leading to the growth or dissolution of conductive filaments. The device conductance change will thus strongly depend on the relative timing of the spikes applied to the top (corresponding to presynaptic neuron inputs) and bottom (corresponding to postsynaptic inputs) electrodes, natively following the STDP learning rule even though the rule is not explicitly enforced. [18] The bioplausible characteristics of the second-order memristor device lay the foundation for the discovery of temporal features hidden in the time series inputs with unsupervised learning.

Memristor Network for Neural Connectivity Reconstruction
The natively exhibited STDP behavior of the device at the initial conductance is shown in Figure 2c, obtained using a dynamic device model developed earlier. [18] Apart from the usual timingdependent plasticity, the conductance update is essentially zero for near-synchronous spikes in the device, as the two spikes (programming pulses) will overlap and cancel out each other. We note that in extracellular recordings near-synchronous spikes are also often undetected due to the fact that sorting algorithms normally fail to identify individual spike from the superposition of multiple spikes. [23] We first examine how neural connection can be obtained from the spike trains, by evaluating the impact of the firing of one neuron on another. A positive synaptic weight represents an excitatory connection, meaning that the firing of the presynaptic neuron i will increase the probability of the firing of postsynaptic neuron j. The impact is reflected in the cross-correlogram as shown in Figure 1b, where a peak is observed after a time lag τ against a relatively flat background, corresponding to a greater number of pre-and postsynaptic spikes arriving in sequence with time interval around τ. Therefore, if an excitatory connection exists between two neurons, the STDP learning rule will induce overall stronger potentiation than depression, gradually enhancing the device conductance over the recording period. On the contrary, for inhibitory connection, the firing of the presynaptic neuron will suppress the activities of the postsynaptic neuron, so the probability to observe a postsynaptic spike following a presynaptic spike will be lower than the baseline probability, represented by the trough in the CC plot. [4b] Consequently, the potentiation portion of STDP will become weaker, causing the device conductance to decrease over the training duration. In the case where no direct connection exists between the neural pair, the significant correlation of firing probability is absent, as represented by the essentially flat black curve in Figure 1b. The potentiation and depression effects induced by STDP will cancel out each other, leading to only small change or almost no change Figure 2. Native biorealistic implementation of STDP in second-order memristor devices. a) The efficacy of synapses can be modulated through the regulation of calcium ions, which affects the insertion or removal of α-amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid receptors (AMPARs) into the postsynaptic membrane. In the second-order memristor, the second-order state variable (e.g., internal temperature) can offer the internal dynamics to modulate ion movements and enable device conductance update. b) Schematics of the pre-and postspikes, the effective pulses applied on the device due to the pre-post pairs, and the device response. The rise and spontaneous decay of the internal temperature naturally encodes the timing information between pre-and postsynaptic spikes. The elevated internal temperature together with the input spike then induces the conductance change. c) STDP learning rule natively implemented in the device, following a physics-based device model. The insets illustrate the pulse scheme used to represent the neural spikes, showing the pre-and postsynaptic spikes (upper) and the effective voltage pulse V pre -V post (lower) on the memristor. in device conductance. Similar logic applies to the reverse direction (connection from neuron j to neuron i) but with opposite sign, where reverse excitatory connection will cause the conductance to decrease, while the reverse inhibitory connection will lead to positive conductance change.
To prove the feasibility of the proposed system, we then test the detection accuracy in a simple, simulated two-layer network of LIF neurons. The true synaptic weights and connection pattern are shown in Figure 1d-e. Apart from excitatory connections, each excitatory neuron also receives the input from one inhibitory neuron in the network. The synaptic connection strengths for the excitatory and inhibitory connection are set as 2 and -5 mV, respectively (characterized by excitatory/inhibitory postsynaptic potential [EPSP/IPSP] amplitudes). The neural spikes are simulated with the NEST simulator, [24] where the neurons fire at approximately 10 Hz and the simulation lasts for 10 min. When the spikes are applied to the memristor network, the device conductance is updated following the physics-based device model and will natively exhibit the STDP learning rule if the relative time interval between pre-and postsynaptic pulses is within AE20 ms. As the internal dynamics in the tantalum oxide memristor occur at a much faster time scale, in our simulation the spike trains are sped up by a factor of 10 000 in time to match the time scale of the memristor dynamics. We expect future developments with memristor devices offering similar internal time scales [9b] and biocomparable low operation voltage (%100 mV) [25] can allow direct coupling without any such preprocessing steps.
For this simple test case, the final memristor conductance matrix identifies both excitatory and inhibitory connections with 100% accuracy, as shown in Figure 1f. The conductance evolution in Figure 1g shows clear separation between three different cases-connected positive (excitatory þ inverse inhibitory), unconnected, and connected negative (inhibitory þ inverse excitatory). Although the conductance of individual devices fluctuates during the learning process, the accumulated effect can accurately reflect the connection information hidden in the spiking sequence as long as a sufficient number of spikes (N spike ) are fed into the memristor network. The stabilization of the system can be attributed to the self-regulating effects in the conductance update. When the device conductance is close to the maximum value, the conductance increase will slow down because the formation of conductive filaments has depleted the available conductive ions. The same applies to the opposite direction. The conductance depression will attenuate when the device is approaching the minimum conductance value because all the conductive ions have been driven back to the reservoir. This behavior resembles the multiplicative STDP learning rule, where the weight update is affected by a soft-bounded term, [26] and reflects the saturation effect of long-term potentiation and long-term depression observed in biological synapses due to the homeostatic effect of calcium ions. [27] We note that with the classical STDP rule, one can predict whether a connection exists or not from the conductance level, but cannot distinguish between excitatory and reverse inhibitory connections because both excitatory and reverse inhibitory (inhibitory and reverse excitatory) connections can lead to positive (negative) conductance changes, as shown in Figure 1f. These connections may be separated in later stages if needed through standard methods such as CC, TE, Granger causality, and generalized linear model.

Connectivity Reconstruction of a Randomly Connected Network
Following the feasibility test using the two-layer neural network, we verify the system performance in the analysis of a larger spiking network with randomly connected neurons. Again, we select the LIF model, which has been shown in many cases to provide a good approximation to the dynamics of more complex neuron models and real neurons. [28] The network is composed of 100 excitatory neurons and 100 inhibitory neurons. Each neuron is randomly connected with ten excitatory and ten inhibitory neurons with corresponding weights of J ex ¼ þ1.0 and J in ¼ À2.0 mV, as shown in Figure 3a. It needs to be noted that even if two neurons are not connected directly, there can be indirect connections in between for these larger networks. The sparsely connected random network of LIF neurons is simulated on the basis of Brunel's model. [29] Following the guidance of Brunel, [29] we set the parameters so that the neurons fire asynchronously with stationary global activity and irregular single cell firing, as can be seen from the raster plot in Figure 3e.
We randomly selected 15 excitatory and 15 inhibitory neurons as inputs to the memristor network, as in real cases one can only have limited access to a small part of the neural circuitry. Figure 3b shows the true weight map and the observed conductance map after applying around 6000 spikes. Again, the memristor conductance map successfully deciphers the direct connections from the spike trains. The histogram plots (Figure 3f ) display the evolution of the device conductance over time, where memristors corresponding to neuron pairs with true connections clearly emerge from the unconnected ones as more spikes are applied. The positive conductance changes are used to determine the presence of a connection because the positive conductance modulation is generally stronger in our devices (e.g., Figure 2c). The threshold value of the conductance change to classify a connection is selected using the Matthews correlation coefficient (MCC) metric to achieve an optimal balance between true positives (TPs) and false positives (FPs). [30] MCC takes account all four components (TPs, true negatives [TNs], FPs, and false negatives [FNs]) in the confusion matrix (METHODS), so it can provide a balanced measure even if the numbers of the positive and negative classes are disproportionate, such as the case studied here. Based on the MCC score, we chose þ7% conductance change as the threshold value and applied this to all the following tasks. Figure 3h shows that for the randomly chosen inputs, all the 76 TPs are detected after %3000 spikes are applied to the memristor network, and the number of FPs remains less than or equal to three over the entire simulation. As the firing of two neurons can be weakly correlated with each other through indirect interactions, FPs can be induced by the accumulated effect of the indirect connections, particularly in larger networks.
Considering the complexity and the heterogeneity of a real biological neural network, it is thus essential to study how different network characteristics would affect the reconstruction accuracy of the memristor-based system. In the following, we evaluate the influence of synaptic weight, firing rate, transmission delay, and www.advancedsciencenews.com www.advintellsyst.com connection density. We also test the robustness to temporal jitter and the capability to capture dynamic synaptic connections.

Influence of Synaptic Weight
In this analysis, we set the synaptic weight to different values that correspond to J ex ¼ þ0.5, þ1, þ2 mV and J in ¼ À1, -2, À4 mV. To make reliable inference on connection, sufficient N spike must be provided. Here, we investigate the performance of detection accuracy, as well as the required N spike for reliable detection at different synaptic weights. Results are shown in Figure 4c,e. It is not surprising to see that the required N spike decreases with the increasing synaptic strength. Analyzing the TP case also shows that most weak connections can be detected as long as the recording duration lasts long enough, e.g., when enough spikes are supplied. However, larger N spike also enhances the overall cumulative effect of indirect connections, resulting in more FP cases for weaker synaptic weights, as shown in the case of J ex ¼ þ0.5 and J in ¼ À1 mV in Figure 4b,d. Despite the undesired influence of indirect connections, the memristor network is able to identify most connections with low error rate even for relatively weak connection, as illustrated by the high MCC score in Figure 4c,e.
Synaptic connections with weight lower than þ0.5 or À1.0 mV are not covered in the simulation. One reason is that the detection of these very weak connections calls for very long duration of stable recording from a static neuron, which is not common in real experiments. In addition, studies [31] have shown that the distribution of synaptic weights in cortical circuitry follows a lognormal distribution, implying that the synaptic weight distribution is dominated by a few strong connections. In general, the local structure inside the brain can be viewed as a skeleton of strong connections immersed in a sea of weaker ones. Lining out the skeleton (strong connections) is therefore critical for the understanding of cortical interaction and dynamics.
We also checked the two FN cases for inhibitory connection at À4 mV. It turns out that these two cases have bidirectional inhibitory connections. The inhibitory effect will thus be partially canceled out by the connection in the opposite direction, making them harder to be detected. As in our tests a subset of neurons is randomly selected from the large neural network, the case at À2 mV happens to not have any bidirectional inhibitory connections and resulted in a slightly higher MCC, as shown in Figure 4e. Despite these small variations due to the random neuron selection process, both cases show high MCC scores after receiving sufficient number of spikes.

Influence of Firing Rate
Firing rate has long been believed to play an important role in information coding and network stability, as well as in neural response to addictive synaptic inputs or external stimuli. [32] To estimate the influence of firing rate, we simulate two kinds of situations: 1) all neurons in the network share roughly the same firing rate ranging from 2.5 to 20 Hz; 2) two groups of neurons fire at different rates in the network.
From the TP number results shown in Figure 5b, one can see that most connections are correctly predicted, except for only one www.advancedsciencenews.com www.advintellsyst.com missed case in the simulation with firing rate of 10-20 Hz. In most cases, the number of FPs is also almost zero, representing excellent classification performance. However, when the firing rate increases to 20 Hz, the number of FPs increases sharply to 10. This effect can be explained by indirect connections. If we assume independent Poisson processes for the spiking events, the baseline probability for the time interval between pre-and postsynaptic spikes to fall within the STDP window is proportional to ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r pre Ã r post p , where r pre and r post stand for the firing rate of pre-and postsynaptic neuron, respectively. The linear dependence on the firing rates can be clearly visualized in the cross-correlogram of two unconnected neurons, as shown in Figure 5a. Therefore, the effect of indirect connections becomes more significant at higher firing rates, generating a higher misclassification error. Despite the influence of indirect connections, the system still preserves high accuracy at 20 Hz with a TP rate of around 88%. The performance is further assessed by MCC together with two other scalar metrics-balanced accuracy (BACC) and F1 score (METHODS) in Figure 5c. All three metrics are found to be over 0.9 in all simulations, demonstrating that our system is generally applicable to neural networks with a wide range of firing rates typically encountered in neural recordings.

Influence of Transmission Delay
When modeling a realistic biological neural network, it is essential to consider the time delays for the signal to propagate  www.advancedsciencenews.com www.advintellsyst.com between two neurons, which has been shown to affect synchronization, neuron response, and spatiotemporal pattern generation. [33] The reliance on spike timing to uncover connectivity patterns in our approach means that it is crucial to assess the sensitivity of the system to transmission delay. In the following, we first simulate networks with different transmission delays in the range of 1-5 ms under the assumption of identical delay time between adjacent neurons within the network. To provide analysis in a more biorealistic situation, we then examine the case where the delay time has a random, i.e., uniform distribution between 1 and 5 ms. The conductance histograms after 4000 spikes for different transmission delays are shown in Figure 6a, showing the separation of connected positives from the rest in all cases. Unsurprisingly, the connected positives shift toward lower values with increasing transmission delay because the amplitude of device response decreases at longer timing delays as reflected by the STDP curve in Figure 2c. However, the delay (up to 5 ms) does not cause a measurable effect on the TP detection accuracy, while the FP numbers show some dependence on delay time. With peak response at around 1 ms, the system prediction has the highest number of FPs at t delay ¼ 1 ms because the higher magnitude of conductance update also leads to greater effect of indirect connections. Despite the minor influence on error rate, the overall system performance maintains high scores (>0.95) for the three metrics, as shown in Figure 6c.
For biological networks with delay up to tens of milliseconds, our systems can still be useable with simple adjustments. As the range of transmission delay inside the recorded network is normally predictable, one can alter the time scaling performed before applying the spikes to the memristor to match the scale of the propagation delay.

Influence of Temporal Jitter
Randomness in the firing of action potentials is inevitable for cortical neurons because of the stochastic nature of synaptic transmission and ion channels, as well as inhomogeneous spike transmission times between consecutive layers. However, both www.advancedsciencenews.com www.advintellsyst.com the dynamics of cortical networks and the active properties of individual neurons are able to overcome these potential problems and produce spikes with high timing accuracy. For example, recordings from frontal lobe reveal that the spiking of cortical neurons has very small temporal jitter on the order of 1 ms. [34] To mimic the neural jitter effect, we add a random noise drawn from normal distribution with standard variance of σ jitter (0.5-1.0 ms) to the spike timing. Figure 7aÀc shows that the peak/trough of excitatory/inhibitory connections in the CC plot spreads out after injecting temporal jitter. The influence of temporal jitter at different synaptic transmission delays is shown in Figure 7dÀf. Almost all TPs can still be correctly classified after the addition of temporal jitter, and only marginal changes in the number of FPs are observed. This effect can be explained by the fact that the overall device conductance modulation is a cumulative effect from a large number of neural spikes; therefore, the temporal perturbation added to each individual spike will be essentially averaged out. The fact that the system maintains high scores of the three scalar metrics further highlights its robustness against temporal jitter (Figure 7f ).

Influence of Connection Density in a Large-Scale Network
As mentioned earlier, the connection pattern of a realistic neuron network is essentially a skeleton consisting of sparse strong connections immersed in the background of weak connections. To study the influence of connection density on the reconstruction performance, we conducted simulations in a large sparsely connected neural network consisting of 5000 LIF neurons (4000 excitatory þ 1000 inhibitory). The synaptic weights are set as J ex ¼ þ0.2, J in ¼ À0.6 mV for weak connections and J ex ¼ þ1.0, J in ¼ À2.0 mV for excitatory and inhibitory connections, respectively. Our goal is to test whether the connection density will affect the identification of the strong backbone.
We first examine the influence of the connection density of weak synapses. The network is constructed in the way that each neuron receives inputs from 0.25% to 20% [number of incoming excitatory synapses per neuron (CE) ¼ 10-800] of excitatory population and 0.5% to 40% [number of incoming inhibitory synapses per neuron (CI) ¼ 5-400] of inhibitory population. The number of strong synapses is fixed at a low value of 50 for both excitatory and inhibitory connections.
The TP curve in Figure 8a shows that the strong connections need longer observation time to emerge from networks with higher density of weak synapses (CE ¼ 400-800, CI ¼ 200-400). After applying 6000 spikes to the device, part of the excitatory connections is still missing, contrary to the quick discovery in the sparse network (CE ¼ 10-100, CI ¼ 5-50). Further extension of the recording time helps to increase the number of TPs. As the connection density increases, neurons are sharing higher number of common inputs from the indirectly connected neurons, which diminishes the system effectiveness in connection detection. At higher connection densities, the assumption of sparse network is no longer strictly valid, such as at CE ¼ 800 and CI ¼ 400, where the dynamic stability of the network is impaired [35] and the autocorrelation plot is not all zero within the STDP window. This means there is a finite probability for a single neuron to generate two spikes within the window at CE ¼ 800 and CI ¼ 400, causing a measurable device conductance change. These two effects at high connection densities reduce the sensitivity of the system to the strong synapses, although the network is still able to make accurate predictions with a longer acquisition time. Figure 8c,d shows that the system maintains a low FP rate under 1% and a high MCC score of around 0.96 at CE ¼ 800 and CI ¼ 400.  We then evaluate the effects of strong connections, by fixing the total number of connections (200 excitatory, 100 inhibitory) and modulating the density of strong connections. Figure 8e,f shows that the density of strong synapse has negligible influence on classification accuracy and only has a minor effect on the required N spike for reliable reference. The FP rate and MCC score are shown in Figure 8g,h and demonstrate that the system offers similar performance when the strong connection density varies from 5% to 20% of fixed total connections. These analyses confirm that the system has good tolerance to connection densities of weak and strong synapses.

Detection of Dynamic Synapses
Activity-dependent modifications of synaptic efficacy are fundamental to the development and cognitive functions of neural circuits such as learning and memory. [36] Most standard statistical methods including CC and TE assume fixed neural connectivity, therefore cannot provide satisfactory performance for networks with underlying synaptic plasticity. To test the proposed system's performance on dynamic synapses, the following experiment is performed.
The simulation is composed of three segments with a total length of 30 min, where in every 10 min two existing connections (marked by red square) will be removed and two new connections (marked by green square) will be created. By comparing the conductance map with the true connection map at each simulation segment, the ability of the memristor system to capture new connection patterns can be tested. Figure 9a-e shows that following the native STDP behavior, the device conductance gradually evolves to match the new connectivity in the neural circuit. To compare with more conventional approaches, we also try to recover the connection with CC and TE. We find that CC generates three FPs and one FN, and TE produces three FNs because they are calculated from the spikes observed during the whole 30 min. All errors can be attributed to no longer existing connections or newly added ones, showing the limitations of these traditional methods to cope with nonstatic synapses. On the contrary, the adaptive feature of the dynamic memristor network offers the system the capability to evolve with new connection patterns, as shown in Figure 9d,e.

Device Nonidealities
The influence of device nonideal factors is then systematically studied. Here to speed up the simulation, we used a phenomenological multiplicative STDP learning rule [26] to numerically simulate the conductance updates, where the parameters are Figure 8. Influence of connection density on the detection of strong connections. The system performance is evaluated in a large neural network with 4000 excitatory and 1000 inhibitory LIF neurons. a-d) Influence of connection density of weak synapses. The evolution of true detection rate of a) excitatory and b) inhibitory connections is plotted against N spike . For higher connection density of weak synapses, larger N spike is required for reliable reconstruction. c) Evolution of FP rate showing its insensitivity to connection density of weak synapses. d) MCC score computed from predicted connections, including both excitatory and inhibitory synapses. eÀh) The influence of the connection density of strong synapses is evaluated by changing the proportion of strong synapses from 5% to 20% while fixing the number of total connections. Evolution of true detection rate of e) excitatory and f ) inhibitory connections shows the impact on detection efficiency. g) Evolution of FP rate demonstrating its insensitive to connection density of strong synapses. h) MCC score at different connection density of strong synapses.
www.advancedsciencenews.com www.advintellsyst.com based on quantitative fitting of the outcome from the physical device model. Regulation on weight update is implemented with a soft bound, where the synaptic change is weight dependent. As a result, synaptic depression becomes stronger for larger weights, in agreement with the device behavior when approaching saturation. The numerical simulations used the following parameters: A p ¼ 0.36, A d ¼ 0.3, τ p ¼ 4.8 ms, τ p ¼ 6.4 ms, and η ¼ 0.01. In addition, the effects of the initial weight distribution, and cycle-to-cycle (C2C) and device-to-device (D2D) variations have also been considered, by adding standard deviations from 5% to 20% to these parameters, corresponding to the range of experimentally observed C2C and D2D variations in physical devices. [37] w ¼ 1.5 ms < Δt < 20 ms (1) As the conductance update is dependent on the current status of the device, we included variations of the initial weight to the simulation. Examining the conductance evolutions in Figure 10a shows that the memristor conductance will eventually evolve to a stable distribution, independent of the initial conductance value as long as a large enough N spike is provided. This adaptive behavior of STDP learning rule enables the device to perform reliably regardless of its initial state, as has been discussed earlier.
Due to the stochastic nature of conductive filament growth, C2C variability is inevitable for memristor devices. We therefore simulate the influence of C2C variation on neural circuit connection inference. Figure 10b shows that no noticeable effects on inference accuracy can be observed, although the weight update processes become noisy. This observation can be explained by the fact that inference prediction is based on the cumulative behavior from thousands of spikes, where the effect of C2C variation is essentially averaged out.
Device-to-device (D2D) variation is then analyzed, by adding a Gaussian distribution to the amplitude (A p and A d ) and relaxation time (τ p and τ d ) of the STDP behavior. The box plot (Figure 10c,d) shows that the distribution of the device conductance in the three groups expands by a significant amount due to D2D variation. The connected-positive group starts to overlap with the unconnected group for larger variations, producing detrimental effect on the accuracy. The MCC score drops below 0.8 when the variation in amplitude (relaxation time) is equal to or greater than 15% (20%). However, when analyzing the more realistic D2D variation cases that keep the ratios between A p /A d and τ p /τ d fixed, we find that the system can effectively absorb the variations with no observable degradation in performance, as shown in Figures 10e,f. Figure 9. Connectivity reconstruction in a dynamic network. Two existing synapses (marked by red square) are pruned and two new synapses (marked by green square) are added to the network every 600 s. The simulation contains three segments, lasting 1800 s in total. a) True weight map and b) connection map at corresponding segments. c) Device conductance maps measured at the end of each simulation segment. Perfect match with the true connection map is obtained in all three segments. d) Ground truth synaptic connection values for the dynamic synapses and e) device conductance evolutions over time, illustrating that the memristor network is capable to learn the new synaptic weights. f,g) Connectivity reconstructions using f ) CC and g) TE mislabels the dynamic synaptic connections (FPs are marked by green outline, and FNs are marked by read outline). The threshold for CC and TE is selected by MCC.
www.advancedsciencenews.com www.advintellsyst.com These results suggest the system is robust to random initializations and C2C variations. D2D variations may significantly impact the classification accuracy, but if the potentiation/ depression ratio is maintained among devices, which appears to be feasible for practical devices, the system performance will not be noticeably affected.

Conclusion
In this work, we presented the concept of coupling a memristor synaptic network to a biological synaptic network to uncover the biological network connectivity in an unsupervised manner. The second-order memristor can natively produce the STDP behavior and decode the temporal dynamics of the spikes. The proposed system is used to decipher the temporal correlation hidden in spike recordings and capture the ground truth connection map in the network. The influence of neural network parameters on the system performance is comprehensively studied with synthetic neural spikes in a randomly connected network. The system reliably achieves high inference accuracy with MCC score higher than 0.9, showing great tolerance to a wide range of network conditions. The system also demonstrates advantages over standard statistical methods when dealing with dynamic synapses, due to the fact that it does not rely on the assumption of stationarity of neurons. The memristor conductance can naturally evolve following the changing firing patterns of the neural network. This ability allows the memristor network to be used to monitor the neural network connectivity evolution, which is crucial for the understanding of neural circuit development, normal and pathological brain functions. [38] Apart from synaptic connection reconstruction, the system can be utilized as coincidence detection [21] and temporal correlation detection between event-based datasets. [39] In earlier studies, [20] waveform engineering or pulse amplitude/duration modulation is required to implement STDP learning rule. The natural implementation of STDP learning rule in the secondorder memristor allows the use of simple, short, and nonoverlapping square pulse, greatly reducing the complexity for the design of the circuit.
We note that in real biological networks, the neuron firing is more complicated and can be influenced by fluctuations produced by external signals. The presence of large background noise may induce temporal correlations between neurons without direct synaptic connections and produce spurious predictions, i.e., FPs. One source of error comes from the indirect interactions via other neurons, especially in a large recurrent network. This might be reflected as a sharp peak in the CC plot, meaning that the neuronal firings are still correlated even in the absence of direct connections. Although the prediction does not indicate a true direct connection, it reflects the functional interactions between neurons; therefore, the interactions between neurons are still worthwhile for further analysis.
Another possible error is the slow fluctuations in the temporal correlation due to variations in the background noise. If the connections are strong enough, the temporal correlation will be significant enough from the baseline and can still be faithfully detected, even when subjected to the fluctuating background. As discussed in Section 2.2.1, the system will still be reliable in lining out the skeleton (strong connections), which is critical for the understanding of cortical interaction and dynamics. If weak synaptic connections also need to be reliably detected, more complicated methods such as spike jittering [40] or generalized linear model, [87c] which take the fluctuating part into consideration, may be used to conduct further examination on the detected neuronal pairs. The proposed system is still valuable www.advancedsciencenews.com www.advintellsyst.com in this scenario because it can remove the large amount of unrelated data and accelerate downstream analysis.

Experimental Section
Biological Neural Network Simulation: Neural simulations are conducted with NEST simulator. [24] The network is composed of excitatory and inhibitory populations. Connections between neurons are created with predefined synaptic weights. Transmission delay time is added to the synapse model when constructing the network. Apart from the synaptic inputs, a Poisson input is fed into the neurons independently as the external synaptic input. We also assign an external current source to the excitatory and inhibitory populations, respectively, to tune the neurons firing rate to the desired value. For randomly connected sparse networks, the parameter setting of the network follows the Brunel's model to ensure all neurons fires in the asynchronous irregular region.
Device Simulation: Device simulations are performed using the SPICE model developed in previous work. [18] The neural spike train pattern is scaled by a factor of 10 000 before being fed into the memristor array to match the time scale of the memristor dynamics. A presynaptic (postsynaptic) spike is composed of a short pulse with 100 ns duration and amplitude of þ1.01 V (À1.01 V), followed by a long pulse with 1.9 μs duration and amplitude of À0.05 V (þ0.05 V). The first short pulse can initiate both programming and heating effects. One purpose of the second, longer pulse is to mimic the refractory period of a biological neuron. The other purpose is to lower the amplitude of the programming pulse to minimize the programming effect induced by a single spike. The device conductance is monitored throughout the simulation. Median of the conductance of all devices is subtracted to eliminate the possible bias that may be induced by the asymmetric potentiation and depression effects.
Scalar Metrics: The MCC [30] is adopted to evaluate the performance of network reconstruction and select the positive threshold. As MCC takes account all four components (TP, FP, TN, and FN) in the confusion matrix, it can provide a balanced measure even when the number of two classes is not in proportion. MCC is computed with the following equation, where a coefficient of þ1 represents a perfect classification.
The BACC [41] is computed based on sensitivity (also known as TP rate or recall) and specificity (also known as TN rate). It simply takes the arithmetic mean of the two metrics. The weighted sum avoids inflated performance estimates on imbalanced datasets.
The F1 score is the harmonic mean of the precision and recall, where the precision is the number of TPs divided by the number of all positive items detected (TPs þ FPs), and the recall is the number of TPs divided by the number of all samples that should be identified as positive. [42] F1 ¼ 2 Recall À1 þ Precision À1 ¼ TP TP þ 1 2 ðFP þ FNÞ (4)