Degree assortativity in networks of spiking neurons

Degree assortativity refers to the increased or decreased probability of connecting two neurons based on their in- or out-degrees, relative to what would be expected by chance. We investigate the effects of such assortativity in a network of theta neurons. The Ott/Antonsen ansatz is used to derive equations for the expected state of each neuron, and these equations are then coarse-grained in degree space. We generate families of effective connectivity matrices parametrised by assortativity coefficient and use SVD decompositions of these to efficiently perform numerical bifurcation analysis of the coarse-grained equations. We find that of the four possible types of degree assortativity, two have no effect on the networks' dynamics, while the other two can have a significant effect.


Introduction
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 pathlength 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 "huband-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 pairsthe 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 postsynaptic 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.

Model description and simplifications
We consider a network of N theta neurons: η j is a constant current entering the jth neuron, randomly chosen from a distribution g(η), K is strength of coupling, k 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 saddlenode on an invariant circle bifurcation (SNIC) as the input current is increased through zero [10,11]. A neuron is said to fire when θ increases through π, and the function (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 2π 0 P q (θ)dθ = 2π independent of q.
The in-degree of neuron j is defined as the number of neurons which connect to it, i.e.
A jn 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 (θ/2) the theta neuron becomes the quadratic integrate-and-fire neuron with infinite thresholds [25,21].

2.
1. An infinite ensemble. As a first step we consider an infinite ensemble of networks with the same connectivity, i.e. the same A jn , but in each member of the ensemble, the value of η j associated with the jth neuron is randomly chosen from the distribution g(η) [2]. Thus we expect a randomly chosen member of the ensemble to have values of η in the ranges with probability g(η 1 )g(η 2 ) . . . g(η N )dη 1 dη 2 . . . dη N . The state of this member of the ensemble is described by the probability density which satisfies the continuity equation If we define the marginal distribution for the jth neuron as where we have now evaluated I j as an average over the ensemble rather than from a single realisation (as in (2)). This is reasonable in the limit of large networks [2]. Multiplying (9) by k =j dθ k dη k and integrating we obtain A network of theta neurons is known to be amenable to the use of the Ott/Antonsen ansatz [29,22,17] so we write The dependence on θ j is written as a Fourier series where the kth coefficient is the kth power of a function α j . Substituting this into (12) and (11) we find that α j satisfies where an overbar indicates complex conjugate, and Assuming that g is a Lorentzian: we can use contour integration to evaluate the integral in (15), and evaluating (14) at the appropriate pole of g we obtain and z j = e iθ j , where the expected value is taken over the ensemble. Now (19) is a set of N coupled complex ODEs, so we have not simplified the original network (1) in the sense of decreasing the number of equations to solve. However, the states of interest are often fixed points of (19) (but not of (1)), and can thus be found and followed as parameters are varied. At this point the network we consider, with connectivity given by A, is arbitrary. If A was a circulant matrix, for example, this would represent a network of neurons on a circle, where the strength of coupling between neurons depends only on the distance between them [17].

2.2.
Lumping by degree. The next step is to assume that for a large enough network, the dynamics of neurons with the same degrees will behave similarly [3]. Such an assumption has been made a number of times in the past [16,15,30]. We thus associate with each neuron the degree vector k = (k in , k out ) and assume that the value of z for all neurons with a given k are similar. There are N k = N k in N k out distinct degrees where N k in and N k out are the number of distinct in-and out-degrees, respectively. We define b s to be the order parameter for neurons with degree k s , where s ∈ [1, N k ], and now derive equations for the evolution of the b s .
Let z be the vector of ensemble states z j , where j ∈ [1, N ] and the degree index of neuron j be d(j), such that k d(j) is its degree. We assume that for all neurons with the same degree k d(j) = k s the ensemble state z j is similar in sufficiently large networks and thus we only care about the mean value z j d(j)=s = b s with s ∈ [1, N k ]. We say that degree k s occurs h s times and thus write where the N k × N matrix C has h s entries in row s, each of value 1/h s , at positions j where d(j) = s and zeros elsewhere, i.e. C sj = δ s,d(j) /h s with δ being the Kronecker delta.
To find the time derivative of b we need to express z in terms of b, which we do with an N × N k matrix B which assigns to z j the corresponding b s value, such that with components B js = δ d(j),s . Note that CB = I N k , the N k × N k identity matrix. Differentiating (21) with respect to time, inserting (19) into this and writing z in terms of b using (22) Considering that for all t there is only a single non-zero entry B jt , equal to 1, the identity holds for any power n. Further we find that Thus, the local term in (23) For the non-local term we writė whereJ s describes the synaptic current of the ensemble equations averaged over nodes sharing the same degree k s . The identity (24) also applies to (20), so that and the current can be written as The effective connectivity between neurons with different degrees is therefore expressed in the matrix E = CAB and we end up with equations governing the b s : These equations are of the same form as (19)- (20) except that A has been replaced by E. Note that the connectivity matrix A is completely general; we have only assumed that neurons with the same degrees behave in the same way. We are not aware of a derivation of this form being previously presented.

Network assembly
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) ∼ 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 η 0 = −2, ∆ = 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 Bb. The real part of z is plotted in orange in Fig. 2. For 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18  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.
3.1. Assortativity. For a given matrix A we can measure its assortativity by calculating the four Pearson correlation coefficients r(α, β) with α, β ∈ [in, out] which read N e being the number of connections and the leading superscript s or r refers to the sending or receiving neuron of the respective edge. For example the sending node's in-degree of the second edge would be s k in 2 . Note that there are four mean values to compute.
We introduce assortativity by randomly choosing two edges and swapping postsynaptic neurons when doing so would increase the target assortativity coefficient [31]. An edge (i, j) is directed from neuron j to neuron i. In order to know whether the pair (i, j) and (h, l) should be rewired or left untouched, we compare their contribution to the covariance in the numerator of (33): If c \ / > c we replace the edges (i, j) and (h, l) by (i, l) and (h, j), respectively, otherwise we do not, and continue by randomly choosing another pair of edges. Algorithm 1 (see Appendix A) demonstrates a scheme for reaching a certain target assortativity coefficient.
We investigate the effects of different types of assortativity (see Fig 1) in isolation. We thus need a family of networks parametrised by the relevant assortativity coefficient. Algorithm 1 is used to create a network with a specific value of one of the assortativity coefficients, but especially for high values of assortativity it may be that in doing so a small amount of assortativity of a type other than the intended one is introduced. Accordingly, it may be necessary to examine all types of assortativity and apply the mixing scheme to reduce other types back to zero, and then (if necessary) push the relevant value of assorativity back to its target value. We do multiple iterations of these mixing rounds until all assortativities are at their target values (which may be 0) within a range of ±0.005. We use Algorithm 1 with a range of target assortativities r, and for each value, store the connectivity matrix A and thus form the parametrised family E(r). We do this for the four types of assortativity.
We have chosen to use the configuration model to create networks with given degree sequences and then introduced assortativity by swapping edges. Although we developed our novel "permutation" method as well (see Appendix A), that method was designed for assembling adjacency matrices with desired multi-edge densities and was applied only for that aspect. By contrast, another common adjacency network assembly technique, that of Chung and Lu [5] together with an analytical expression for assortativity (as in [3]), proved inadequate. We found that the latter approach significantly alters the degree distribution for large assortativity, whereas the configuration model combined with our mixing algorithm does not change degrees at all. For our default network this approach allows us to introduce assortativity of any one kind up to r = ±0.5.

