Small-worldness favours network inference in synthetic neural networks

A main goal in the analysis of a complex system is to infer its underlying network structure from time-series observations of its behaviour. The inference process is often done by using bi-variate similarity measures, such as the cross-correlation (CC) or mutual information (MI), however, the main factors favouring or hindering its success are still puzzling. Here, we use synthetic neuron models in order to reveal the main topological properties that frustrate or facilitate inferring the underlying network from CC measurements. Specifically, we use pulse-coupled Izhikevich neurons connected as in the Caenorhabditis elegans neural networks as well as in networks with similar randomness and small-worldness. We analyse the effectiveness and robustness of the inference process under different observations and collective dynamics, contrasting the results obtained from using membrane potentials and inter-spike interval time-series. We find that overall, small-worldness favours network inference and degree heterogeneity hinders it. In particular, success rates in C. elegans networks – that combine small-world properties with degree heterogeneity – are closer to success rates in Erdös-Rényi network models rather than those in Watts-Strogatz network models. These results are relevant to understand better the relationship between topological properties and function in different neural networks.


Introduction
Network Neuroscience seeks to unravel the complex relationship between functional connectivity in neural systems (i.e., the correlated neural activity) and their underlying structure (e.g., the brain's connectome); among other goals [1][2][3][4][5] .The functional connectivity is responsible for many tasks, such as segregation, transmission, and integration of information 6,7 .These (and other) tasks are found to be optimally performed by neural structures that show small-world properties 2,4,6 , which range from human brain connectomics 5,8,9 to the Caenorhabditis elegans (C.elegans) nematode neural networks [10][11][12] .In general, structural networks have been revealed by tracing individual neural processes, e.g., by diffusion tensor imaging or electron-microscopy -methods that are typically unfeasible for large neural networks.On the other hand, functional networks are revealed by performing reverse engineering on the time-series measurements of the neural activity, e.g., by using EEG or EMG recordings.These methods are known as network inference and are mainly affected by data availability and precision.
In general, network inference has been approached by means of bi-variate similarity measures, such as the pair-wise Cross-Correlation 9,13,14 , Granger Causality [15][16][17] , Transfer Entropy [18][19][20] , and Mutual Information [21][22][23][24] , to name a few.The main idea behind the similarity approach is that, units sharing a direct connection (namely, a functional or structural link exists that joins them) have particularly similar dynamics, whereas units that are indirectly connected (namely, a functional or structural link joining them is absent) are less likely to show similar dynamics.Although intuitive, this approach has found major challenges in neural systems due to their complex behaviour and structure connectedness, resulting in highly correlated dynamics from indirectly connected units and loosely correlated dynamics for directly connected units.Moreover, because most works have focused on maximising the inference success (in relation to its ability to discover the structural network) and/or optimising its applicability 22,[24][25][26] , we are still unaware of which are the main underlying mechanisms that affect the inference results.Namely, differentiating the underlying structure with the functional connectivity -particularly with respect to establishing which of the different network properties are mainly responsible for hindering inference success rates.
In this work, we reveal that the degree of small-worldness is directly related to the success of correctly inferring the network of synthetic neural systems when using bi-variate similarity analysis.Our neural systems are composed of pulse-coupled Izhikevich maps and connected in network ensembles with different small-worldness levels -but statistically similar to the C. elegans neural networks.The inference process is done by means of the pair-wise cross-correlation between the neurons' activity.We assess the inference effectiveness by means of receiver operating characteristic (ROC) analysis and, in particular, the true positive rate (T PR), which we show is the only relevant quantity under our inference framework.Our findings show that the T PR peaks around a critical coupling strength where the system transitions from synchronous bursting dynamics to a spiking incoherent regime.Specifically, we find that the highest T PR is for networks with significant small-worldness level.We analyse these results in terms of different topology choices, collective dynamics, neural activity observations (inter-spike intervals or membrane potentials), and time-series length.We expect that these results will help to understand better the role of small-worldness in brain networks, but also in other complex systems, such as climate networks [27][28][29] .

Results
We infer the underlying network of a synthetic neural system by creating a binary matrix of 1s and 0s from the pair-wise cross-correlation (CC) matrix of the signal-measurements.The resultant binary matrix represents the inferred connections that the neurons composing the system share, which we obtain by applying a threshold to the CC matrix.The threshold assumes that a strong [weak] similarity in the measured signals, i.e., a CC value above [below] the threshold, correspond to a 1 [0] in the inferred adjacency matrix, suggesting that a direct [indirect] structural connection exists.In spite of this (seemingly) over-simplification, this binary process is broadly used in network inference 4,5,9,22,24 and it tends to keep the most relevant information from the underlying connectivity.Moreover, when the underlying network is known, it allows to quantify how poorly or efficiently the bi-variate method performs in terms of the receiving operation characteristic (ROC) analysis [30][31][32] .In particular, we set the threshold such that the inferred network has the same density of connections, ρ = 2M/N (N − 1), as the underlying structure, where N is the network size and M is the number of existing links.This means that we assume an a priori (minimal) knowledge about the underlying structure, namely, we require knowing ρ in order to choose the threshold.
The true positive rate, also known as sensitivity, is the proportion of correctly identified connections with respect to the total of existing connections 30 , i.e., T PR ≡ T P/ (T P + FN), where T P is the number of true positives and FN is the number of false negatives.This quantity is part of the ROC analysis, which includes the true negative rate, T NR, false positive rate, FPR, and false negative rate, FNR = 1 − T PR.Taken together, these variables quantify the performance of any method.However, when fixing the inferred network's density of connections, ρ, to match that of the underlying network, we can show that the T PR is the only relevant variable in the ROC analysis -all remaining quantities can be expressed in terms of the T PR and ρ.For example, the T NR = 1 − FPR, also known as specificity, can be expressed in terms of the T PR and ρ by where we use the fact that T P + FN = M is the number of existing connections, is the density of connections, and T P + FP = M is the number of connections we keep fixed for the inferred matrix in order to maintain ρ invariant.Hence, as a result of our threshold choice, FN = FP, implying that the FNR ≡ FN/ (FN + T P) is identical to the false discovery rate, FDR ≡ FP/ (FP + T P), and that the precision, PPV ≡ T P/ (T P + FP) is identical to the T PR.Overall, these relationships mean that in order to quantify the inference success or failure, we can solely focus on studying how the T PR changes as the dynamical parameters and network structure change.
The following results are derived from time-series measurements of pulse-coupled Izhikevich maps interacting according to different network structures and coupling strengths, where each map's uncoupled dynamic is set to bursting (see Methods for details on the map and network parameters).Pulse coupling is chosen because of its generality, which has been shown to allow the representation of several biophysical interactions [33][34][35] , and single parameter tuning, i.e., the coupling strength, ε.In particular, we register the neurons' membrane potentials (signals coming from the electrical impulses) and inter-spike intervals (time windows between the electrical pulses) of 10 randomly-set initial conditions and T = 7 × 10 4 iterations, from which we discard the first 2 × 10 4 iterations as transient (we also analyse the effects of keeping shorter time-series).
Without losing generality, we restrict our analyses to connecting the neurons (maps) in symmetric Erdös-Rényi (ER) 36 and Watts-Strogatz (WS) 37 network ensembles that have identical size, N, and sparse density of connections, ρ, to that of the C. elegans frontal and global neural network [10][11][12] .Specifically, when constructing the ensembles we set N f = 131 with ρ f 0.08 (which corresponds to setting a mean degree, k f 10.5), or N g = 277 with ρ g 0.05 (which corresponds to setting a mean degree, k g 13.8), respectively.The reason behind this choice is that the C. elegans neural networks are one of the most cited examples of real-world small-world networks 5,[10][11][12]37 , showing small average shortest paths connecting nodes and high clustering. Naely, the C. elegans neural networks have a high small-worldness coefficient σ , defined as the normalised ratio between the clustering coefficient and average path length 38,39 , but also show an heterogeneous degree distribution.More importantly, these network ensembles constitute a controlled setting where to compare and distinguish the main topological factors favouring or hindering the inference success, providing us with a reproducible framework to modify the network properties within each realisation.
We find that the resultant average sensitivity from these network ensembles is more significant, robust, and reliable on WS ensembles than on ER ensembles, pointing to a fundamental importance of the underlying small-worldness for a successful inference.In particular, Fig. 1 shows the resultant success rates for ER (dotted lines with unfilled squares) and WS (dotted lines with unfilled diamonds) ensembles, plus, a comparison with the results we obtain when using the C. elegans (CE) neural frontal (left panel) and global (right panel) network structure (continuous lines with filled circles).Specifically, Fig. 1(a) and (b) show the ensemble-averaged T PR results for N = 131 and N = 277 pulse-coupled Izhikevich maps, respectively, as a function of the coupling strength, ε, between the maps.From both panels we also note that the CE overall resultant success rates are closer to the ER ensemble-averaged T PR results than to the WS ensemble-averaged T PR results -in spite of the CE small-worldness coefficient for the N = 131 networks being the same as the WS, σ = 2.80.The results in Fig. 1 show how important the underlying degree distribution and small-worldness are in the generation of collective dynamics that can be analysed by means of a bi-variate inference method with a sufficiently high success rate.map parameters are set such that the isolated dynamics is bursting (see Methods).The T PR values for the ER and WS come from averaging the results over 10 initial conditions and 20 network realisations with similar topological properties to that of the CE (i.e., number of nodes, average degree, and density of connections).For the CE, the results are averaged only on the initial conditions.The T PR is found by comparing the true underlying network with the binary matrix obtained from cross-correlating the membrane potential time-series (T = 5 × 10 4 iterations) and fixing a threshold such that the inferred density of connections ρ f matches that of the CE, with ρ f 0.08 in panel (a) and ρ f 0.05 in panel (b).The horizontal dashed line in both panels is the random inference T PR, namely, the null hypothesis.
In order to critically explore the significance that the underlying small-worldness has on the resultant inference, we fix the degree distribution and density of connections as we increase [decrease] σ in each of the 20 underlying ER [WS] network realisations using the rewiring method proposed in Ref. 40 .Figure 2 shows the resultant network inference -quantified by the T PR(ε, σ ) -after we make the isolated changes in the small-worldness coefficient, σ , of the underlying structure for the N = 131 pulse-coupled Izhikevich maps (similar results are found for N = 277).The ensemble-averaged inference results (colour coded curves) that we get from making this topological change to σ on the underlying ER and WS networked system are shown in Fig. 2(a) and (b), respectively.We can see from these panels that the highest T PR values are achieved for the largest σ values, meaning that the best inference happens for networks with large σ .Also, we can see that there is a broad coupling strength interval (0.25 ε 0.33) for both network classes that allows us to infer better than making blind random inference (dashed horizontal lines).From these panels, we note that network inference effectiveness increases robustly (namely, regardless of parameter changes) and significantly (namely, reliably across ensembles and random initial conditions) as the small-worldness, σ , of the underlying structure is increased -whilst keeping its density of connections and degree distribution invariant.Consequently, in order to increase the inference success rates in the sparse ER networks, we need to increase the local clustering inter-connecting the maps.On the contrary, WS networks show optimal inference efficiency without modifying their clustering because of their inherent large small-worldness coefficient.
The coupling strength ε * 0.26, which maximises the T PR, implies an average impulse per map of ε * / k 0.025.This coupling is associated to a collective regime with a loosely coupled dynamic, as we show in Fig. 3, since it corresponds to an average 2.5% synaptic increment of the membrane potential range (i.e., the difference between maximum and minimum membrane potential values) due to the action from neighbouring maps.The fact that our inference method recovers approximately 50% of all existing connections at ε ∼ 0.26 (T PR = 0.5), implies that the FNR = 0.5, and from Eq. ( 1), this also implies that T NR 0.96 and that FPR 0.04 for ρ 0.08.This means that for ε ∼ 0.26, the inference method is highly efficient in detecting true nonexistent connections (T NR → 1) and falsely classifying these connections as existing ones (FPR → 0).This efficiency is a consequence of sparse networks having more non-existing connections than existing connections (N (N − 1)/2 M); as it happens in our ensembles.Hence, in sparse networks the challenge is to correctly identify  In spite of the similar results in Fig. 2 panels, we can distinguish a slight advantage in the ER networks ensemble-averaged T PR values [Fig.2(a)] over the WS T PR values [Fig.2(b)], which also appear when N = 277.We understand that differences in the inference results have to appear because of the dependence on the underlying degree distribution (as well as in the small-worldness), as it is observed in the results for the CE and ER networks shown in Fig. 1.In Fig. 2, the degree distributions correspond to those of the ER and WS networks respectively, which are kept invariant as the small-worldness of the underlying network is changed.However, the similarity in the results from Fig. 2(a) and (b) (for similar small-worldness values) can be explained due to the finite size systems, which make the ER and WS degree distributions similar (a similar behaviour is also observed for N = 277 -not shown here).We also note that the modified ER networks T PR results narrowly outperform WS inference results, where success rates reach values higher than 50% for coupling strengths close to ε * 0.26.These T PR values are significantly higher than making a blind random inference of connections (dashed horizontal lines), which successfully recover only 8% of the existing connections.
All T PR curves share an abrupt increase in the success rate around a critical coupling strength of ε * ≈ 0.26.Figures 1 and  2 show this abrupt jump in the inference success for all networks analysed -though the exact value of ε * may vary slightly for each topology realisation.This sudden increase in the success rates points to a drastic change in the systems' collective dynamics as the coupling strength is increased beyond ε * .In order to analyse the maps' collective dynamics as a function of ε and the underlying topology, we compute the order parameter, R, defined as the time-average of the squared difference between two inter-spike intervals (ISIs) time-series (i.e., the series of time differences between two consecutive spikes) summed over all pairs of ISIs.Specifically, R = ∑ j<i R i j , with R i j = (T i − T j ) 2 t , where T i is the i-th neuron ISI time-series and • t is the time-average.This means that the R value is high [low] when the time series are different [similar].
In Fig. 3(a) we show how the order parameter R changes with the coupling strength, ε, for different systems with N = 131 neurons.Namely, the results for the C. elegans frontal neural network structure are shown by the filled (black online) circles, whereas for the ER and WS networks, the ensemble-averaged and initial-condition averaged R values are shown by unfilled (blue online) diamonds and unfilled (red online) circles, respectively.The ensemble averages are found from 20 realisations and from 10 different initial conditions for each topology realisation.For the CE network each R value represents an average over 10 initial conditions.We can observe that close to ε ≈ 0.25, R decreases abruptly for all network structures falling 2 orders in magnitude.This drop corresponds to a switch from a collective bursting regime to a lightly disordered spiking regime which has higher collective coherence.This spiking regime appears as disordered when compared to the apparently synchronous bursting dynamics, but it is in fact partially coherent -a particularly suitable condition to perform a successful network inference 22 .For example, Fig. 3  WS degree distributions but with different small-worldness levels (as previously described).This means that the collective dynamics' abrupt change also happens for the modified networks, namely, the networks modified by our rewiring process to increase or decrease their overall small-world coefficient.Panels (e) and (f) show the resultant raster plots for ε = 0.23 and 0.26, respectively, for a realisation of an ER network with N = 131 maps and σ = 2.1.Furthermore, we can see from Figs. 1 and 2 that the sensitivity falls rather smoothly for all networks as we increase ε beyond the critical value ε * .The reason behind this smooth change is that, as ε increases beyond ε * , the neurons gradually begin to fire in a more ordered spiking, namely, achieving synchronisation.Thus, partial coherence between the time-series vanishes and inference becomes impossible.The smooth decrease in sensitivity can be observed by the rate in which the order parameter decreases for ε larger than ε * , as in Figs. 3 panels (a) and (d).In general, the neural systems we analyse stay in a partially coherent spiking regime for an interval of coupling strength values (approximately between ε ≈ 0.25 and 0.30), where network inference T PR values remain above the random line.
So far we have shown that increasing small-worldness favours network inference, obtaining success rates that appear robust to changes in the degree-distribution type (e.g., ER, WS, and CE), initial conditions, and similar for a broad coupling strength region.However, we can see from Fig. 1 that as the N increases from 131 (panel (a)) to 277 (panel (b)), the T PR drops significantly for the C. elegans networks.The reason for this drop comes from the broadness in the global CE neural-network's degree-distribution.As we can see from Fig. 4, when N = 131 (panel (a)), all degree distributions are somewhat similar and narrow, but when N = 277 (panel (b)), the CE topology shows the presence of hubs and a long tailed distribution.This is why on Fig. 1(b), the T PR results for the N = 277 WS network ensemble are extremely similar to those T PR values when N = 131 in Fig. 1(a).Similarly, we can see the same resemblance in the T PR results for ER networks, which also hold a narrow degree distribution, as shown by the dashed curves in Fig. 4. On the contrary, the significant differences emerging from the CE degree distributions for N = 131 and N = 277 impact directly into the inference success rates.This leads us to believe that heterogeneity in the node degrees hinders network inference.We find that the former results are also robust to changes in the time-series length.Moreover, our conclusions regarding small-worldness favouring inference still hold if one chooses a different time-series representation for the neural dynamics of each map.Namely, our findings hold for inter-spike intervals (ISI) as well as membrane potentials (see Supplementary Information for further details).Remarkably, we see that when using ISIs to measure the neural activity, inference success rates are significantly lower than when using membrane potentials -regardless of the particular topology or collective dynamics.However, these T PR values are more robust to changes in the topology realisation (i.e., different structures constructed with the same canonical models and parameters).This leads us to conclude that, using ISIs instead of membrane potentials, allows us to achieve worst inference success rates but with a more reliable outcome.

Discussion
In this work we have shown that the network inference methods based on the use of cross-correlation (CC) to measure similarities between components are more effective when inferring small-world structures than other types of networks.This conclusion has broad implications, since cross-correlation is widely used to reveal the underlying connectivity of neural systems, such as the brain, and to gain information about long-range interaction in climate networks.
We have shown that, for networks with similar degree distributions, the small-worldness level is the main topological factor affecting the inference success rates.The results shown in Figs. 1 and 2 account for the reliability and robustness of this conclusion.Also, in Fig. 1 we observe that our inferring method is consistently and robustly more successful when inferring Watts-Strogatz networks than the C. elegans (CE) neural structure, despite both having the same small-worldness level.This points to the effect that degree-distribution range has on the inference success, namely, the broader the degree distribution, the less effective the inference is.
Our results show that the appearance of highly connected nodes (hubs), such as in the CE global network, is another important factor, hindering successful inference.In other words, we find that success rates are generally lower when inferring networks which have higher degree heterogeneity.This finding is relevant because real small-world networks, such as the C. elegans neural structure, often combine small-world properties with other complex features -such as the presence of hubs, hierarchies or rich clubs -, resulting in a higher degree heterogeneity with respect to the canonical Watts-Strogatz network model.In particular, hubs have been related to the phenomenon of hub synchronisation in brain networks 41 and scale-free networks 42 .Hub synchronisation is particularly detrimental for network inference, since it leads to strong correlations between non-connected nodes and weak correlations between the hubs and their neighbouring nodes.
Regardless of the underlying structure, the coupling strength range that allows for a successful network inference is located after a critical value, in which the system transitions from a collective dynamics with an apparently synchronous bursting regime to collective asynchronous spiking.This transition has been reported to take place in C. elegans neural networks 43 and corresponds to a partially-coherent state, which has been shown to be necessary for having successful network inference 22 .Figure 3 shows an abrupt change in the systems' order parameter around the critical coupling strength, revealing a coherence-loss after the transition.Although we can only perform successful network inferences when the systems are in partial coherence, it is reasonable to assume that real neural systems are in such states, since they consistently transition between synchronous and asynchronous states in order to perform different tasks and cognitive functions 1-6, 43, 44 .

Synthetic neural network model
Our synthetic neural model is the Izhikevich map [33][34][35] , which belongs to the bi-dimensional quadratic integrate-and-fire family.This map consists of a fast variable, v, representing the membrane potential, and a slow variable, u, modelling changes in the conductance of the ionic channels.One of the main advantages of using Izhikevich maps is that it combines numerical efficiency (inherent to map-based models) with biological plausibility 34 .The isolated map equations of motion are given by When different values of I, c, and d, are fixed, the Izhikevich map can show extensive dynamical regimes, which have been observed in real neurons.We set the parameter values such that the regime exhibits bursting dynamics.Namely, d = 0, c = −58, and I = 2.However, when Izhikevich maps interact, the resulting single-neuron dynamics can differ significantly from the bursting regime.
The interactions are set to be pulsed via the fast variable, v, and controlled by a global coupling-strength order parameter, ε.This pulse-coupling type is able to represent many real neural interactions.With this model, every time a neuron spikes it sends a signal to the adjacent neurons (i.e., to all neurons that are connected to it), instantly advancing their membrane potentials by a constant value.Specifically, the dynamics for the n-th neuron is given by , if v i,n < 30 mV, and where k i is the i-th node degree (i.e., its number of neighbours), A i j is the i j-th entry of the network's adjacency matrix, ε is the coupling strength, and δ (x) is the Kronecker's delta-function.We iterate Eq. (3) 7 × 10 4 steps from 10 random initial conditions for each topology and coupling strength, removing a transient of 2 × 10 4 steps.In particular, the coupling term, 30), acts as follows.If a connection between neurons i and j exists, then A i, j = 1 and neuron i receives an input of value ε/k i every time neuron j reaches the threshold v j,n = 30.Otherwise, the neuron remains unchanged.

