Degree assortativity in networks of spiking neurons

Our nervous system consists of a vast network of interconnected neurons. The network structure is dynamic and connections are formed or removed according to their usage. Much effort has been put into creating a map of all neuronal interconnections; a so-called brain atlas or connectome. Given such a network there are many structural features and measures that one can use to characterise it, e.g. betweenness, centrality, average path-length and clustering coefficient [27].

Obtaining these measures in actual physiological systems is challenging to say the least; nevertheless, insights into intrinsic connectivity preferences of neurons were observed via their growth in culture [8, 35]. Neurons with similar numbers of processes (e.g., synapses and dendrites) tend to establish links with each other – akin to socialites associating in groups and vice-versa. Such an assortativity, typically referred to as a positive assortativity, or a tendency of elements with similar properties to mutually connect, emerges as a strong preference throughout the stages of the cultured neuron development. Furthermore, this preferential attachment between highly-connected neurons is suggested to fortify the neuronal network against disruption or damage [35]. Moreover, a similar positive assortativity is inferred in human central nervous systems as well [8] at both a structural and functional level, where a central “core” in the human cerebral cortex may be the basis for shaping overall brain dynamics [14]. It seems that in no instance, however, is the directional flow of information (e.g., from upstream neuron via axon to synapse and downstream neuron) observed – either in culture or in situ.

Little is known about if and how this positive assortativity or its negative analogue (a tendency for highly connected elements to link with sparsely connected – akin to a “hub-and-spoke” system) influences the dynamics of the neurons networked together. We here consider these effects of assortativity with a network of idealised mathematical neurons known as theta neurons, and represent their connections of downstream (outbound) and upstream (inbound) synapses in a directed graph of nodes and edges. Each neuron then is a node with in- and out-degrees depicting the number of such connections within the network. Assortativity in this context refers to the probability that a neuron with a given in- and out-degree connects to another neuron with a given in- and out-degree. If this probability is what one would expect by chance, given the neurons’ degrees (and this is the case for all pairs of neurons), the network is referred to as neutrally assortative. If the probability is higher (lower) than one would expect by chance — for all pairs — the network is assortative (disassortative). Interchangeably, we will use the term postive assortativity (negative assortativity).

Assortativity has often been studied in undirected networks, where a node simply has a degree, rather than in- and out-degrees (the number of connections to and from a node, respectively) [30, 27, 26]. Since neurons form directed connections, there are four types of assortativity to consider [12]: between either the in- or out-degree of a presynaptic neuron, and either the in- or out-degree of a postsynaptic neuron (Figure 1).

We are aware of only a small number of previous studies in this area [31, 16, 1, 7]. Kähne et al. [16] considered networks with equal in- and out-degrees and investigated degree assortativity, effectively correlating both in- and out-degrees of pre- and post-synaptic neurons. They mostly considered networks with discrete time and a Heaviside firing rate, i.e. a McCulloch-Pitts model [24]. They found that positive assortativity created new fixed points of the model dynamics. Schmeltzer et al. [31] also consider networks with equal in- and out-degrees and investigated degree assortativity. These authors considered leaky integrate-and-fire neurons and derived approximate self-consistency equations governing the steady state neuron firing rates. They found, among other things, that positive assortativity increased the firing rates of high-degree neurons and decreased that of low-degree ones. Positive assortativity also seemed to make the network more capable of sustained activity when the external input to the network was low. De Fransciscis et al. [7] considered assortative mixing of a field of binary neurons, or Hopfield networks. They concluded that assortativity of such simple model neurons exhibited associative memory (similar to bit fields of a magnetic storage medium), and robustly so in the presence of noise that negatively assortative networks failed to withstand. Avalos-Gaytan et al. [1] considered the effects of dynamic weightings between Kuramoto oscillators — effectively a dynamically evolving network — on assortativity. They observed that if the strength of connections between oscillators increased when they were synchronised, a strong positive assortativity evolved in the network, suggesting a potential mechanism for the creation of assortative networks, as observed in cultured neurons mentioned above, and as we study here.

