Dynamics of Structured Networks of Winfree Oscillators

Winfree oscillators are phase oscillator models of neurons, characterized by their phase response curve and pulsatile interaction function. We use the Ott/Antonsen ansatz to study large heterogeneous networks of Winfree oscillators, deriving low-dimensional differential equations which describe the evolution of the expected state of networks of oscillators. We consider the effects of correlations between an oscillator's in-degree and out-degree, and between the in- and out-degrees of an “upstream” and a “downstream” oscillator (degree assortativity). We also consider correlated heterogeneity, where some property of an oscillator is correlated with a structural property such as degree. We finally consider networks with parameter assortativity, coupling oscillators according to their intrinsic frequencies. The results show how different types of network structure influence its overall dynamics.


INTRODUCTION
The behavior of networks of coupled oscillators is a topic of ongoing interest (Strogatz, 2000;Pikovsky et al., 2001;Arenas et al., 2008). While an individual oscillator may have very simple behavior, it is the emergent behavior such as synchronization that has gained much attention (Winfree, 2001;Strogatz, 2003). Networks of coupled oscillators provide insights into physiological systems such as neuronal or cardiac systems, where synchrony or lack thereof can have profound implications (Fenton et al., 2002;Milton and Jung, 2013).
One of the first models for interacting oscillators was the Winfree model (Winfree, 1967;Ariaratnam and Strogatz, 2001;Pazó and Montbrió, 2014;Ha et al., 2015;Gallego et al., 2017;Pazó et al., 2019;Pazó and Gallego, 2020). Each Winfree oscillator is described by a single angular variable and when uncoupled is assumed to undergo periodic oscillations. Each oscillator is assumed to have a phase response curve, a function of its own phase, which can be measured from individual neurons, for example Schultheiss et al. (2011) and Netoff et al. (2005). This describes how an oscillator's phase changes as the result of input from other oscillators. The output from an oscillator is assumed to be in the form of a non-negative pulsatile function of its own phase, and the inputs to an oscillator are assumed to be additive.
A number of authors have studied networks of Winfree oscillators, but as far as we are aware, only in the all-to-all coupled case. Although straightforward to assemble, such networks do not reproduce complex network structures observed in real-world systems such as assortativities between individual neurons (de Santos-Sierra et al., 2014;Teller et al., 2014). We are interested in networks with far more varied structure, not just randomly connected. These networks are directed, i.e., edges connect one oscillator to another, without necessarily having a reciprocal connection, as occurs in networks of neurons. There are many ways to create structured networks and here we consider the following: correlating the in-and out-degrees of an oscillator (i.e., the number of inputs and the number of outputs of an oscillator, section 3), inducing degree assortativity (i.e., connecting two oscillators based on their in-and outdegrees, section 4), correlating some local property of an oscillator with either its in-or out-degree (section 5), and inducing parameter assortativity (i.e., connecting two oscillators based on the similarities of an intrinsic property of the two oscillators, such as their free-running frequency, section 6).
Our main tool is the derivation and then numerical analysis of moderately large sets of coupled ordinary differential equations (ODEs). The derivation utilizes the Ott/Antonsen ansatz Antonsen, 2008, 2009), an exact technique for dimension reduction in large networks of sinusoidally coupled phase oscillators, of which Winfree oscillators are an example. Some of the computational techniques used here have been presented before Laing and Bläsche, 2020) but for networks of theta neurons (Ermentrout and Kopell, 1986) rather than Winfree oscillators. In section 2, we present the general model and its reduction using the Ott/Antonsen ansatz. The results are presented in sections 3-6 and we conclude in section 7.