C. elegans neural structure
We use data from Dynamic Connectome Lab 10,11 to construct the C. elegans' frontal (N = 131 nodes) an global (N = 277 nodes) neural networks.These networks are represented by weighted and directed graphs, that we simplify by considering the unweighted non-directed versions.The reason behind this choice is to bring up-front the role of structure alone into our systems' collective behavior.Under these symmetric considerations, we find that the network's mean degree is k f = 10.5 (k g = 13.8) for the frontal (global) connectome, with a sparse edge density of ρ f = 0.08 (ρ g = 0.05).The average shortest-path length for our C. elegans frontal (global) neural network is l f = 2.5 (l g = 2.6) and its clustering coefficient is C f = 0.25 (Cg = 0.28).Both neural structures have short average path lengths (similar to Erdös-Rényi networks with equal edge density) and high clustering coefficients (similar to Watts-Strogatz networks with equal average degree), hence, they show small-world properties.The small-worldness coefficient of our C. elegans frontal (global) neural network is σ f = 2.8 (σ g = 5.1), which falls within the expected small-world range (σ > 2).

Network ensembles
In order to study the role that small-worldness and degree heterogeneity have in the network inference results, we build two ensembles of 20 Erdös-Rényi (ER) adjacency matrices with N = 131 and N = 277 nodes respectively, and two ensembles of 20 Watts-Strogatz (WS) adjacency matrices.We choose the number of nodes in our network ensembles to match the sizes of the C. elegans (CE) frontal (N = 131 nodes) and global (N = 277) neural structures.We also tune the algorithms to build these networks such that they have the same edge densities as the CE frontal and global neural networks, namely, ρ = 0.08 and ρ = 0.05, respectively.In addition, the WS ensemble is also tuned such that the algorithm parameters produce networks with similar small-worldness levels to that of the CE networks.In what follows, k denotes average among network ensembles, while k expresses the average among nodes of a single network.
Our ER ensemble is built with a probability to linking nodes in each network of p = 0.08 (0.05) for N = 131 (N = 277) -values which are above the percolation transition.These probabilities yield mean degrees of k = 10.5 ( k = 13.8) for the N = 131 (N = 277) networks.In both cases, the variability within the ensemble of these mean degrees is σ k = 0.3.We can corroborate that the nematode's neural networks also have mean degrees falling within one standard deviation of the ER ensemble-averaged mean degrees.The clustering coefficient in the ER model is usually low (C = p in the thermodynamic limit), being C = 0.08 ( C = 0.05) for our N = 131 (N = 277) network ensembles.The ER networks also hold a small average shortest-path length.In our ensembles, the shortest-path lengths are l = 2.3 ( l = 2.4) for the N = 131 (N = 277) networks.Correspondingly, the averaged small-worldness levels of our ER ensembles are σ = 1.0 in both cases -as expected -, indicating the absence of small-world effect.
The Watts-Strogatz (WS) algorithm takes an initial ring configuration in which all nodes are linked to K/2 neighbours to each side, and then rewires all edges according to some probability p 37 .Using this model, we construct a network ensemble with link density and small-worldness levels similar to the CE neural networks.In particular, we choose a mean degree and rewiring probability that yields similar average path lengths and clustering coefficients as the CE neural structures.Specifically, we fix the mean degree at an integer value, namely, k = 10 ( 14) for N = 131 (277) nodes.Then, for each rewiring probability p, we generate 20 adjacency matrices and calculate the mean average path length and clustering coefficients.Thus, for each rewiring probability p we have a point in the [C, l ] space.This allows us to choose the rewiring probability, p , which holds the closest point in the [C, l ] space to the CE networks values.For such p , our ensembles have shortest-path lengths of l = 2.5 ( l = 2.6), and clustering coefficients of C = 0.4 (C = 0.28) when N = 131 (N = 277).Correspondingly, our networks' small-worldness levels are σ = 2.8 ( σ = 5.3) when N = 131 (N = 277).These values are similar to the CE small-worldness levels and indicate the presence of the small-world effect.neural networks; with map parameters set such that the isolated dynamics is bursting (see Methods).The T PR values for the ER and WS come from averaging the results over 10 initial conditions and 20 network realisations with similar topological properties to that of the CE (i.e., number of nodes, average degree, and density of connections).For the CE, the results are averaged only on the initial conditions.The T PR in each case is found by the same process as Fig. 1, but using the interspike intevals time-series (considering 2000 spikes) instead of the membrane potentials.The horizontal dashed line in both panels is the random inference T PR, namely, the null hypothesis.
In order to gain insight on the differences between the use of interspike intervals and membrane potentials to perform network inference process, another key factor is the robustness of the inference sensitivity among different realisations of the topology.We can assess this aspect by comparing the standard deviation of the sensitivity, σ T PR , in our N = 131 and N = 277 ER and WS network ensembles, when we perform inference processes using membrane potentials and interspike intervals.Fig. 6 shows σ T PR as a function of coupling strengths in the N = 131 and N = 277 ER and WS network ensembles.We observe that, in all cases, the sensitivity's standard deviation is considerably higher when using membrane potentials than when employing interspike intervals, for most coupling strengths.Combining these results with those shown in Fig. 5, we conclude that while using the interspike intervals (instead of the full membrane potentials information) hinders our network inference process, it provides a higher robustness.Hence, there is a trade-off between efficiency and robustness when choosing different neural behaviour representations.