To briefly summarise our results, we find that only two out of the four types of degree assorativity have any influence on the network’s dynamics: those when the in-degree of a presynaptic neuron is correlated with either the in- or out-degree of a postsynaptic neuron. Of these two, (in,in)-assortativity has a greater effect than (in,out)-assortativity. For both cases, negative assortativity widens the parameter range for which the network is bistable (for excitatory coupling) or undergoes oscillations in mean firing rate (for inhibitory coupling), and positive assortativity has the opposite effect.

Our work is similar in some respects to that of Restrepo and Ott [30] who considered degree assortativity in a network of Kuramoto-type phase oscillators. They found that for positive assortativity, as the strength of connections between oscillators was increased the network could undergo bifurcations leading to oscillations in the order parameter, in contrast to the usual scenario that occurs for no assortativity. However, their network was undirected, and thus there is only one type of degree assortativity possible.

The outline of the paper is as follows. In Sec. 2 we present the model and then derive several approximate descriptions of its dynamics. In Sec. 3 we describe the method for creating networks with prescribed types of degree assortativity, and in Sec. 4 we discuss aspects of the numerical implementation of the reduced model. Results are given in Sec. 5 and we conclude with a discussion in Sec. 6. Appendix A contains the algorithms we use to generate networks with prescribed assortativity.

We consider a network of $N$ theta neurons:

for $j=1,2,\dots N$ where

$\eta_{j}$ is a constant current entering the $j$th neuron, randomly chosen from a distribution $g(\eta)$, $K$ is strength of coupling, $\langle k\rangle$ is mean degree of the network, and the connectivity of the network is given by the adjacency matrix $A$, where $A_{jn}=1$ if neuron $n$ connects to neuron $j$, and zero otherwise. The connections within the network are either all excitatory (if $K>0$) or inhibitory (if $K<0$). Thus we do not consider the more realistic and general case of a connected population of both excitatory and inhibitory neurons, although it would be possible using the framework below.

The theta neuron is the normal form of a Type I neuron which undergoes a saddle-node on an invariant circle bifurcation (SNIC) as the input current is increased through zero [10, 11]. A neuron is said to fire when $\theta$ increases through $\pi$, and the function

in (2) is meant to mimic the current pulse injected from neuron $n$ to any postsynaptic neurons when neuron $n$ fires. $a_{q}$ is a normalisation constant such that $\int_{0}^{2\pi}P_{q}(\theta)d\theta=2\pi$ independent of $q$.

The in-degree of neuron $j$ is defined as the number of neurons which connect to it, i.e.

while the out-degree of neuron $n$ is the number of neurons it connects to, i.e.

Since each edge connects two neurons, we can define the mean degree

Networks such as (1) have been studied by others [22, 3, 17], and note that under the transformation $V=\tan{(\theta/2)}$ the theta neuron becomes the quadratic integrate-and-fire neuron with infinite thresholds [25, 21].

We are interested in the effects of degree assortativity on the dynamics of the network of neurons. We will choose a default network with no assortativity and then introduce one of the four types of assortativity and investigate the changes in the network’s dynamics. Our default network is of size $N=5000$ neurons where in- and out-degrees $k$ for each neuron are independently drawn from the interval $[750,2000]$ with probability $P(k)\sim k^{-3}$ (i.e. a power law, as found in [9] and used in [3]). We create networks using the configuration model [27], then modify them using algorithms which introduce assortativity and then remove multiple connections between nodes (or multi-edges) (described in Appendix A). Further, to observe any influence of multi-edges on the dynamics investigated here, we also developed a novel network assembly technique permitting introduction of very high densities of multi-edges also described in the Appendix; we refer to this novel technique as the “permutation” method. We choose as our default parameters $\eta_{0}=-2,\Delta=0.1,K=3$, for which a default network approaches a stable fixed point. The sharpness of the synaptic pulse function is set to $q=2$ for all simulations.

We first check the validity of averaging over an infinite ensemble. We assemble 20 different default networks and for each, run (19)-(20) to a steady state and calculate the order parameter $z$, the mean of $B\bm{b}$. The real part of $z$ is plotted in orange in Fig. 2. For each of these networks we then generated 50 realisations of the $\eta_{i}$’s and ran (1)-(2) for long enough that transients had decayed, and then measured the corresponding order parameter for the network of individual neurons

and plotted its real part in blue in Fig. 2. Note that the orange circles always lie well within the range of values shown in blue. The fact that deviations within the 50 realisations are small relative to the value obtained by averaging over an infinite ensemble provide evidence for the validity of this approach, at least for these parameter values.