MODEL
We consider the network version of the model as presented in Pazó and Montbrió (2014) for j = 1, . . . N where the ω j are chosen from a distribution g(ω), ǫ is the strength of coupling, k is the mean degree of the network, and the connectivity of the network is given by the adjacency matrix A, where A jn = 1 if oscillator n connects to oscillator j and zero otherwise. The function U is known as the phase response curve and we choose it to be U(θ ) = sin β − sin (θ + β) so that U(0) = 0. If β < π/2 then this function describes a type-II oscillator whereas β = π/2 describes a type-I oscillator (Tsubo et al., 2007). We consider only type-II oscillators in this work. The pulsatile function T is given by where q is a positive integer and a q = 2 q (q!) 2 /(2q)! so that 2π 0 T(θ )dθ = 2π. The in-degree of oscillator j is and the out-degree of oscillator n is We consider large networks with all oscillators having large inand out-degrees. Following Chandra et al. (2017) and Laing and Bläsche (2020), we assume that the network can be characterized by two functions: the degree distribution P(k), where k = (k in , k out ), and k in and k out are the in-and out-degrees of an oscillator, respectively, and an assortativity function a(k ′ → k) giving the probability that an oscillator with degree k ′ connects to one with degree k, given that such oscillators exist. Note that we follow (Chandra et al., 2017;Laing and Bläsche, 2020) and normalize P(k) such that k P(k) = N.
In the limit N → ∞ the network is described by the probability density function f (θ , ω|k, t) where f (θ , ω|k, t)dθ dω is the probability that an oscillator with degree k has phase in [θ , θ + dθ ] and value of ω in [ω, ω + dω] at time t. This function satisfies the continuity equation where and The nature of this system [specifically, having U(θ ) being a single sinusoidal function of θ ] means that it is amenable to the use of the Ott/Antonsen ansatz Antonsen, 2008, 2009). We assume that where is the half-width-at-half-maximum and ω 0 the median of the distribution of intrinsic frequencies. Using standard techniques (Chandra et al., 2017;Laing, 2017) which rely on the Ott/Antonsen theory, one can show that the long-time dynamics of the network is described by where overline indicates complex conjugate and The quantity is the complex-valued order parameter for oscillators with degree k. Equation (11) is the general equation describing the dynamics of the network and we use it as a base for analysing a number networks with different types of structure. In section 3, we consider correlations between an individual oscillator's in-degree and its out-degree, as described by the degree distribution P(k). In section 4, we consider correlations between the degrees of connected oscillators, effectively modifying the function a(k ′ → k). In section 5, we investigate the results of one of the parameters intrinsic to an oscillator (ω 0 , , or β) being correlated with a network property of that oscillator (its in-or out-degree). Section 6 considers the case when all oscillators have the same in-and out-degrees, and the assortativity function a(k ′ → k) is replaced by a function describing the probability of connecting oscillators based on the values of one of their intrinsic parameters-in this case, ω 0 . We conclude in section 7.

WITHIN OSCILLATOR CORRELATIONS
We first consider the effects of correlating an oscillator's inand out-degree. This general question has been considered by a number of authors studying different types of oscillators (LaMar and Smith, 2010;Vasquez et al., 2013;Martens et al., 2017;Nykamp et al., 2017;Vegué and Roxin, 2019) and experimental evidence for within-neuron degree correlations is given in Vegué et al. (2017). Our derivation follows Laing and Bläsche (2020).
Assuming neutral assortativity we have (Restrepo and Ott, 2014) where we have assumed that the largest in-and out-degrees are significantly smaller than N, so that a(k ′ → k) ≤ 1. We will write P(k in , k out ,ρ) instead of P(k)/N, whereρ is a parameter controlling the correlation between k in and k out , explained in detail below. Substituting (15) into (8) we have This is clearly independent of k out , thus v must also be independent of k out , the state of an oscillator with degree (k in , k out ) must be independent of k out , and thus G must be independent of k ′ out . So we can write where Thus, the model equations of interest are where G is given by (12) but with the degree dependence being on only k in . Note that the model equations are independent of N, the total number of oscillators.

Generating Correlated Degrees
The correlations between an oscillator's in-and out-degree are controlled by the function P(k in , k out ,ρ) and we now describe how to generate these correlations. For simplicity we assume that the distributions of the in-and out-degrees are the same, namely uniform distributions between m and M, i.e., We introduce correlations between the in-and out-degree of an oscillator while retaining these marginal distributions, using a Gaussian copula (Nelsen, 2007). The correlated bivariate normal distribution with zero mean is whereρ ∈ (−1, 1) is the correlation between x and y. The variables x and y have no physical meaning and we use the copula just as a way of deriving an analytic expression for P(k in , k out ,ρ) for which the correlations between k in and k out out can be varied systematically. The cumulative distribution function for x is and the cumulative distribution function for degree k is The joint degree distribution for k in and k out is where the superscript "−1" indicates the inverse of the corresponding function. Now and Note the special case P(k in , k out , 0) = p(k in )p(k out ), as expected. Several plots of this function are shown in Figure 1.
We also need to relate the parameterρ to ρ, the Pearson correlation coefficient between the in-and out-degrees of a neuron. We have (29) where˜ indicates a sum over all k in and k out . This relationship is numerically determined and shown in Figure 2A, and it is nearly the identity. Note that the sums in (29) are over m + 1 ≤ k ≤ M − 1, since P(k in , k out ,ρ) is undefined for k = m, M.
We can also calculate the function Q(k in ,ρ) (Equation 18) where P(k in , k out ,ρ) is given in (28). This function is shown in Figure 2B, where we see that increasingρ gives more weight to high in-degree nodes and less to low in-degree nodes and vice versa. This can be understood by realizing that Q(k in ,ρ) is the "weight" given to outputs from oscillators with in-degree k in . If, for example,ρ > 0, then oscillators with high in-degree will be likely to have high out-degree, and thus their output should be weighted more.