Different time-series length
In all our network inference processes we used a fixed time-series length (T = 5 × 10 4 when using membrane potentials and T = 2 × 10 3 when using interspike intervals).To check the robustness of our results against changes in the length of the  12/13 time-series used to calculate cross-correlations between pairs of neurons, we calculate the same T PR(ε) curves shown in Fig. 1, but using different time-series lengths.Fig. 7 shows a comparison between the sensitivity as a function of the coupling strength for Erdös-Rényi networks, using different time-series lengths.We observe a great similarity between all the T PR(ε) curves, with only quantitative differences when the time-series is too short.We obtain similar results in the Watts-Strogatz and C. elegans case.Therefore, we can confirm that our findings are independent of the time-series length considered in our calculations.

ROC Analysis
The analysis of the Receiver Operating Characteristic (ROC) curve is a potent tool for evaluating the efficiency of a classification algorithm [30][31][32] .It has been used, for example, to assess the efficiency of clinical test for distinguishing between infected and healthy patients 47 , or for testing machine-learning-based classification algorithms 48 .In the network inference scene, ROC analysis is often used to test and optimise inference methods 22 .To assess the efficiency of an inference technique, ROC analysis relies on four fundamental quantities: the Sensitivity, or True Positive Rate (T PR, the fraction of present links that were correctly identified as such), Specificity, or True Negative Rate (T NR, the fraction of absent links that were correctly identified as such), False Positive Rate (FPR, the fraction of absent links that were incorrectly identified as present), and the False Negative Rate (FNR, the fraction of existing links that were incorrectly identified as absent).
In the network inference context, these four quantities are naturally related by two constrains: the total number of links (N[N − 1]/2 in a network with N nodes) and the number M of links conserved by the inference method.It can be readily shown that these constrains relate the four ROC analysis quantities according to Eq. 4.
Hence, all the information that ROC analysis provides about the inference method is contained in two of the four fundamental variables.Furthermore, as we have shown in the Results section, if we have an estimate on the link density we can condense all the ROC analysis information in just one variable (e.g. the T PR). 13/13