We also investigate the influence of multi-connections (i.e. more than one connection) between neurons on the network dynamics. The configuration model creates a network in which the neuron degrees are exactly those specified by the choice from the appropriate distribution, but typically results in both self-connections and multiple connections between neurons. We have an algorithm for systematically removing such connections while preserving the degrees, and found that removing such edges has no significant effect (results not shown). We also have an algorithm (see “Permutation Method” in Appendix A) for increasing the number of multi-edges from the number found using the default configuration model. This novel network assembly method meets specified neuron degrees and also produces specified densities of multi-connections ranging from none to 97%; see Fig. 3 for the results of such calculations. We see that only when the fraction of multi-edges approaches $90\%$ do we see a significant effect. However, in our simulations we use simple graphs without multi-edges.

For networks of the size we investigate it is impractical to consider each distinct in- and out-degree (because $E$ will be very large and sparse). Due to the smoothness of the degree dependency of $b(\bm{k})$ we coarse-grain in degree space by introducing “degree clusters” — lumping all nodes with a range of degrees into a group with dynamics described by a single variable. Let there be $N_{c^{\text{in}}}$ clusters in in-degree and $N_{c^{\text{out}}}$ clusters in out-degree, with a total of $N_{\bm{c}}=N_{c^{\text{in}}}\cdot N_{c^{\text{out}}}$ degree clusters. The matrix $C$ then is an $N_{\bm{c}}\times N$ matrix and constructed as previously, except that $d(j)$ is not the degree index of neuron $j$, but the cluster index and $s$ is the cluster index running from 1 to $N_{\bm{c}}$. Similarly for the matrix $B$. There are multiple options for how to combine degrees into a cluster. The cluster index of a neuron can be computed linearly, corresponding to clusters of equal size in degree space. However, with this approach, depending on the degree distribution, some of the clusters may be empty or hardly filled, resulting in poor statistics. To overcome this issue, the cumulative sum of in- and out-degree distribution can be used to map degrees to cluster indices. Thus, clusters are more evenly filled and at the same time regions of degree space with high degree probability are more finely sampled. The dynamical equations (30)-(31) are equally valid for describing degree cluster dynamics with $s,t\in[1,N_{\bm{c}}]$ and $E=CAB$, where $C$ and $B$ are cluster versions of their previous definitions.

To check the effect of varying the number of clusters we generate 20 default matrices and then generate the corresponding matrix $E$ with varying numbers of clusters ($N_{c^{\text{in}}}$ and $N_{c^{\text{out}}}$ are equal), then run (30)-(31) to a steady state and plot the real part of $z$ in Fig. 4. We see that the order parameter is well approximated using as little as about $N_{c^{\text{in}}}=N_{c^{\text{out}}}=10$ degree clusters. Beyond that, fluctuations between different network realisations exceed the error introduced by clustering. In our simulations we stick to the choice of 10 degree clusters per in- and out-degree space.

Having performed this clustering, we find that it is possible to represent $E$ using a low-rank approximation, calculated using singular value decomposition. Thus for a fixed $r$ we have

where $S$ is a diagonal matrix with decreasing entries, called singular values, and $U$ and $V$ are unitary matrices. In Fig. 5 we plot the largest 6 singular values of $E$ as function the assortativity coefficient, for the 4 types of assortativity. Even for large $|r|$ the singular values decay very quickly, thus a low-rank approximation is possible. We choose a rank-3 approximation, so approximate $E$ by

where $u_{i}$ is the $i$th column of $U$, similarly for $v_{i}$ and $V$, and $s_{i}$ is the $i$th singular value. We have such a decomposition at discrete values of $r$ and use cubic interpolation to evaluate $E(r)$ for any $r$. This decomposition means that the multiplication in (31) can be evaluated quickly using 3 columns of $U$ and $V$ rather than the full $N_{\bm{c}}\times N_{\bm{c}}$ matrix $E$.

We note that the components for the approximation of $E(r)$ are calculated once and then stored, making it very easy to systematically investigate the effects of varying any of the parameters $\eta_{0},\Delta,K$ and $q$ (governing the sharpness of the pulse function (3)).