Results
We set q = 4 (so a q = 8/35 and C 0 = 35/8, C 1 = 7/2, C 2 = 7/4, C 3 = 1/2, C 4 = 1/16), ω 0 = 1, and consider four different values of β: 0, 0.5, 0.7, and 1 (all corresponding to type-II oscillators). There are two types of behavior typically seen in such a network: synchronous and asynchronous (Pazó and Montbrió, 2014), although the fraction of oscillators actually oscillating can vary in the asynchronous states. Increasing ǫ (the strength of coupling) tends to destroy synchronous behavior through a saddle-node-on-invariant-circle (SNIC) bifurcation, as many of the oscillators "lock" at an approximate fixed point. Increasing (the spread of intrinsic frequencies) tends to destroy synchronous behavior through a Hopf bifurcation, as the oscillators become too dissimilar to synchronize (Pazó and Montbrió, 2014). Examples of typical behavior in a default network are shown in Figure 3. The global order parameter for a network of N phase oscillators is a measure of their synchrony, and is defined as (Strogatz, 2000) We see that its magnitude has large, nearly periodic oscillations in the synchronous state, but is approximately constant in the asynchronous state-note the different vertical scales in Figures    The network whose behavior is shown in Figure 3 was created using the configuration model (Newman, 2003). Such a network typically has both self-connections (i.e., an oscillator is connected to itself) and multiple connections from one particular oscillator to another. We remove these in a random way as shown in Figure 4. For a self-connection we randomly choose another connection and reconnect as in the top panel. For a double connection we randomly choose another connection and reconnect as in the bottom panel.
We now investigate the effects of varyingρ and thus ρ on the dynamics of Equation (19). As mentioned, it is FIGURE 4 | (Top) We remove a self-connection to oscillator a (left) by rewiring the randomly chosen connection from oscillator b to c, giving the configuration at the right. (Bottom) we remove the double connection from oscillator a to b by rewiring the randomly chosen connection from oscillator c to d, giving the configuration at the right.
known that increasing (making the intrinsic frequencies more diverse) destroys the synchronous state in a supercritical Hopf bifurcation (Pazó and Montbrió, 2014). In Figure 5A, we show how the value of at which this bifurcation occurs varies as a function of ρ, for four different values of β. We varyρ but interpolate the data shown in Figure 2A in order to plot the curves in Figure 5A as functions of ρ. We see that increasing ρ increases the value of at which the bifurcation occurs, at least for small β, and vice versa, but the effect is small compared with that of varying β. Put another way, for a fixed value of , increasing ρ can cause macroscopic oscillations within the network (at least for β close to zero).
We now fix = 0.05 and consider the effects of varying both ρ and ǫ (the strength of coupling between oscillators). It is known that for an all-to-all coupled network increasing ǫ destroys the synchronous state in a SNIC bifurcation (Pazó and Montbrió, 2014). For our network this is also what happens for β = 0, as shown in Figure 5B (blue circles joined by line). However, for β = 0.5, 0.7, and 1, there is instead a supercritical Hopf bifurcation as ǫ increases, in contrast with the situation for all-to-all coupled network (for these values of ω 0 , β and ), illustrating a nontrivial effect of network structure: even the type of bifurcation occurring is changed. These curves of Hopf bifurcations are also shown in Figure 5B and we see that increasing ρ decreases the value of ǫ at which the synchronous solution is destroyed and vice versa. Note that between β = 0 and β = 0.5, guided by the results for the fully-connected network (Pazó and Montbrió, 2014), we expect there to be several curves of Hopf, homoclinic, and saddle-node bifurcations in Figure 5B organized around a Takens-Bogdanov and a saddle-node separatrix-loop point (Gallego et al., 2017), but we have not investigated them here. The results in Figure 5 have been compared with those from simulation of the full network (Equation 1) and found to agree very well (results not shown).