Figure 1 .
Figure 1.Network inference success rates for different networks, coupling strengths, and sizes.Panel (a) [Panel (b)] shows the true positive rate, T PR, as a function of the coupling strength, ε, between N = 131 [N = 277] pulse-coupled Izhikevich maps connected in Erdös-Rényi (ER), Watts-Strogatz (WS), or C. elegans (CE) frontal [global] neural networks;map parameters are set such that the isolated dynamics is bursting (see Methods).The T PR values for the ER and WS come from averaging the results over 10 initial conditions and 20 network realisations with similar topological properties to that of the CE (i.e., number of nodes, average degree, and density of connections).For the CE, the results are averaged only on the initial conditions.The T PR is found by comparing the true underlying network with the binary matrix obtained from cross-correlating the membrane potential time-series (T = 5 × 10 4 iterations) and fixing a threshold such that the inferred density of connections ρ f matches that of the CE, with ρ f 0.08 in panel (a) and ρ f 0.05 in panel (b).The horizontal dashed line in both panels is the random inference T PR, namely, the null hypothesis.

Figure 2 .
Figure 2. Network inference success rates as a function of coupling strength and small-worldness coefficient.Using map and network parameters set as in Fig. 1, panels (a) and (b) show the ensemble and initial-condition averaged T PR as function of ε for N = 131 pulse-coupled Izhikevich maps in Erdös-Rényi (ER) and Watts-Strogatz (WS) network ensembles, respectively.A successive rewiring process 40 is done to each network realisation in order to change its small-worldness coefficient, σ , whilst maintaining the underlying density of connections and degree distribution invariant.The colour code indicates the resultant σ for each rewiring step that increases [panel (a)] or decreases [panel (b)] the network small-worldness.
panels (b) and (c) show the raster plots for ER networks with N = 131 maps coupled with ε = 0.23 and ε = 0.26, respectively.Similarly, Fig. 3(d) shows the same behaviour for the averaged R parameter in networks with ER and 4/13

