Synchronization from Second Order Network Connectivity Statistics

We investigate how network structure can influence the tendency for a neuronal network to synchronize, or its synchronizability, independent of the dynamical model for each neuron. The synchrony analysis takes advantage of the framework of second order networks, which defines four second order connectivity statistics based on the relative frequency of two-connection network motifs. The analysis identifies two of these statistics, convergent connections, and chain connections, as highly influencing the synchrony. Simulations verify that synchrony decreases with the frequency of convergent connections and increases with the frequency of chain connections. These trends persist with simulations of multiple models for the neuron dynamics and for different types of networks. Surprisingly, divergent connections, which determine the fraction of shared inputs, do not strongly influence the synchrony. The critical role of chains, rather than divergent connections, in influencing synchrony can be explained by their increasing the effective coupling strength. The decrease of synchrony with convergent connections is primarily due to the resulting heterogeneity in firing rates.

results from specifying only the connection probabilities, as it adds minimal additional structure beyond what is required to by the connection probabilities, i.e., the independent model maximizes entropy (Jaynes, 1957) subject to constraints set by the connection probabilities.
Evidence is emerging that structure in neuronal networks is not captured by the independent model. Some studies (Holmgren et al., 2003;Song et al., 2005;Oswald et al., 2009) but not all (Lefort et al., 2009) found reciprocal connections between pairs of neurons to be more likely than predicted by the independent model. Studies have also found that the frequencies of other patterns of connections differ from the independent model prediction (Song et al., 2005;Perin et al., 2011).
The independent random network model is based only on the first order statistics of the network connectivity, i.e., the marginal probability that any given connection exists. In the independent model, higher order statistics, such as the covariances among connections, are set to zero. We have developed a model framework extending the independent model to allow non-zero second order statistics, which determine the covariance among connections (Zhao et al., in preparation). We refer to these network models as SONETs. SONETs retain the principle from the independent model of adding minimal structure beyond what is required to match the specified constraints (in this case, matching both the first order connection probabilities and the second order connectivity statistics).
The first class of SONETs that we have developed are for the special case where the connection probability among neurons is fixed to a single constant p (Zhao et al., in preparation). In this way, these SONETs are a generalization of the Erdős-Rényi random network model. Although the second order connectivity statistics are specified as covariances among connections, it is equivalent, but more intuitive, to determine them using the probability of observing patterns of two connections (or motifs). There are four possible motifs defining SONETs: reciprocal, convergent, divergent, and chains, as illustrated in Figure 1. A given SONET is determined by specifying the connection probability p and the probability of each of these motifs.
To define the SONETs and the second order connectivity statistics more precisely, we introduce notation for the connectivity matrix W which defines the connectivity among N neurons. The component W ij = 1 indicates a connection from neuron j onto neuron i, and W ij = 0 if that connection is absent. In SONETs self coupling is not allowed; therefore W ii = 0.
In this paper we consider only homogeneous networks where the probability of any combination of connections is independent of the neurons' labels. Thus, there is only one first order connectivity statistic: the connection probability p that determines the marginal probability of the existence of any connection W ij , (1) For the independent random network model the probability of any pair of connections would be p 2 . We introduce second order connectivity statistics that measure deviations from independence among pairs of connections that share a neuron. Since the network is homogeneous, we define a single second order connectivity statistic for each motif of Figure 1. The statistics a recip , a conv , a div , and a chain specify how the probabilities of the reciprocal, convergent, divergence, and chain motifs, respectively, deviate from independence: for any triplet of distinct neurons i, j, and k. In this way, a non-zero a yields a correlation between the connections of its corresponding motif. The SONET model is a five-dimensional probability distribution (parametrized by p and the four a's) for the network connectivity W that satisfies conditions (1) and (2). As the probability distribution adds minimal structure (i.e., higher order correlations) beyond these constraints, the SONET model reduces to the Erdős-Rényi model when all a's are zero.
For technical reasons, the SONET model only approximately satisfies the maximum entropy condition for minimal structure (Jaynes, 1957), but instead uses dichotomized Gaussian random variables analogous to those of Macke et al. (2009). The dichotomized Gaussian random variables satisfy the maximum entropy condition for continuous random variables, which are then dichotomized to be either zero or one, forming the W ij . Hence, to generate a SONET, one simply samples a N(N − 1)-dimensional vector of joint Gaussian random variables and thresholds them into the Bernoulli random variables W ij . The main technical challenge is calculating the covariance matrix of the joint Gaussian so that W satisfies the second order statistics (2). See (Zhao et al., in preparation) for details.

estiMating connectivity statistics
The above definitions of the first and second order connectivity statistics are based on a probability distribution for the connectivity matrix W. From a given network with connectivity matrix W, one can estimate these statistics directly from W. We use p and â to denote these descriptive statistics.
The average connectivity p is simply the ratio of the number of connections N conn present in W with the total possible: ˆ/ ( ) p N N N = − conn 1 , where N conn = ||W|| 1 ≡Σ i,j |W ij |. The average second order connectivity statistics can likewise be determined by the total number of corresponding two-connection motifs that are present in the network, divided by the total possible:  Hence, a conv and a div are nearly equal to the variance of the in-degree and out-degree distributions, normalized by the mean degree squared (with a small offset 1/E(d)). a chain is nearly equal to the covariance between in-degree and out-degree, normalized by the mean degree squared.