BETWEEN NEURON DEGREE CORRELATIONS
We now turn to the question of correlations between connected oscillators based on their degrees, often referred to as degree assortativity (Foster et al., 2010;Bläsche et al., 2020). Assortativity has often been studied in undirected networks, where a node simply has a degree, rather than in-and out-degrees (Restrepo and Ott, 2014). Here we consider directed networks, which a small number of previous authors have considered (De Franciscis et al., 2011;Avalos-Gaytan et al., 2012;Schmeltzer et al., 2015;Kähne et al., 2017), although they have often imposed other structure on the network such as equal in-and out-degrees for each model neuron.
Because we consider directed networks there are four possible types of degree assortativity, between either the in-or out-degree of the "upstream" (sending) oscillator and either the in-or outdegree of the "downstream" (receiving) oscillator (see Figure 6).
Roughly speaking, degree assortativity can be thought of in this way: given an upstream oscillator with specific inand out-degrees, and a downstream oscillator with specific in-and out-degrees, one can calculate the probability of a connection from the upstream to the downstream oscillator.
FIGURE 6 | Assortativity in undirected and directed networks. An undirected network (left) is assortative if high degree nodes are more likely to be connected to high degree nodes, and low to low, than by chance (top left). Such a network is disassortative if the opposite occurs (bottom left). (Here, the degree of a node is given by k.) In directed networks (right) there are four possible kinds of assortativity. The probability of a connection (red) depends on the number of red shaded links of the sending (left) and receiving (right) node. (Here, either the in-degree, k in , or out-degree, k out , of a node is the relevant quantity).
If this probability-averaged over the network-is other than that expected by chance, and is further dependent on the degrees of the oscillators, the network shows degree assortativity. One can use this idea to create networks with assortativity, by creating connections where they would typically not occur.
A measure of assortativity for a network with a given connectivity matrix A is by way of calculating the four Pearson correlation coefficients r(α, γ ) with α, γ ∈ [in, out] given by where N e being the number of edges and the leading superscript u or d refers to the "upstream" or "downstream" oscillator of the respective edge . For example the upstream node's in-degree of the second edge would be u k in 2 . Note that there are four mean values to compute.
To induce assortativity within a network we start by randomly choosing in-degrees and out-degrees from the distribution given in Equation (20). If the total number of out-degrees does not equal that of the in-degrees (i.e., the network cannot be created; Anstee, 1982) we choose again until it does. We then use the configuration model (Newman, 2003) with these prescribed degrees to create the network, and utilize the same procedure as described earlier for removal of self-and multiple-connections (see Figure 4).
To induce assortativity of the form (α, γ ) we randomly choose two edges, one connecting oscillator j to oscillator i and another connecting oscillator l to oscillator h. We calculate their contribution to the numerator of (31) and the contribution if we replaced these two edges with one connecting oscillator j to oscillator h and another connecting oscillator l to oscillator i: If c \ / > c we make the swap, otherwise we do not. We then repeat this process many times, storing A, and calculating the value of r(α, γ ) at regular intervals.
We now discuss how to implement the system Equation (11). Choosing m = 100, M = 400, k in , and k out take on values in {100, 101, 102, . . . 400} and thus there are 301 × 301 possible values of k. Considering that we use a network of size N = 2, 000 it is clear that there may be many values of k for which there is not even one oscillator in the network. Thus, we coarse-grain by degree: we divide the interval [100, 400] into 15 equal-size bins with centersk in,1 ,k in,2 , . . . ,k in,15 and describe the state of an oscillator by the value of b associated with the 2D bin it is in (there are 15 × 15 of these 2D bins). We can think of Equation (11) as being a matrix-valued ODE, with the (i, j)th element of the matrix being b(k out,i ,k in,j , t). We can easily convert this to a vectorvalued ODE by stacking the columns of b(k out ,k in , t), from left to right, into a vector,b(t), where the sth entry isb s (t) = b(k out,i ,k in,j , t) and s = i + 15(j − 1). Note that i, j ∈ {1, 2, . . . 15} and s ∈ {1, 2, . . . 225}.
Dropping the hat on b we have We need to calculate R s (t) from G s (t) using the equivalent of (8).
We can write the analog of (8) as where E(s, s ′ ) encodes the connectivity from the 2D bin with index s ′ to that with index s. Given the connectivity matrix A it is straightforward to calculate E(s, s ′ ) as explained in Bläsche et al. (2020). E can be thought of as a 225 × 225 matrix, with (i, j)th entry E(i, j), so we can write Equation (37) as where R and G are vector-valued variables and Equation (35) is just the sth component of a vector-valued ODE.
Since we have recorded A at discrete values of the correlation coefficient r, we can also calculate E at these values. To form a parameterized family, E(r), we fit a quadratic to each entry of E as a function of r, i.e., we write E ij (r) = B ij r 2 + C ij r + D ij for i, j ∈ [1, 225], using linear least-squares. We can then efficiently evaluate an approximation of E(r) as where the (i, j)th entry of B is B ij etc. In summary, we have a parameterized set of ODEs, where r is one of the parameters. Note that we only vary one of the four r(α, γ ) at a time.