Figure 3 .
Figure 3. Collective dynamics for different coupling strengths and network structures.In panel (a) we show the ensemble-averaged order parameter, R, for the inter-spike intervals time-series of N = 131 maps connected using the C. elegans frontal neural network (CE, with small-worldness coefficient σ = 2.8), Erdös-Rényi (ER, with σ = 1.0) and Watts-Strogatz (WS, with σ = 2.8) ensembles.For the ER networks, panels (b) (ε = 0.23) and (c) (ε = 0.26) show raster plots indicating the firing pattern of the coupled neuron maps before and after the abrupt drop in panel (a)'s R values.Panel (d) shows R for the rewired ER and WS networks, such that all these ensembles have σ = 2.1 (for comparison, black dots show R for the CE network).Panels (e) (ε = 0.23) and (f) (ε = 0.26) show the corresponding raster plots.

Figure 5 .
Figure 5. Network inference success rates for different networks, coupling strengths, and sizes, using the interspike intervals time-series.Panel (a) [Panel (b)] shows the true positive rate, T PR, as a function of the coupling strength, ε, for N = 131 [N = 277] pulse-coupled Izhikevich maps connected in Erdös-Rényi (ER), Watts-Strogatz (WS), or C. elegans (CE) frontal [global] neural networks; with map parameters set such that the isolated dynamics is bursting (see Methods).The T PR values for the ER and WS come from averaging the results over 10 initial conditions and 20 network realisations with similar topological properties to that of the CE (i.e., number of nodes, average degree, and density of connections).For the CE, the results are averaged only on the initial conditions.The T PR in each case is found by the same process as Fig.1, but using the interspike intevals time-series (considering 2000 spikes) instead of the membrane potentials.The horizontal dashed line in both panels is the random inference T PR, namely, the null hypothesis.

Figure 6 .
Figure 6.Standard deviation of the sensitivity as a function of the coupling strength in the Erdös-Rényi and Watts-Strogatz network ensembles.In panels (a) and (c) [(b) and (d)] we show the standard σ T PR , as a function of the coupling strength, ε, in the N = 131 [N = 277] ER and WS network ensembles, respectively.All network ensembles consist of 20 adjacency matrices, and all points result from an average over 10 different initial conditions.

Figure 7 .
Figure 7. Effect of the time-series length used on the efficiency when inferring Erdös-Rényi networks.Panel (a) [(b)] shows the sensitivity, T PR, as a function of the coupling strength, ε, when inferring N = 131 [N = 277] ER networks using membrane potentials with different time-series length.All points represent an average among 20 different realisations and 10 different initial conditions.