Implementation
For networks of the size we investigate it is impractical to consider each distinct inand out-degree (because E will be very large and sparse). Due to the smoothness of the degree dependency of b(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 in clusters in in-degree and N c out clusters in out-degree, with a total of N c = N c in · N c out degree clusters. The matrix C then is an N c × 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 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 ∈ [1, N 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 in and N c 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 in = N c 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 ith column of U , similarly for v i and V , and s i is the ith 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 c × N 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 η 0 , ∆, K and q (governing the sharpness of the pulse function (3)).

Excitatory coupling.
We take K = 3 to model a network with only excitatory connections. To study the dynamical effect of assortativity we generate positive and negative (r = ±0.2) assortative networks of the four possible kinds and follow fixed points of (30)-(31) as a function of η 0 , and compare results with a neutral (r = 0) network. We use pseudo-arc-length continuation [18,13].
To calculate the mean frequency over the network we evaluate z = Bb and then use the result that if the order parameter at a node is z, then the frequency of neurons at that node is [19,25] (39) 1 π Re 1 −z 1 +z Averaging these gives the mean frequency.
Results are shown in Figure 6, where we see quite similar behaviour in each case: apart from a bistable region containing two stable and one unstable fixed point, there is only a single stable fixed point present. Further, the two assortativity types (out,in) and (out,out) apparently do not affect the dynamics, whereas the saddle-node bifurcations marking the edges of the bistable region move slightly for (in,out) and significantly for (in,in) assortativity. Following the saddle-node bifurcations for the latter two cases we find the results shown in Figure 7. We have performed similar calculations for different networks with the same values of assortativity and found similar results.

5.2.
Inhibitory coupling. We choose K = −3 to model a network with only inhibitory coupling. Again, we numerically continue fixed points for zero, positive and negative assortativity (r = 0, ±0.2) as η 0 is varied and obtain the curves shown in Figure 8. Consider the lower left plot. For large η 0 the system has a single stable fixed point which undergoes a supercritical Hopf bifurcation as η 0 is decreased, creating a stable periodic orbit. This periodic orbit is destroyed in a saddle-node bifurcation on an invariant circle (SNIC) bifurcation at lower η 0 , forcing the oscillations to stop. Decreasing η 0 further, two unstable fixed points are destroyed in a saddle-node bifurcation. In contrast with the case of excitatory coupling, oscillations in the average firing rate are seen. These can be thought of as partial synchrony, since some fraction of neurons in the network have the same period and fire at similar times to cause this behaviour. The period of this macroscopic oscillation tends to infinity as the SNIC bifurcation is approached, as shown in the inset of the lower left panel in Fig. 8.
As in the excitatory case, we see that assortativities of type (out,in) and (out,out) have no influence on the dynamics in this scenario. However, type (in,out) does have a small effect, slightly moving bifurcation points (top right panel in Fig. 8). Type (in,in) has the strongest effect, resulting in a qualitative change in the bifurcation scenario for large enough assortativity: there is a region of bistability between either two fixed points or a fixed point and a periodic orbit. This is best understood by following the bifurcations in the top panels of Fig. 8 as r is varied, as shown in Figure 9. There is one fixed point in regions A, B and D, and three in region C. For (in,out) assortativity there is a stable periodic orbit in region B and never any bistability.
We now describe the case for (in,in) assortativity. For negative and zero r the scenario is the same as for the other three types, but as r is increased there is a Takens-Bogdanov bifurcation where regions C,D,E and F meet, leading to the creation of a curve of homoclinic bifurcations, which is destroyed at another codimension-two point where there is a homoclinic connection to a non-hyperbolic fixed point [4]. There are stable oscillations in region E, created or destroyed in supercritical Hopf or homoclinic bifurcations. In region F there is bistability between two fixed points.

Discussion
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 indegrees 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 η 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 η 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 (k) is the degree distribution and a(k → k) is the assortativity function, which specifies the probability of a link from a node with degree k to one with degree 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 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, This term contributes to the input current to a neuron with degree 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 = c T d where c = (k 1 in , k 2 in . . . k N k in in )/N and d = (P i (k 1 in ), P i (k 2 in ) . . . , P i (k N k in 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 η 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 (η 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.
Algorithm 1: Assortative mixing. Randomly pair up all N e edges of the network with adjacency matrix A and reconnect them at once where preferable with respect to target assortativity r target . Repeat the process until the assortativity coefficient lies within the tolerance. Once overshooting the target coefficient, interpolate the length of a shortened list of edge pairs and reconnect those. showing arrangement of edge entries (all solo connections: ρ + m = 0) for each row aligned left where row sums add up to k in . If multi-edges are desired, we simply distribute them in the rows of A (0) satisfying the proportion, ρ + m = 0 and the row sum. Bottom: permutation of rows in A (0) into this example A (1) . Note row sums still add up to k in , with column sums adding to a current k (1) out -likely violating the designated k out .
Algorithm 3: Permutation Method, Phase (2). Pseudo-code for permuting initial matrix A (0) into final matrix A (n) satisfying all constraints, k in , ρ + m and k out .