Results
The results are shown in Figure 7, where we vary and the four r(α, γ ) for four different values of β, and Figure 8, where we vary ǫ and the four r(α, γ ) for the same four values of  β. As was seen in Bläsche et al. (2020), assortativities of the type r(out,out) and r(out,in) have no discernable effect on the bifurcations, whereas the other two types do. We can understand this by realizing that the dynamics of an oscillator depend only on its inputs. Since an oscillator's dynamics are independent of its downstream oscillators, neither the r(out,out) nor the r(out,in) assortativities influence the overall network dynamics as shown in all the traces of Figures 7C,D. Note, this dynamic interplay is quite different for a network with strong preferential attachment between high in-degree and high out-degree oscillators as when r(in,out) is positive (Figure 7B). The influence of the upstream oscillator (with high in-degree, receiving multiple inputs) is amplified or "passed on" to more oscillators via its downstream companion with high out-degree. This pair with high input and high output is thus far more influential than, say, a pair of oscillators preferentially attached according to the upstream node's out-degree. In that scenario, the upstream node of an attached pair may only integrate a small number of inputs (low in-degree), whose behavior is strikingly distinct from an oscillator with many inputs (high in-degree).
Analogously, a positive r(in,in) assortativity demonstrates preferential attachment between high in-degree upstream and downstream pairs of oscillators. In this case, they are relatively potent integrators and concentrators of upstream impulses. We see in Figure 7A, the influence of high r(in,in) where the parameter space in which stable periodic orbits exist shrinks, increasing sensitivity to the destructive influence of on synchrony.

CORRELATED HETEROGENEITY
We have so far assumed that the parameters ω 0 and (the mean and width, respectively, of the distribution of intrinsic frequencies, see Equation 10) and β (the parameter in the phase response curve, see Equation 2) are the same for each oscillator, but now consider the case of them being correlated with a structural property of an oscillator such as its indegree or out-degree. Correlating an oscillator's frequency with its degree is known to cause "explosive" synchronization in undirected networks of coupled phase oscillators, for example Liu et al. (2013), Gómez-Gardeñes et al. (2011), andBoccaletti et al. (2016), and we are interested in whether similar effects occur in networks of Winfree oscillators. For simplicity we will use linear relationships between a parameter and its relevant degree.