generating networks with correlated power law degree distribution
We developed an algorithm based on the approach by Chung and Lu (2002) to generate directed networks with expected power law degree distributions with a range of second order connectivity statistics (3). We let p d (x, y) be this degree distribution, with x = d in and y = d out . We set the marginal distributions of the in-and outdegrees to truncated power laws with exponentb i : where the cutoff parameters satisfy 0 < a i < b i < N, and the normalization constants c i and d i are chosen to make the marginal distributions continuous at a i and integrate to 1. Below the lower cutoff a i , we let the distribution increase smoothly with exponent g i . For a given set of networks, we fixed the right cutoff parameters b i and the rising exponents g i . The left cutoff parameters a i and the exponent b i were then chosen for each network to match the desired average connectivity p as well as the second order connectivity statistics a conv (for i = in) or a div (for i = out) by numerically solving the equations [cf.
for (a in , b in ) and analogous equations for (a out , b out ) in terms of a div .
To vary the correlation between in-and out-degrees and thus a chain , we let the full degree distribution p d (x,y) be determined by the marginal distributions p in and p out and a Gaussian copula with correlation coefficient r. In this way, we can sample from p d (x,y) by a two step process. We generate a joint Gaussian with correlation coefficient r and compose each component with the cumulative where N k for k ∈ {recip, conv, div, chain} is the number of times motif k appears in the network. The N k can be quickly calculated from W via N recip = Tr(W 2 )/2, N conv = (||W T W|| 1 −||W|| 1 )/2, N div = (||WW T || 1 −||W|| 1 )/2, N chain = (||W 2 || 1 −Tr(W 2 ), where W T denotes matrix transpose and Tr(W) denotes trace.
In our network simulations, we compare the measured α to the network synchrony even when we know the underlying probability distribution and hence the parameters α from the definition (2). Because these networks are generated probabilistically, the actual α will vary slightly from α. We find that α predicts synchrony better than α.

relating second order connectivity statistics to degree distribution
A common method to describe network structure is through the degree distribution (Newman, 2003;Boccaletti et al., 2006). In this section we derive a relationship between the second order connectivity statistics and degree distribution. This relation will be useful in understanding the effect of the connectivity statistics and synchrony.
The in-degree of a neuron is the total number of connections directed toward a neuron, d W i j i j in = ∑ for neuron i, and the outdegree is the total number of connections directed away from the neuron, d W i j j i out = ∑ . The mean in-degree is equal to the mean out-degree, so we just refer to a mean degree.
Using the definitions (2) of a, one can calculate that the variance of the in-degree depends on a conv , the variance of the out-degree depends on a div , and the covariance between in-degree and outdegree depends on a recip and a chain : where i is an arbitrary neuron index. We can simplify the relationships between the a and the second order statistics of the degree distribution in the limit of large network size N. Neglecting terms that go to zero in this limit 1 and rewriting in terms of the expected mean degree E(d) = (N − 1)p, we find the following relationships between the connectivity statistics and the degree distribution: However, a tighter spread of the eigenvalues will make it easier to scale them by S to fit within the critical region. Therefore, changes in network structure that decrease the spread of the eigenvalues will favor synchronizability. One measure of this eigenvalue spread for asymmetric matrices, as suggested by Nishikawa and Motter (2010), is the variance s m 2 of the eigenvalues normalized by the mean degree squared, where m m = ∑ − = 1 1 2 N i N i is the mean of the eigenvalues and d is the mean degree (4). This variance ignores the zero eigenvalue m 1 = 0 present in every Laplacian matrix.

Modifications for chemical synapses
The master stability function requires the subtraction of H(x i ) in (8) to ensure that the synchronous state x i t t ( ) ( ) = ∧ x exists and satisfies the equation d . The H(x i ) term is suitable for neuron models with electrical synapses but is not an appropriate model for chemical synapses because it implies that a connection from neuron j onto neuron i leads to an equivalent negative effect from neuron i onto itself. To model chemical synapses, we must remove the offending H(x i ) term, yielding: The rows of the connectivity matrix W may not sum to the same value, in contrast to the Laplacian matrix L of (9) which must have rows that sum to zero. Hence, if we substitute the synchronous state x x i =ˆ, we would get different, inconsistent equations that x must satisfy: where the ith row sum of W is d i in , the incoming degree of neuron i. If all neurons do not have the same incoming degree, then we cannot solve for x independent of i and the synchronous state does not even exist.
Therefore, for chemical synaptic networks (11), rather than analyzing the stability of the synchronous state, we seek to determine how far the network is from having a synchronous solution. In the spirit of the master stability function analysis, we desire a measure of synchronizability analogous to (10) based only on network properties. We propose that network synchrony will decrease with the variance of the in-degrees, normalized by their mean squared

Connecting stability measures to network topology
Fortunately, we can merge the two normalized variances s m 2 and s d 2 into the same measure. It turns out that, for many large networks, the eigenvalues of the Laplacian matrix L = D − W are dominated by the diagonal entries D, which are the in-degrees d i in . If the connectivity matrix were normal (WW T = W T W), then this distribution function of the Gaussian to obtain a sample (u, v) from the Gaussian copula. Next, we compose the result with the inverse of the cumulative distribution functions F to obtain the expected in-degree and out-degree ( , ) ( ( ), ( )) x y F u F v = − − in out 1 1 as a sample from p d (x,y).
Finally, to generate the network, we sample N points (x i ,y i ) from p d (x,y), which give the expected in-degree and out-degree, respectively, of each neuron i. We let the probability of each network connection W ij from neuron j to neuron i be proportional to x i ,y j , where the proportionality constant is chosen so that the average connectivity is Pr(W ij = 1)=p. As described in (Chung and Lu, 2002), as long as max i,j x i y j ≤N(N−1)p, this approach yields a network whose distribution of expected degrees is given by p d (x,y). This condition was not satisfied for some networks so that their degree distribution deviated from the prescribed p d (x,y). Since we measure the second order statistics directly from the generated connectivity matrix via (3), such deviations were not a concern.
analysis of coMplete synchrony To examine the influence of the second order connectivity statistics on network synchrony, we will analyze the stability of the two extreme cases of synchrony, the perfectly synchronous state and the completely asynchronous state. Through linearizations around these states, we derive conditions between the connectivity statistics and synchrony that are valid near these states. We begin by analyzing complete synchrony.

Stability of synchrony determined using master stability function
To determine the stability of the perfectly synchronous case, we linearize the equations for network dynamics around synchrony using the master stability function approach by Pecora and Carroll (1998). This approach is based on determining synchrony in a system of identical, noise-free, coupled oscillators: where the state of oscillator i at time t is described by the vector x i (t). The dynamics of this system are determined by the scalar coupling strength S and two vector-valued functions of the same dimension as the interval variable vector x i : the internal dynamics function F and the coupling function H. We rewrite (8) in terms of the Laplacian matrix L = D − W, where D is the diagonal of row sums of the connectivity matrix W, Using the master stability function approach, one can determine a critical stability region in the complex plane that is based solely on the neuron dynamics (given by F and H) and is independent of network connectivity. The influence of the network on the stability enters through the eigenvalues m i of L, as the eigenvalues scaled by the coupling strength Sm i must lie in the critical region for the synchronous state to be stable (Arenas et al., 2008).
If the neuron model is not known, it is not possible to know the critical region in the complex plane, reflecting the fact that the topology alone cannot determine whether a network will synchronize. stability of the asynchronous state

The emergence of synchrony
An approach to perform a similar stability analysis for the completely asynchronous state of heterogeneous, noisy oscillators has been recently developed by Restrepo et al. (2006a,b). The asynchronous state is where the oscillators are as spread out evenly as though they were uncoupled (in a sense determined by the oscillator dynamics). The analysis begins with a coupled oscillator model similar to the one used in the master stability function approach (8): where j i (t) is a noise term and the model is specified by the functions F, H, and k. One notable difference in this model is that the oscillators can be heterogeneous since their dynamics depend on the parameter vectors h i . Two important conditions of the analysis are that the parameters cannot be correlated with the network structure and that each neuron must receive many inputs. For the asynchronous state to be an approximate solution of this system, the factor multiplying W must be approximately zero when the neurons are asynchronous. To satisfy this condition, the equation has an extra term 〈H(x)〉, which is the average of H(x) over the oscillators when they are spread over the asynchronous state. Since 〈H(x)〉 is just a constant vector, it is less problematic than the H(x i ) of (8). Restrepo et al. (2006a,b) analyze the coupling strength where synchrony just begins to emerge in the network (i.e., where the asynchronous state becomes unstable). Their key result is that this critical coupling strength depends on the network structure only through the largest eigenvalue l max of the connectivity matrix W. Clearly, this coupling strength also depends on the oscillator models and the manner in which they are coupled, but these influences are a separate factor multiplied by l max . Given two networks with different values of l max but with neurons whose dynamics are governed by the same oscillator model, synchrony will emerge at a lower coupling strength for the network with the larger l max . This separation of the oscillator dynamics from the network structure is in the same spirit as the master stability function analysis and allows us to incorporate l max into our notion of the synchronizability of the network. Restrepo et al. (2007) also derive approximate expressions for the largest eigenvalue l max . If we ignore any correlation between the degrees of neighboring neurons, l max is a simple function of the degree distribution l max i n out

Linking emergence of synchrony to network topology
Given the relationship (6c) between a chain and the covariance of the degree distribution, the maximum eigenvalue can be written in terms of a chain : Hence, l max will be nearly zero when a chain = −1 (its minimum value) and will increase linearly with a chain with slope equal to the mean degree. fact would follow from the Wielandt-Hoffman theorem (Hoffman and Wielandt, 1953), as this theorem shows that (with a suitable reordering of the eigenvalues m i ) Therefore, reproducing an argument from (Zhan et al., 2010), the relative average squared deviation is which will tend to be small. The worst case would be when all indegrees d i in are the same (and equal to d), in which case the bound would be 1/d. Variation in the degrees makes the bound smaller and the eigenvalues m i even closer to the in-degrees d i in . Although our matrices W are not necessarily normal, we have observed that the eigenvalues m i are indeed close to the in-degrees d i in . An example is shown in Figure 2. Thus, we infer that our two measures based on (8) and (11)  ≈ d , and will simply use s m 2 . Given (6a), we can go one step further, and connect the synchronizability measure s m 2 directly with the network statistic a conv . Comparing (6a) to (13), we see that where we have used the mean degree d of the network as a proxy for its expected value E(d). Hence, we expect a higher frequency of the convergent connection motif to decrease synchrony.
Our model diverges from the original Morris-Lecar model, which used Ca ++ rather than Na + , but essentially retains the same framework. We have also adjusted the model parameters as follows to make the resulting phase response curve look more realistic: C = 0.1 nF, g L = 800 pS, E L = −53.24 mV, g Na = 1.822 nS, E Na = 60 mV, g K = 400 pS, E K = −95.52 mV, V 1/2,m = −7.37 mV, k m = 11.97 mV, V 1/2,n = −16.35 mV, k n = 4.21 mV, E syn = 0 mV, t s = 0.5 ms, and s = 20 pA. The synaptic conductance g syn was varied in the individual simulations.

spike rate norMalization using an integral controller
When neurons receive different numbers of synaptic inputs, the disparity in level of synaptic drive can cause neurons to fire at different rates. In some simulations, we desired to compensate for this effect to make the neurons have similar firing rates. Rather than attempt to numerically solve for the input currents to achieve a desired firing rate, we developed an auto-tuning network by implementing an integral controller that adjusts the input current for each neuron until it fires at the desired firing rate. The controller calculates the input current for each neuron at each firing from the error between the interspike interval and that of the desired rate: For all times t during the neuron's j + 1st interspike interval, the current to neuron k is held at the constant (Ogata, 2010). When neuron k spikes again, a new e k (j + 1) is calculated and added to the sum determining I k (t). For our purposes, we set K i = 1, which produced long transients, but results in stable network behaviors. Simulation were run to steady state with the integral controller on and then the current values were held constant (virtually identical results were obtained when controllers were left on through the entire simulation).

Measuring synchrony with the order paraMeter
We use a simple quantification of network synchrony, the Kuramoto order parameter, to measure synchrony in the network simulations. The measure is based on defining a phase u j (t) of each neuron, which ranges from 0 to 2p to indicate the relative state of the neuron through its spiking cycle. For models that are not directly based on the phase, we define the phase as the time since the previous spike, normalized by the average period over the previous five spikes.
To calculate the order parameter, we represent the phase of each neuron by a vector on the unit circle at angle u j (t). The population vector is the average of these unit vectors. The order parameter is the length of the population vector (Strogatz, 2000) r where i = −1 and |·| represents the absolute value of the complex vector. The order parameter ranges between 0 and 1, with r = 1 representing perfect synchrony. The minimum value r = 0 occurs when the network is asynchronous and the phases are evenly distributed around the unit circle. However, r = 0 whenever neuron Models We simulated many neuron models to explore the extent to which our results were model-independent.

Kuramoto model
The Kuramoto model (Kuramoto, 1984), originally proposed by Yoshiki Kuramoto, is a network of continuously coupled oscillators with phase u i (modulo 2π), where ω is the natural frequency and J i (t) is white noise, , which we scale by s. We specify the connectivity strength S relative to the number of oscillators N and the connection probability p of the connectivity matrix W.

PRC model
The phase response curve networks are pulse-coupled networks, meaning that the effect of the pre-synaptic neuron on the postsynaptic neuron's phase only occurs when the pre-synaptic neuron crosses the zero phase (modulo 2π), to simulate the effect of synaptic release at the time of an action potential. At the time of the action potential, the phases u i of the post-synaptic neurons are advanced according to their phase response curves f i (u i ). If a neuron fires, it is unresponsive to any other inputs it receives at that same instant, but its phase remains at u = 0. The PRC model uses the same parameters as the Kuramoto model, except in some simulations we allow the natural frequency v i to depend on neuron i. T j k is the time of kth spike of neuron j. The shape of the phase response curve f i is determined by choice of a i . c(a i ) is chosen to make the maximum value of phase response curve f i be one. We make the average of a i be two to make PRCs that would promote synchrony in a network coupled through excitatory synapses.

Morris-Lecar model
The Morris-Lecar model is a reduced Hodgkin-Huxley like conductance based oscillator (Morris and Lecar, 1981) with two dynamic variables voltage V and channel inactivation n determined by the equations where . The neurons are coupled by synapses that include a synaptic conductance wave-the relative frequency of the chain connection motif (increasing a chain ) should increase synchrony, as it decreases the stability of the asynchronous state. We demonstrate in simulations that the frequency of chains and convergent motifs can influence network synchrony even when a network is not at or near the extremes of perfect synchrony and asynchrony. We present the results in four parts. First, we generate SONETs with a range of second order connectivity statistics and test the predictions of the analysis through simulations of a simple neuron model. Second, to demonstrate that these results can be explained by the analysis, we confirm that critical eigenvalue quantities from the analysis do have the predicted relationships with the network statistics. Third, to assure that our findings are not model dependent, we test networks with several different single neuron models. Fourth, we test the results with a larger range of networks to evaluate the influence of higher order network statistics on the predictions based on second order statistics.

synchrony deterMined by second order connectivity statistics
To test the effects of the second order connectivity statistics on synchrony, we generate 186 SONETs with N = 3000 neurons and 10% connectivity (p = 0.1). In addition to using the standard Erdős-Rényi network (all α's zero), we randomly sample the second order connectivity statistics in the ranges a recip ∈[−1,4], a conv ∈[0,1], a div ∈[0,1], and a chain ∈[−1,1]. These ranges are chosen to include the a's used to reproduce the motif spectrum obtained from neuronal networks published by Song et al. (2005): (a recip , a conv , a div , a chain )≈(3, 0.4, 0.3, 0.2). As hinted by (6c), the valid range of a chain depends on the a conv and a div (Zhao et al., in preparation). We maximize its range by sampling some networks where a chain is set to its maximum and minimum value conditioned on a conv and a div .
We measure synchrony in network simulations of 3000 pulsecoupled PRC neuron models (17) in these 186 SONETs. Three example networks are shown in Figure 4. For each network, we calculate synchrony measured by the Kuramoto order parameter (19) as a function of time, and calculate its average once the network reaches steady state.
In Figure 5 the steady state synchrony is plotted against the measured second order statistics a calculated by (3). No simple relationship between any single connectivity statistic and synchrony is obvious, though it appears the Erdős-Rényi network (black) is among the most synchronous. There appears to be a bimodal distribution of synchrony, suggesting that synchrony rapidly jumps up when some threshold is crossed.
When synchrony is plotted against pairs of connectivity statistics, as shown in Figure 6A, the pattern of dependence on connectivity emerges. The synchrony seems to be a function of (a conv , a chain ) alone, as the level of synchrony varies more or less smoothly in the (a conv , a chain ) projection (lower left of Figure 6) despite the fact that a recip and a div are varying widely across the whole plot. For a given value of a conv , there appears to be a threshold of a chain above which synchrony jumps up, supporting the trend of higher synchrony with a chain predicted by our analysis. Similarly, synchrony jumps down at a threshold value of a conv for a given value of a chain . the population is balanced around the circle, which could also occur when two or more equal clusters are evenly spread around the circle.
Since we are looking at oscillating neurons, we obtain equivalent results if we measure spike count correlation rather than synchrony. As long as one counts spikes in bins sufficiently smaller than the oscillation period, the synchrony induces spike correlations.

results overview
We examine the influence of network structure on the synchrony of neuronal networks. We generate networks using the framework of SONETs which allows us to systematically vary second order connectivity statistics, which are the frequency of reciprocal, convergent, divergent, and chain connections shown in Figure 1 (see Materials and Methods). SONETs extend the commonly used Erdős-Rényi random network model (Erdős and Rényi, 1959), retaining the feature of adding minimal structure beyond what is required to match the connectivity statistics.
The effect of the second order connectivity statistics underlying SONETs are illustrated by the small networks of Figure 3. The networks illustrate how, if one keeps the connection probability fixed, increased convergence or divergence arises from having some neurons with large in-or out-degree, respectively (Figures 3A,B). The networks have many chains when the large in-degree neurons also have large out-degree ( Figure 3D), but they have few chains when neurons with large in-degree have small out-degree and neurons with large out-degree have small in-degree ( Figure 3C).
The analyses of perfect synchrony and asynchrony (see Materials and Methods) predict that two second order statistics of connectivity should play key roles in determining synchrony. First, increasing the relative frequency of the convergent connection motif (increasing a conv ) should decrease synchrony, as it pushes the network state further from complete synchrony. Second, increasing A B C D

Figure 3 | illustrations of networks with different connectivity statistics.
All networks contain N = 9 neurons with N conn = 22 connections, so that the connection probability is ˆ. p ≈ 0 3. For each neuron, darker inner color indicates larger in-degree and darker outer color indicates larger out-degree. (A) Network with large convergence. There is large variance in in-degree. N recip = 4, N conv = 42, N div = 23, and N chain = 45 so that ˆ.
α chain ≈ 0 0 . (B) Network with large divergence (the same network as in panel A, but with connections reversed). There is large variance in out-degree. N recip = 4, N conv = 23, N div = 42, and N chain = 45 so that ˆ.
α div ≈ 0 8, and ˆ. α chain ≈ 0 0. (C) Network with few chains. Neurons with large in-degree tend to have small out-degree, and neurons with large out-degree tend to have small in-degree. N recip = 2, N conv = 33, N div = 35, and N chain = 29 so that ˆ.
α div ≈ 0 4, and ˆ. α chain ≈ 0 4. Note the dependence of the valid range of a chain on the other statistics (bottom row of Figure 6). The relationship between a conv /a div / a chain and the variances/covariance of the in-and out-degree distribution shown in (6) implies that the valid range of a chain increases with a conv and a div (left two panels of last row of Figure 6). The range of a chain is also affected by a recip (lower-right panel of Figure 6). The dependence of the range of a chain and the other statistics creates the appearance of spurious correlations between synchrony and a recip / a div in the top two rows of Figure 6.

eigenvalue analysis
The influence of the connectivity statistics on synchrony was predicted based on analysis of eigenvalues of the network connectivity matrix W and Laplacian matrix L (see Materials and Methods). An analysis of the asynchronous state predicts synchrony will increase with the largest eigenvalue l max of the connectivity matrix, which should increase linearly with a chain [see (15)]. An analysis of the synchronous state predicts synchrony will decrease as the normalized variance σ µ 2 (10) of the Laplacian eigenvalues increases, which should be nearly equal to a conv [see (14)].
In Figure 7, sample spectrum from a few connectivity matrices demonstrate these trends: l max increases with a chain , and only a conv modulates the variance s m 2 of the eigenvalues of the Laplacian L. These observations are confirmed by the plots of l max and s m 2 versus each connectivity statistic for the 186 generated networks (Figure 8). The right top panel shows how l max is nearly completely determined by a chain according to the relationship (15) plotted in red. The second panel in the bottom row of Figure 8 shows that the normalized variance σ µ 2 of the eigenvalues is nearly identical to the a conv and follows the relationship (14) plotted in red. The other relationships between the eigenvalues and connectivity statistics simply reflect how the range of sampled values of a chain depend on the other statistics, as seen in Figure 6.
The tight relationships between a chain and l max and between a conv and σ µ 2 confirm that the influence of these two connectivity statistics on synchrony is indeed through their influence on the key spectral properties of the connectivity matrix W and Laplacian matrix L.  (19) calculated as a function of time, where colors correspond to the associated rastergrams. The networks were chosen to illustrate the range of observed synchrony. The top network (in black) is the Erdős-Rényi random network (a recip = a conv = a div = a chain = 0) which reached high synchrony (steady state r = 0.8). The middle network (in green) that reached moderate synchrony (r = 0.5) was generated from the second order connectivity statistics a recip = −0.2, a conv = 0.7, a div = 0.6, a chain = 0.6. The bottom network (in red) that stayed fairly asynchronous (r = 0.1) was generated from a recip = 0.1, a conv = 0.9, a div = 0.9, a chain = −0.6. For all networks, average connectivity was p = 0.1. PRC model parameters were S = 6, s = 3, a i = 2, and v i = 60 for all neurons i. The connectivity statistics change the average firing rate of the network, as shown in Figure 6B. Little difference in synchrony is observed if we adjust the intrinsic frequency v for each network to keep the average firing rate fixed (not shown). We explore below the effect of fixing not only the firing rate averaged across each network, but also making the firing rate of each neuron be the same.
on linearizations around those states, there is no guarantee those predictions will apply to other networks states, such as intermediate levels of synchrony. Non-linearities in neuron models could lead to model-dependent deviations from those predictions.

dependence of synchrony on single neuron Model
The synchrony analysis, which was based around the completely synchronous and asynchronous states, predicts that the effect of topology on synchrony should be independent of the neuron model. Since these analyses were based All three networks had substantial convergence and divergence (a conv = a div = 0.5) to allow large a chain . Parameters: a chain = 0.5 (increased chains),

synchrony with other networks
All the above simulations were based on SONETs of N = 3000 neurons where we fixed the first order connectivity statistics (p = 0.1) and varied the second order connectivity statistics a. We did not separately manipulate higher order connectivity statistics, as the definition of SONETs does not add additional structure into the connectivity beyond that determined by the second order statistics.
To explore how well a conv and a chain determine synchrony under a wider range of networks, we simulate networks with different first order statistics and networks with different higher order statistics.
To test the effect of changing first order statistics, we decrease the connection probability by a factor of 10 to p = 0.01. We generate 130 SONETs over the same ranges of a as the original p = 0.1 networks and simulate the PRC model of Figure 9B. As shown in the left panel of Figure 10A, the results are similar to those with p = 0.1. We also increase the network size to N = 10000 and set p = 0.03 so that the average degree p(N − 1) is approximately the same as in the original networks. Simulations from 130 SONETs demonstrate that a conv and a chain determine synchrony in the same manner as before (Figure 10A, right).
To manipulate higher order connectivity statistics of the network, we generate networks with power law degree distributions, i.e., with scale-free connectivity. The algorithm for generating the power law networks is described in Materials and Methods. The power law degree distribution has different higher order statistics from the degree distribution of SONETs, for fixed first (p) and second (a) order statistics. By comparing power law networks with SONETs that have the same first and second order statistics, we can evaluate the effect of the higher order statistics.
We exploit relationships (6) between the a's and the degree distribution to generate networks across a range of second order connectivity statistics. We vary the power law exponents of the in-degree and out-degree distributions (7) independently to allow a conv and a div to range in (0,1). To maximize the range of a chain , we let the correlation between in-degree and out-degree range in (−1,1). We test for model-dependence of the results using the three neuron models described in Methods. The simplest model is the Kuramoto oscillator (16). Second, we run additional simulations with pulse-coupled PRC oscillator models (17) under two conditions: one with noise as in Figure 6 and the other with heterogeneity in the model parameters. The third model is a simple conductance based neuronal model based on the . We run all models on the above 186 SONETs and observe how synchrony is affected by changing the coupling strength.
The results shown in Figure 9 demonstrate that for each model, the influence of the network structure on synchrony is indeed determined only by a conv and a chain , as in each case, synchrony appears to be a function of those parameters. The relative influence of a conv and a chain does vary across models. For the Kuramoto model, synchrony is more strongly influenced by a chain than a conv (Figure 9A). For the noisy homogeneous PRC model, the result from Figure 6 is maintained at different connectivity strengths, as both convergence and chains strongly modulated the synchrony (Figure 9B). The noise-free heterogeneous PRC model has a similar dependence on the connectivity parameters ( Figure 9C); replacing noise with variability across the population in intrinsic frequency and PRC shape only alters the connectivity strength at which synchrony emerges. The Morris-Lecar model ( Figure 9D) has a similar dependence on a chain and a conv as the PRC models.
In all cases, increasing a chain or decreasing a conv tends to increase synchrony. The increase in synchrony appears to be concentrated at an effective threshold line through the (a conv, a chain ) plane. Whether a particular network will synchronize or not depends on many factors, but relative changes in a conv or a chain affect synchrony in the same way. Since the effects of changing these two connectivity statistics combine in a model-dependent manner, we obtain a two-dimensional index (a conv, a chain ) for the synchronizability of a network. chrony (Chow, 1998;Denker et al., 2004). We investigate whether or not this heterogeneity is the only mechanism by which increasing a conv decreases synchrony. We repeat the simulation of the noisy PRC model on the SONETs (Figure 9B), where we virtually eliminate the variance in firing rates across neurons by adjusting the intrinsic frequency v i of each neuron to make its average firing rate be close to 10 Hz. Rather than attempting to solve for the 3000 v i , we use an integral controller to creating an auto-tuning network that adjusts the v i (see Materials and Methods, where we refer to v more generically as an input current I). Using this approach, we can reduce the SD in the firing rate across neurons in each network to below 0.1 (compared to firing rate SD than ranged from 0.1 to upward of 10 for the high a conv networks of Figure 9B).
As shown in Figure 11, the effect of a conv is greatly reduced. Synchrony decreases only slightly with increased a conv . The threshold value of a chain where synchrony increases rapidly with a conv depends only weakly on a conv . The maximum value of synchrony for any a chain does still decrease with a conv . This residual effect of a conv cannot be explained by its effect on the heterogeneity of firing rates, as we have eliminated that heterogeneity.

discussion linking synchrony to network Motifs
In this paper we present a framework to relate network connectivity to synchrony in homogeneous neuronal networks. We define network structure through the frequency of second order motifs, which form the basis for four second order statistics of connectivity. These second order connectivity statistics define conditional probabilities of connections between neurons, that is, how a connection from one neuron to a second neuron affects the the probability of (1) a synapse back onto the first neuron (reciprocal), (2) the second neuron receiving inputs from a third neuron (convergent), (3) the first neuron sending out connections to a third neuron (divergent), or (4) either the first neuron receiving inputs or the second neuron sending out connections (chain), as illustrated in Figure 1.
The main result of this paper is that the primary influence of second order network connectivity on synchrony is through the chain and convergent connection motifs. We obtain a two-dimensional index (a conv, a chain ) of synchronizability that captures the effect of the network structure on the synchrony of the neuronal network. Synchrony tends to increase with the relative frequency of chains and decrease with the relative frequency of convergence, where the majority of this change in synchrony occurs abruptly We generate power law networks with N = 3000 neurons and average connectivity p = 0.01. We first generate 160 such power law networks where we set the rising exponents of the degree distributions (7) to g in = g out = 1 so that the degree distributions initially increase linearly as shown by the green trace of Figure 10C. Synchrony for the PRC model and these power law networks is shown in Figure 10B, left. Even for these scale-free networks, synchrony does appear to be a function of just a conv and a chain , as the synchrony varies smoothly in the (a conv, a chain ) projection as before. For the most part, the same qualitative trends persist that were predicted by the synchrony analysis and were observed in the SONETs; synchrony tends to increase with a chain and decrease with a conv . It is clear, though, that the higher order connectivity statistics are modifying the synchrony. The synchrony pattern as a function of second order connectivity statistics in the left panels of Figures 10A,B are not identical despite the model and first order connectivity statistics being the same.
If we increase the steepness of the rise of the degree distribution by setting g in = g out = 10 (red trace of Figure 10C), the dependence of synchrony on the a changes dramatically (Figure 10B, right). Synchrony is still a function of just (a conv , a chain ), but for larger a conv , synchrony no longer increases monotonically with a chain . Although the eigenvalue l max still has the same linear relationship (15) with a chain , this eigenvalue no longer accurately predicts the emergence of synchrony according to the analysis of the stability of the asynchronous state.
We found large values of g in also lead to a non-monotonic dependence of synchrony on a chain for other degree distribution parameters. We hypothesize that this breakdown of the asynchronous state analysis is caused by the presence of many neurons with small incoming degree. One assumption of the stability analysis of the asynchronous state is that all neurons receive many inputs. When g in is large, the number of neurons with small incoming degree (red trace of Figure 10C) is large enough to invalidate the analysis. The results indicate that the moderate number of neurons with even smaller incoming degree in the SONETs (blue trace of Figure 10C) does not invalidate the conclusions of the analysis.

heterogeneity due to convergence
Since increasing a conv increases the variance of the in-degree distribution, neurons in networks with high a conv will have large variability in the number of synaptic connections that they receive. This heterogeneity in the input strength will lead to heterogeneity in the firing rates across neurons in the network, thus decreasing syn-  Figure 9B were repeated where the v i were adjusted by an integral controller to fix the firing rate of each neuron at 10 Hz. The panels correspond to coupling strengths S=3, 6, and 9 of the PRC model (17). modifying the first order statistics of connectivity rather than the second order statistic a conv . They do not introduce the heterogeneity of increased in-degree variance.

independence of synchrony froM coMMon input
One surprising result is that the synchronizability of a network appears to be unaffected by the relative frequency of divergent connection motifs as captured by a div . Intuitively, one would imagine that the common input represented by the divergent connection motif would tend to synchronize neurons, as the common input leads to correlations in their inputs (Shadlen and Newsome, 1998).
In feedforward networks, such correlations are amplified and can lead to synchrony in deeper layers of the network (Diesmann et al., 1999;van Rossum et al., 2002;Reyes, 2003;Tetzlaff et al., 2003;Liu and Nykamp, 2009). To see the relationship between a div and shared input, define b to be the fraction of shared input, i.e., the ratio between the expected number of shared connections between a pair of neurons and the expected number of connections to a single neuron. For SONETs, we can calculate from (2) how b depends on divergence where the last approximation is in the limit of large N.
As an extreme test to see if the common input determined by a div would lead to increased synchrony, we generated SONETs with connection probability p = 0.1 where a div varied over its full range a div ∈[0, 1/p − 1] = [0,9], leading to shared input in the range b∈[0.1, 1]. Such large values of a div are not realistic, as shown by the outdegree distributions of Figure 12A. When a div is near its maximum, the out-degrees are nearly all zero or N − 1: most neurons do not have any projections to others neurons while a small fraction (around pN) project to nearly every other neuron. Nonetheless, we wanted to test to see if such extreme cases, where neurons share almost all of their pre-synaptic connections, would lead to higher synchrony. However, as shown in Figure 12B, even large values of a div had little influence on synchrony. An ensemble average over many networks does unmask a small increase of synchrony with a div (not shown), but this effect is dwarfed by increases in synchrony due to a chain .
For our networks, we obtain similar results if we measure spike correlations rather than synchrony (see Materials and Methods). Other researchers have also found that high common input may not necessarily lead to large spike correlations. Roxin (2011) found that changing common input by manipulating the outgoing degree distribution typically had little impact on correlations in networks of excitatory and inhibitory neurons. Renart et al. (2010) similarly showed that common input had little influence on spiking correlation, but since they used Erdős-Rényi networks where b≈p, they could not isolate changes in common input b from changes in the first order connection probability p. Since our networks do not have inhibition, our results cannot be explained by the opposing effect of inhibition canceling out correlations in input currents (Renart et al., 2010). Instead, we are left to conclude that in these recurrent networks, as opposed to feedforward networks, common input due to divergence only weakly influences correlation, when it is not combined with convergence and chains. along a line in (a conv, a chain ) parameter space that depends on connectivity strength and the neuron model. There is little dependence of synchrony on divergence or reciprocal connections. The result was robust to changes in first order connectivity statistics or network size.
The value of the synchronizability index (a conv, a chain ) stems from the fact that the influence of these network structures on synchrony was similar for a wide variety of neuron models. While there were model specific differences in the simulations, the qualitative dependency of synchrony on the proportion of chains and convergent connections in the networks was the same for most single neuron models and network models. Since the influence of chains and convergence on synchrony was predicted from linearizations around synchrony and asynchrony, it is surprising how well those statistics determined the level of synchrony over the whole range of synchrony. We expect that these linear predictions will not hold for very complicated neuron models, especially with strong coupling. The predictions did break down for networks where many neurons received few inputs, as we discovered parameter regimes where synchrony decreased with a chain .
One result of the two-dimensional synchronizability index is that one cannot arrange networks in order of increasing synchronizability, and hence it may be impossible to predict if some network changes will increase or decrease synchrony without knowing the neuron model. If one network has both a higher a chain and a conv than a second network, it may be that either network may elicit higher synchrony, depending on the neuron model. For example, the Erdős-Rényi network with all a's zero usually synchronized better than other networks, even those with higher a conv and a chain . However, this trend was not model-independent, as we observed some exceptions for low connectivity strength with the Kuramoto model ( Figure 9A) or when we eliminated the firing rate heterogeneity caused by positive a conv (Figure 11).

the decrease of synchrony with convergence
The decrease of synchrony with convergence is not surprising, given that the convergence is correlated with the variance of the in-degree distribution (6a) and therefore results in increased heterogeneity in firing rates. The effect of this in-degree heterogeneity can be reduced by making the synaptic strengths inversely proportional to the in-degree as done by Motter et al. (2005) or by adjusting the current to each cell so that the firing rates are the same, as we have done in our auto-tuned networks. Reducing the heterogeneity greatly reduced the effect of a conv on synchrony but some residual effect of a conv on synchrony remained indicating higher order effects of a conv on synchrony that cannot be accounted for by firing rate heterogeneity alone.
At first glance, the decrease of synchrony with a conv seems to stand in contradiction of recent results by Rosenbaum et al. (2010), where they find that the build up of synchrony in feedforward networks is primarily due to the level of convergence in the network. They show that when neurons receive many inputs (have high d in ), any small correlation among their inputs is amplified. This pooling of many inputs averages out independent fluctuations but not the correlations. Thus, increasing convergence increases correlations and synchrony. The apparent contradiction is resolved by observing that Rosenbaum et al. (2010) change the level of convergence by  (17) was simulated using 7 different coupling strengths S. Since the same networks are used, fluctuations due to the particular network structure are similar for each connectivity strength. Values of α div can exceed the theoretical maximum of 1/p − 1 = 9 as they are estimates from the generated connectivity matrix using (3). PRC model parameters are the same as used in Figure 9B.

the increase of synchrony with chains
Instead of increasing with common input determined by a div , synchrony increases dramatically with chains in the network. This result is consistent with the recent work by LaMar and Smith (2010), who investigated how correlation between in-degree and out-degree decreases the time until complete synchrony in noiseless pulse-coupled PRC models. We hypothesize that the primary influence of chains is that they increase the effective coupling strength of the network. This effect can be seen by the linear dependence of the largest eigenvalue in (15). Intuitively, one can think that signals are propagated along chains in the network, so that increasing the number of chains increases the effect of the connectivity. When a chain is large, neurons with large in-degree tend to have a large out-degree; since in this case, neurons that listen to more neurons tend to talk to more neurons, communication within the network is enhanced. Thus, synchrony emerges at a lower connectivity strength when a chain is increased.
The analysis of the influence of chains can be extended to networks that include inhibition in order to investigate interactions within and between inhibitory and excitatory populations. We plan to study such heterogeneous SONETs to examine the effect of chains involving different patterns of excitatory and inhibitory neurons. We hypothesize that these chains will play a critical role in determining how synchrony is amplified or suppressed.
The frequency of chains in the network may have additional effects beyond the effective coupling strength that contribute to synchrony. For example the influence of chains on loops might play an additional role, as increasing a chain increases the frequency of short loops in the network (Bianconi et al., 2008).

higher order network structure
In the SONETs, the synchronizability of the networks was determined by the second order connectivity statistics of convergence and chains. However, higher order network structure can also influence synchronizability. We introduced higher order connectivity statistics of one particular type: the additional structure contained in the power law in-degree and out-degree distributions of scalefree networks. Adding this higher order structure did alter the synchrony, even if we kept the second order connectivity statistics fixed. This simple manipulation demonstrated that knowing the second order connectivity statistics may not be sufficient to fully know the effect of the network on synchrony. Nonetheless, within a particular class of networks, we did see the same qualitative dependence of synchrony on the second order statistics, where only convergence and chains had an influence. As long as there weren't too many neurons that received few inputs, synchrony tended to increase with chains and decrease with convergence as it did with the SONETs.
The differences observed from the in-degree distribution of the power law networks could likely be captured by adding one-third-order connectivity statistic: the probability of three connections onto a single neuron, which would presumably be linked to the skew of the in-degree distribution in analogy to (6). However, there are many more thirdorder connection motifs (12 in all) so that dipping into third-order statistics would greatly increase the complexity. We expect, though, that only a subset of those third-order statistics would influence synchrony.
Other higher order network statistics that have captured a lot of attention are the clustering coefficient and mean path length often examined in the context of small world networks Watts and Strogatz (1998). A number of studies have looked at the influence of network the frequency of convergent and chain connections. Given this analysis, one can predict how modification of specific connections between neurons in the network affects synchrony. By identifying the critical network structures of convergent and chain connections, one can focus efforts to identify structural changes that may play a role in enhancing or reducing susceptibility to different network states.
One application of this theory may be to understand how changes in neuronal wiring leads to pathological population activity such as abnormal synchrony. Second order statistics may provide a formalism by which to interpret changes in the network structure that may occur in a diseased state and how these changes affect network synchrony.

acknowledgMents
We thank Michael Buice, Chin-Yueh Liu, and Ofer Zeitoni for helpful discussions on the developments of the theory. We thank Alex Roxin for a careful reading of the manuscript. This research was supported by the University of Minnesota's Grant-in-aid (Theoden Netoff and Bryce Beverlin II) and National Science Foundation grants DMS-0719724, DMS-0847749 (Duane Q. Nykamp), and CBET-0954797 (Theoden Netoff). connectivity on neuronal activity and synchrony in the context of small world models (Netoff et al., 2004;Kitano and Fukai, 2007;Kriener et al., 2009). In the classical small world models, all the second order network statistics a are zero except a recip , so one cannot observe the effect of a conv or a chain in these models. For the networks parameters we used, the influence of these second order statistics on synchrony was substantially larger than the influence of the small world rewiring parameter.
While the brain may have higher order connectivity statistics, there are several advantages to using second order statistics to describe a network. In general, lower order statistics require less data to estimate than do higher order statistics; it may be possible to estimate the five SONET parameters from experimental data but harder to estimate the higher order statistics. Furthermore, higher order statistics should be compared to the null hypothesis that they can be explained by the second order statistics observed. Finally, the low-dimensional framework of SONETs provides a mechanism for systematically exploring a wide range of network structures and quantifying their differences.

conclusion
In the current study, SONETs facilitated our discovery that influence of network structure on the level of synchrony in a network of homogeneous excitatory neurons is primarily mediated through