We investigated the effects of degree assortativity on the dynamics of a network of theta neurons. We used the Ott/Antonsen ansatz to derive evolution equations for an order parameter associated with each neuron, and then coarse-grained by degree and then degree cluster, obtaining a relatively small number of coupled ODEs, whose dynamics as parameters varied could be investigated using numerical continuation. We found that degree assortativity involving the out-degree of the sending neuron, i.e. (out,in) and (out,out), has no effect on the networks’ dynamics. Further, (in,out) assortativity moves bifurcations slightly, but does not lead to substantial differences in dynamical behaviour. The most significant effects were caused by creating correlation between in-degrees of the sending and receiving neurons. For our excitatorially coupled example, positive (in,in) assortativity narrows the bistable region, whereas negative assortativity widens it (see Fig. 7). In the inhibitory case introducing negative assortativity increased the amplitude of network oscillations and extended their range to slightly larger $\eta_{0}$. On the contrary, positive (in,in) assortativity in this network has an opposite effect and eventually stops oscillations (see Fig. 9).

The most similar work to ours is that of [3]. These authors also considered a network of the form (1)-(2) and by assuming that the dynamics depend on only a neuron’s degree and that the $\eta_{j}$ are chosen from a Lorentzian, and using the Ott/Antonsen ansatz, they derived equations similar to (30)-(31). The difference in formulations is that rather than a sum over entries of $E$ (in (31)), [3] wrote the sum as

where $P({\bf k})$ is the degree distribution and $a({\bf k}^{\prime}\rightarrow{\bf k})$ is the assortativity function, which specifies the probability of a link from a node with degree ${\bf k}^{\prime}$ to one with degree ${\bf k}$ (given that such neurons exist). They then chose a particular functional form for $a$ and briefly presented the results of varying one type of assortativity (between $k_{in}^{\prime}$ and $k_{out}$). In contrast, our approach is far more general (since any connectivity matrix $A$ can be reduced to the corresponding $E$, the only assumption being that the dynamics are determined by a neuron’s degree). We also show the results of a wider investigation into the effects of assortativity.

This alternative presentation also explains why $E$ can be well approximated with a low-rank approximation. If the in- and out-degrees of a single neuron are independent, $P({\bf k}^{\prime})=P_{i}(k_{in}^{\prime})P_{o}(k_{out}^{\prime})$, and with neutral assortativity, $a({\bf k}^{\prime}\rightarrow{\bf k})=k_{out}^{\prime}k_{in}/(N\langle k\rangle)$. Thus

This term contributes to the input current to a neuron with degree ${\bf k}=(k_{in},k_{out})$, but is independent of $k_{out}$. Thus the state of a neuron depend only on its in-degree, so

Comparing with (31) we see that $E={\bf c}^{T}{\bf d}$ where ${\bf c}=(k_{in}^{1},k_{in}^{2}\dots k_{in}^{N_{k^{in}}})/N$ and ${\bf d}=(P_{i}(k_{in}^{1}),P_{i}(k_{in}^{2})\dots,P_{i}(k_{in}^{N_{k^{in}}}))$, i.e. $E$ is a rank-one matrix. Varying assortativity within the network is then a perturbation away from this, with the effects appearing in the second (and third) singular values in the SVD decomposition of $E$.

A limitation of our study is that we considered only networks of fixed size with the same distributions of in- and out-degrees, and a specific distribution of these degrees. However, our approach does not rely on this and could easily be adapted to consider other types of networks, although we expect it to become less valid as both the average degree and number of neurons in the network decrease. We have also only considered theta neurons, but since a theta neuron is the normal form of a Type I neuron, we expect similar networks of other Type I neurons to behave similarly to the networks considered here. The approach presented here could also be used to efficiently investigate the effects of correlated heterogeneity, where either the mean or width of the distribution of the $\eta_{j}$ is correlated with a neuron’s in- or out-degree [33, 34, 6]. We could also consider assortativity based on a neuron’s intrinsic drive ($\eta_{j}$) [32] rather than its degrees, or correlations between an individual neuron’s in- and out-degree [36, 20, 23, 37, 28]. We are currently investigating such ideas.

Acknowledgements: This work is partially supported by the Marsden Fund Council from Government funding, managed by Royal Society Te Apārangi. We thank the referees for their helpful comments.