In-Degree
We first consider the case of correlation with in-degree. Assuming neutral assortativity and independence between an oscillator's in-and out-degree, following the reasoning in section 3 the dynamics of b depend only on in-degree and are governed by where G is given by (12) but with the degree dependence being on only k in . We define a scaled in-degreê which varies linearly from −1 to 1 as k in goes from m to M, respectively. We first consider the case where ω 0 is a function of k in . We write ω 0 (k in ) = 1 + σk in where σ controls the strength of dependence between k in and ω 0 . (Recall that we previously set ω 0 = 1.) Setting β = 0, ǫ = 0.2, and = 0.05 we find that when σ = 0 the network is attracted to a stable periodic orbit. However, increasing or decreasing σ causes the oscillations to cease through a Hopf bifurcation as shown in Figure 9A. To visualize the oscillations we define the complex order parameter for (40) as This is an appropriate definition since the distribution of indegrees is uniform; if it were not we would have to weight the contributions from different k in values. Figure 9A shows the maximum and minimum over one period of Im(Z) for oscillatory solutions, and just Im(Z) for fixed points. Simulations of a finite network are shown in the lower panels of Figure 9 (transients not shown) which confirm the results in Figure 9A.
The small amplitude oscillations seen in Figures 9A,C are a result of finite size fluctuations about the fixed point of Equation (40), the linearization about which has complex eigenvalues. The amplitude of these oscillations decreases as the number of oscillators used increases (not shown). One might think that having ω 0 depend on in-degree broadens the distribution of intrinsic frequencies in the network, which is equivalent in some sense to increasing . However, it is not completely equivalent for several reasons. Firstly, the distribution of all intrinsic frequencies is no longer Lorentzian (although for each oscillator we choose the frequency from a Lorentzian), and depends on both the form of dependence of ω 0 on k in (linear in this case) and the distribution of the k in (uniform in this case). Secondly, the intrinsic frequency of each oscillator now depends on a structural property: its in-degree. But for comparison, the oscillations seen in Figure 9 for σ = 0 are destroyed in a Hopf bifurcation as is increased through ∼ 0.085 (not shown).
Next consider β being a function of k in . In order to not have negative β we set β = σ (k in + 1). We choose ω 0 = 1, ǫ = 0.8, = 0.05. For these parameters the network is attracted to a stable fixed point. However, increasing σ first induces oscillations through a SNIC bifurcation and then destroys them through a Hopf bifurcation, as shown in Figure 10A. Simulations of a finite network are shown in the lower panels of Figure 10 and these are consistent with the results in Figure 10A.
As a third possibility we let depend on k in . (the width of the distribution of intrinsic frequencies) cannot be negative so we set = 0.09 + σk in and consider only −0.09 ≤ σ ≤ 0.09. We set other parameters ω 0 = 1, β = 1 and ǫ = 0.6. A Hopf bifurcation occurs as σ is increased as shown in Figure 11A. Simulations of a finite network are shown in the lower panels of Figure 11. Significant oscillations are seen for σ = 0, and the amplitude of oscillations for σ = 0.09 is less than expected. However, we repeated this type of simulation with N = 5, 000 and found that the amplitude of oscillations with σ = 0.09 better matched the results in Figure 11A (i.e., were bigger than seen for N = 2, 000) and that the amplitude of oscillations for σ = 0 were slightly smaller than seen for N = 2, 000 (results not shown), suggesting that this apparent disagreement is a finite-size effect.

Out-Degree
Now consider the possibility that one of ω 0 , , and β are correlated with an oscillator's out-degree. From Equation (11), it is clear that even for neutral assortativity, b will depend on both k in and k out . Thus, the relevant equations are where b is a function of both k in and k out , but R is a function of k in only: where G is given by Equation (12).

Computational Approach
If J = M − m + 1 is the number of distinct in-degrees (and outdegrees) then b can be thought of as a J × J matrix with J 2 entries. This is too many to deal with computationally, so we discretize in degree space. In the same way that one can approximate a definite integral using Gaussian quadrature, it is possible to approximate a double sum like that in (45) using a double sum over far fewer points (Engblom, 2006). The theory is explained in Laing and Bläsche (2020), but put briefly we define an inner product on either degree space and assume that there is a corresponding set of orthogonal polynomials {q n (k)} 0≤n associated with this product. We choose a positive integer s and let {k i } i=1,...s be the roots of q s , found using the Golub-Welsch algorithm, and {w i } be the weights associated with these roots. The approximation of the double sum in (45) is then (47) Note that the k i are not integer-valued. We thus solve (44) on the non-uniform 2D grid of s 2 "virtual" degree. An example for s = 10 is shown in Figure 12.
The convergence with s is observed to be geometric (not shown) and we use s = 20 to calculate the results below.

Results
We writek Setting ω 0 (k out ) = 1 + σk out and varying σ we obtain the results in Figure 13. Two Hopf bifurcations are seen, as in Figure 9A but at different values of σ from in that figure. Writing β(k out ) = σ (k out + 1) and varying σ we obtain the results in Figure 14.
The bifurcations are the same as in Figure 10A, but again, at different values of σ . Using the parameters shown in Figure 11A, setting = 0.09 + σk out and varying σ ∈ [−0.09, 0.09] (as cannot be negative) the fixed point was always stable (not shown). Simulations of a discrete network of N = 2, 000 oscillators confirmed all of the results in this section (not shown).

PARAMETER ASSORTATIVITY
We now consider assortativity by a parameter other than degree, in this case ω 0 value. We first describe how to create a network with such assortativity, then derive the relevant continuum equations. We follow Skardal et al. (2015) in our derivation.
To create a particular network we first create a network where the in-and out-degrees of all oscillators are the same, in order that degree not affect the dynamics. To do this we use the configuration model (Newman, 2003), then remove all self-connections and multi-edges as before. With N oscillators we randomly choose N target values of ω 0 from a distribution p(ω 0 ), which is non-zero only if ω 0 ∈ [ω 0 , ω 0 ], i.e., ω 0 is the minimum value of ω 0 and ω 0 is the maximum, and assign these to oscillators. We can calculate the assortativity of the network using similar ideas as those in section 4. We calculate the Pearson correlation coefficient where ω ′ 0,e is the value of the target ω 0 associated with the oscillator at the start of edge e and ω 0,e is the value of the target ω 0 associated with the oscillator at the end of edge e, and N e is the number of edges. The means are Initially the network will have r ≈ 0. We induce assortativity in a similar way to that described in section 4. We randomly chose two edges, one connecting oscillator j to oscillator i and another connecting oscillator l to oscillator h. We calculate their contribution to the numerator of Equation (49) and the contribution if we replaced these two edges with one connecting oscillator j to oscillator h and another connecting oscillator l to oscillator i. If performing this swap increases r we make the swap, otherwise we do not. We then repeat this process many times, storing A and calculating the value of r at regular intervals. (To decrease r from its initial value of 0 we just consider whether making the swap decreases r). As a last step, in order to use the Ott/Antonsen ansatz, we then randomly assign to oscillator i a value of ω i chosen from a Lorentzian with mean equal to the target ω 0 for that oscillator and with half-width-at-halfmaximum . This will result in the creation of a network in which all oscillators have the same in-and out-degree, but those with high ω 0 are more likely to connect to those also having high ω 0 and vice versa.
To derive the continuum equations we see that the state of an oscillator can only depend on its ω 0 value. We discretize the range of ω 0 values, [ω 0 , ω 0 ], into m equal-sized bins, and thus we have for s = 1, 2, . . . m, where ω s is the value of ω 0 in the center of the sth bin. The analog of Equation (8) is where k is the degree of each oscillator, and the matrix E encodes the connectivity of the network, i.e., E su is proportional to the number of oscillators in the uth bin which connect to oscillators in the sth bin, which can be determined from the connectivity matrix A. As in section 4, we record A at discrete values of the correlation coefficient r, so can construct E(r) at those values. We fit a quadratic through each entry of E as a function of r and thus write E(r) = Br 2 + Cr + D where B, C, and D are m × m constant matrices.
As an example we choose β = 0, = 0.01, and p(ω 0 ) to be the uniform distribution on [0, 2]. (p(ω 0 ) must have bounded support so we can discretize its domain into a finite number of bins.) We compare the results of simulating a full network from Equation (1) with those from the reduced model in Equation (51). We use a network of N = 2, 000 with each oscillator having degree k = 100. We have stored the connectivity matrix A at 101 values of r, and vary both ǫ and r. At each point in this parameter space we solve Equation (1) for 100 time units, discard the first 50 as transients, then calculate the order parameter using Equation (30). The difference between the maximum of |Z| over the final 50 time units and the minimum of |Z| over this time is shown in Figure 15.
When this difference is close to zero, most of the oscillators are "locked" at zero frequency, but for ǫ = 0.8 there is a transition at r ≈ 0 where some the oscillators start unlocking, with those having largest ω 0 unlocking first. Note that this is not a "classical" bifurcation, as the system is not at fixed point before this transition. However, solving the reduced Equations (51) we find that there is a stable fixed point to the left of the red curve in Figure 15 which is destroyed in a Hopf bifurcation, leading to periodic and then quasiperiodic behavior as r is increased. Thus the reduced model provides an explanation for the observed behavior of the full model (1).
The results in Figure 15 are an example of the types of results we can obtain using the framework presented here. We could vary parameters other than ǫ, or introduce assortativity by another intrinsic parameter, β. In this case we would have to use a different measure of correlation between the β values for connected oscillators, as β is an angular variable (Fisher and Lee, 1983).
FIGURE 15 | Difference between the maximum of |Z| over the last 50 time units out of 100 and the minimum, having already discarded the first 50 as transient. p(ω 0 ) is uniform on [0, 2]. The red curve shows the Hopf bifurcation of the steady state of (51) which is stable to the left of this curve. Other parameters: β = 0, = 0.01, k = 100, N = 2, 000. We use m = 20 bins to calculate the blue curve.

CONCLUSION
We studied large directed networks of Winfree oscillators under the assumption that the expected dynamics of an oscillator in such a network is determined by its degree: either its indegree, out-degree, or both (apart from the homogenous degree networks in section 6). Using the Ott/Antonsen ansatz we find that the dynamics are given by Equation (11). Correlations between the in-and out-degree of an oscillator were introduced using a Gaussian copula in section 3, where we investigated the influence of these correlations on the position of bifurcations destroying stable periodic orbits. In section 4, we investigated four types of degree assortativity, as in Bläsche et al. (2020), and found similar results, viz. two types of assortativity have no effect on the network dynamics, while the other two do. Correlations between an oscillator's intrinsic parameter and either its in-or out-degree were examined in section 5. Parameter assortativity was considered in section 6. The framework presented here is quite general, and we believe it to be a powerful method for investigating the general issue of the influence of a network's structure on its dynamics.
The main tool used was numerical continuation and bifurcation analysis of a large number of coupled ODEs (Laing, 2014), which enabled the determination of bifurcation points as parameters were varied-including correlations between various network properties. Following such bifurcations shows the influence of network properties in their dynamics.
The influence of these correlations and assortativities on network synchrony is complex, nuanced, and multi-faceted. Introducing degree correlations within oscillators subtly shapes sensitivity of the network to oscillator parameters such as variability of intrinsic frequencies, , and coupling strength, ǫ. Assortativities between the degrees of connected oscillators can have similar effects for in-degree correlations-r(in,·)-or none at all-r(out,·). More dramatically, for the parameters considered, inducing correlations between oscillator degree (in or out) and intrinsic frequency, ω, destroys oscillatory synchrony. Similarly, if the degree and phase offset, β, are correlated, this may cause or destroy synchronized oscillations. The theme continues with assortativities between intrinsic frequencies, where if they are excessively assortative, oscillators in the network unlock from the population-requiring a higher coupling level to stay locked. Conversely, correlating an oscillator's degree with the width of the distribution from which its intrinsic frequency is chosen, , has little effect.
Network structure such as preferential attachment between similar (or dissimilar) oscillators and the influence we have observed here in idealized systems may reflect structural influences in physiological networks of neurons. Intrinsic connectivity preferences observed of neurons grown in culturee.g., similar numbers of synaptic or dendritic processes connected to each other in groups-results in strong assortativity patterns (de Santos-Sierra et al., 2014;Teller et al., 2014) further inferred in the human cerebral cortex (Hagmann et al., 2008). Our observations of network structure influencing the overall synchrony of a network may be a structural means of calibrating the dynamics of physiological neurons.
Regarding section 5, correlating degree with intrinsic frequency is known to cause explosive synchronization, characterized by bistability between asynchronous and partially synchronized states, in undirected networks of Kuramoto phase oscillators (Gómez-Gardeñes et al., 2011;Liu et al., 2013). We did not observe such behavior but we only considered uniform degree distributions (not power law; Gómez-Gardeñes et al., 2011;Liu et al., 2013) and have directed connections, not undirected. Also, there are many ways to correlate an intrinsic parameter with a degree (Skardal et al., 2013); our form of modification keeps the parameter for nodes with mean degree the same and increases/decreases the parameter for those with degrees above/below mean (or vice versa) in a linear way.
We certainly do not yet have a full understanding of the possible dynamics of the network (1). Possible extensions of the work here include simultaneously having more than one type of structure present in the network (for example, both within-oscillator degree correlations and degree assortativity) or correlating an oscillator's intrinsic parameter with some other network property such as the oscillator's centrality (Newman, 2018) or local clustering coefficient (Watts and Strogatz, 1998). More detailed knowledge about the connectivity in networks of neurons of interest would provide motivation to study these extensions, and help verify some of our results.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.