Dynamics of Structured Networks of Winfree Oscillators

1. 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 out-degrees, 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 (Ott and 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 (Bläsche et al., 2020; 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.

2. Model

We consider the network version of the model as presented in Pazó and Montbrió (2014)

$\begin{array}{ll}\frac{d{\theta }_j}{dt}={\omega }_j+U({\theta }_j)\frac{ϵ}{\langle k\rangle }{\displaystyle \sum _{n=1}^N}{A}_{jn}T({\theta }_n)& (1)\end{array}$

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

$\begin{array}{ll}U(\theta )=sin\beta -sin(\theta +\beta )& (2)\end{array}$

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

$\begin{array}{ll}T(\theta )=a_q{(1+cos\theta )}^q& (3)\end{array}$

where q is a positive integer and $a_q=2^q{(q!)}^2/(2q)!$ so that ${\int }_0^{2\pi }T(\theta )d\theta =2\pi$. The in-degree of oscillator j is

$\begin{array}{ll}{k}_{in,j}={\displaystyle \sum _{n=1}^N}{A}_{jn}& (4)\end{array}$

and the out-degree of oscillator n is

$\begin{array}{ll}{k}_{out,n}={\displaystyle \sum _{j=1}^N}{A}_{jn}& (5)\end{array}$

We consider large networks with all oscillators having large in- and 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 $\sum _{\text{k}}P(\text{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

$\begin{array}{ll}\frac{\partial f}{\partial t}+\frac{\partial }{\partial \theta }(vf)=0& (6)\end{array}$

where

$\begin{array}{ll}v(\theta ,\omega ,\text{k},t)=\omega +ϵU(\theta )R(\text{k},t)& (7)\end{array}$

where

$\begin{array}{ll}R(\text{k},t)=\frac{1}{\langle k\rangle }{\displaystyle \sum _{{\text{k}}^{\prime }}}P({\text{k}}^{\prime })a({\text{k}}^{\prime }\to \text{k})G({\text{k}}^{\prime },t)& (8)\end{array}$

and

$\begin{array}{ll}G({\text{k}}^{\prime },t)={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}f({\theta }^{\prime },{\omega }^{\prime }|{\text{k}}^{\prime },t)T({\theta }^{\prime })d{\theta }^{\prime }d{\omega }^{\prime }& (9)\end{array}$

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 (Ott and Antonsen, 2008, 2009). We assume that

$\begin{array}{ll}g(\omega )=\frac{\Delta /\pi }{{(\omega -{\omega }_0)}^2+{\Delta }^2}& (10)\end{array}$

where Δ is the half-width-at-half-maximum and ω₀ 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

$\begin{array}{l}\frac{\partial b(\text{k},t)}{\partial t}=\frac{ϵe^{-i\beta }R(\text{k},t)}{2}+\left[i{\omega }_0-\Delta +iϵsin\beta R(\text{k},t)\right]b(\text{k},t)\\ -\frac{ϵe^{i\beta }R(\text{k},t)}{2}{\left[b(\text{k},t)\right]}^2& (11)\end{array}$

where

$\begin{array}{ll}G(\text{k},t)=a_q\left[C_0+{\displaystyle \sum _{j=1}^q}{C}_j\left\{{\left[b(\text{k},t)\right]}^j+{\left[\bar{b}(\text{k},t)\right]}^j\right\}\right]& (12)\end{array}$

where overline indicates complex conjugate and

$\begin{array}{ll}{C}_j={\displaystyle \sum _{k=0}^q}{\displaystyle \sum _{m=0}^k}\frac{q!{\delta }_{k-2m,j}}{2^k(q-k)!m!(k-m)!}& (13)\end{array}$

The quantity

$\begin{array}{ll}b(\text{k},t)={\displaystyle {\int }_{-\infty }^{\infty }}{\displaystyle {\int }_0^{2\pi }}f(\theta ,\omega |\text{k},t)e^{i\theta }d\theta d\omega & (14)\end{array}$

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 (ω₀, Δ, 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, ω₀. We conclude in section 7.

3. Within Oscillator Correlations

We first consider the effects of correlating an oscillator's in- and 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)

$\begin{array}{ll}a({\text{k}}^{\prime }\to \text{k})=\frac{k_{out}^{\prime }{k}_{in}}{N\langle k\rangle }& (15)\end{array}$

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},\widehat{\rho })$ instead of P(k)/N, where $\widehat{\rho }$ is a parameter controlling the correlation between k_in and k_out, explained in detail below. Substituting (15) into (8) we have

$\begin{array}{l}R(k_{in},k_{out},t)=\frac{N}{\langle k\rangle }{\displaystyle \sum _{k_{in}^{\prime }}}{\displaystyle \sum _{k_{out}^{\prime }}}P(k_{in}^{\prime },k_{out}^{\prime },\widehat{\rho })a({\text{k}}^{\prime }\to \text{k})G(k_{in}^{\prime },k_{out}^{\prime },t)\\ =\frac{k_{in}}{{\langle k\rangle }^2}{\displaystyle \sum _{k_{in}^{\prime }}}{\displaystyle \sum _{k_{out}^{\prime }}}P(k_{in}^{\prime },k_{out}^{\prime },\widehat{\rho })k_{out}^{\prime }G(k_{in}^{\prime },k_{out}^{\prime },t)& (16)\end{array}$

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}^{\prime }$. So we can write

$\begin{array}{ll}R(k_{in},t)=\frac{k_{in}}{{\langle k\rangle }^2}{\displaystyle \sum _{k_{in}^{\prime }}}Q(k_{in}^{\prime },\widehat{\rho })G(k_{in}^{\prime },t)& (17)\end{array}$

where

$\begin{array}{ll}Q(k_{in}^{\prime },\widehat{\rho })\equiv {\displaystyle \sum _{k_{out}^{\prime }}}P(k_{in}^{\prime },k_{out}^{\prime },\widehat{\rho })k_{out}^{\prime }& (18)\end{array}$

Thus, the model equations of interest are

$\begin{array}{l}\frac{\partial b(k_{in})}{\partial t}=\frac{ϵe^{-i\beta }R(k_{in})}{2}+\left[i{\omega }_0-\Delta +iϵsin\beta R(k_{in})\right]b(k_{in})\\ -\frac{ϵe^{i\beta }R(k_{in})}{2}{\left[b(k_{in})\right]}^2& (19)\end{array}$

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.

3.1. Generating Correlated Degrees

The correlations between an oscillator's in- and out-degree are controlled by the function $P(k_{in},k_{out},\widehat{\rho })$ 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.,

$\begin{array}{ll}p(k)=\{\begin{array}{ll}\frac{1}{M-m}\hfill & m\le k\le M\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}& (20)\end{array}$

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

$\begin{array}{ll}f(x,y,\widehat{\rho })=\frac{1}{2\pi \sqrt{1-{\widehat{\rho }}^2}}{e}^{-(x^2-2\widehat{\rho }xy+y^2)/[2(1-{\widehat{\rho }}^2)]}& (21)\end{array}$

where $\widehat{\rho }\in (-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},\widehat{\rho })$ for which the correlations between k_in and k_out out can be varied systematically. The cumulative distribution function for x is

$\begin{array}{ll}C(x)=\left[1+\text{erf}(x/\sqrt{2})\right]/2& (22)\end{array}$

and the cumulative distribution function for degree k is

$\begin{array}{ll}{C}_k(k)={\displaystyle {\int }_m^k}\frac{1}{M-m}ds=\frac{k-m}{M-m}& (23)\end{array}$

The joint degree distribution for k_in and k_out is

$\begin{array}{l}P(k_{in},k_{out},\widehat{\rho })={\{C^{-1}[C_k(k_{in})]\}}^{\prime }{\{C^{-1}[C_k(k_{out})]\}}^{\prime }f\{C^{-1}[C_k(k_{in})],\\ C^{-1}[C_k(k_{out})],\widehat{\rho }\}& (24)\end{array}$

where the superscript “−1” indicates the inverse of the corresponding function. Now

$\begin{array}{ll}{C}^{-1}[C_k(k_{in})]=\sqrt{2}{\text{erf}}^{-1}\left(\frac{2(k-m)}{M-m}-1\right)& (25)\end{array}$

and

$\begin{array}{l}{\{C^{-1}[C_k(k_{in})]\}}^{\prime }=\sqrt{\frac{\pi }{2}}exp\left[{\left\{{\text{erf}}^{-1}\left(\frac{2(k-m)}{M-m}-1\right)\right\}}^2\right]\frac{2}{M-m}\\ =\sqrt{\frac{\pi }{2}}\frac{2}{M-m}exp\left[\frac{{(C^{-1}\left[C_k(k_{in})\right])}^2}{2}\right]& (26)\end{array}$

So

$\begin{array}{l}P(k_{in},k_{out},\widehat{\rho })=\frac{1}{{(M-m)}^2\sqrt{1-{\widehat{\rho }}^2}}\exp \left\{\frac{\widehat{\rho }{C}^{-1}[C_k(k_{in})]C^{-1}[C_k(k_{out})}{1-{\widehat{\rho }}^2}\right\}\\ \times exp\left[\frac{-{\widehat{\rho }}^2\left({\left\{C^{-1}\left[C_k(k_{in})\right]\right\}}^2+{\left\{C^{-1}\left[C_k(k_{out})\right]\right\}}^2\right)}{2(1-{\widehat{\rho }}^2)}\right]& (27)\end{array}$

$\begin{array}{l}=\frac{p(k_{in})p(k_{out})}{\sqrt{1-{\widehat{\rho }}^2}}exp\left\{\frac{\widehat{\rho }{C}^{-1}\left[C_k(k_{in})\right]C^{-1}\left[C_k(k_{out})\right]}{1-{\widehat{\rho }}^2}\right\}\\ \times exp\left[\frac{-{\widehat{\rho }}^2\left({\left\{C^{-1}\left[C_k(k_{in})\right]\right\}}^2+{\left\{C^{-1}\left[C_k(k_{out})\right]\right\}}^2\right)}{2(1-{\widehat{\rho }}^2)}\right]& (28)\end{array}$

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.

FIGURE 1

Figure 1. Joint degree distribution $P(k_{in},k_{out},\widehat{\rho })$ for (A) $\widehat{\rho }=0.5$ and (B) $\widehat{\rho }=-0.5$. The log of P is shown, with red corresponding to higher values and blue to lower. Parameters: m = 100, M = 400.

We also need to relate the parameter $\widehat{\rho }$ to ρ, the Pearson correlation coefficient between the in- and out-degrees of a neuron. We have

$\begin{array}{ll}\rho =\frac{\tilde{\Sigma }P(k_{in},k_{out},\widehat{\rho })(k_{in}-\langle k\rangle )(k_{out}-\langle k\rangle )}{\sqrt{\tilde{\Sigma }P(k_{in},k_{out},\widehat{\rho }){(k_{in}-\langle k\rangle )}^2}\sqrt{\tilde{\Sigma }P(k_{in},k_{out},\widehat{\rho }){(k_{out}-\langle k\rangle )}^2}}& (29)\end{array}$

where $\tilde{\Sigma }$ 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},\widehat{\rho })$ is undefined for k = m, M.

FIGURE 2

Figure 2. (A) Correlation coefficient between the in- and out-degrees of an oscillator, ρ, as a function of the parameter $\widehat{\rho }$ used in the Gaussian copula. (B) The function $Q(k_{in},\widehat{\rho })$ (Equation 18) for different values of $\widehat{\rho }$. Parameters: m = 100, M = 400.

We can also calculate the function $Q(k_{in},\widehat{\rho })$ (Equation 18) where $P(k_{in},k_{out},\widehat{\rho })$ is given in (28). This function is shown in Figure 2B, where we see that increasing $\widehat{\rho }$ 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},\widehat{\rho })$ is the “weight” given to outputs from oscillators with in-degree k_in. If, for example, $\widehat{\rho }>0$, then oscillators with high in-degree will be likely to have high out–degree, and thus their output should be weighted more.

3.2. Results

We set q = 4 (so a_q = 8/35 and C₀ = 35/8, C₁ = 7/2, C₂ = 7/4, C₃ = 1/2, C₄ = 1/16), ω₀ = 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)

$\begin{array}{ll}Z=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{e}^{i{\theta }_j}.& (30)\end{array}$

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 3A,C,E. Note as well the high |Z| value reflects a large fraction of quiescent oscillators in Figures 3C,D—a “trivial synchrony.”

FIGURE 3

Figure 3. Dynamics of the system (1) with uncorrelated degrees. (A,B) Correspond to (ϵ, Δ) = (0.2, 0.05) (synchronous state), (C,D) to (ϵ, Δ) = (0.8, 0.05), and (E,F) to (ϵ, Δ) = (0.2, 0.5) (asynchronous states). The left panels show the magnitude of the global order parameter, and the right show sin θ_j. Other parameters: ω₀ = 1, β = 0, m = 100, M = 400, N = 2000.

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.

FIGURE 4

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.

We now investigate the effects of varying $\widehat{\rho }$ and thus ρ on the dynamics of Equation (19). As mentioned, it is 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 $\widehat{\rho }$ 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).

FIGURE 5

Figure 5. (A) Hopf bifurcation curves of the fixed point of (19). A stable periodic orbit exists below the curve. Other parameters: ϵ = 0.2, ω₀ = 1, m = 100, M = 400. (B) SNIC bifurcation curve (β = 0, blue circles joined by line) and Hopf bifurcation curves (β = 0.5, 0.7 and 1) for (19). For each value of β a stable periodic orbit exists below the curve. Other parameters: Δ = 0.05, ω₀ = 1, m = 100, M = 400.

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 ω₀, β 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).

4. 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 out-degree of the “downstream” (receiving) oscillator (see Figure 6).

FIGURE 6

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ⁱⁿ, or out-degree, k^out, of a node is the relevant quantity).

Roughly speaking, degree assortativity can be thought of in this way: given an upstream oscillator with specific in- and 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. 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

$\begin{array}{ll}r(\alpha ,\gamma )=\frac{{\displaystyle {\sum }_{e=1}^{N_e}{(}^u{k}_e^a-\langle {}^u{{\displaystyle k}}^{\alpha }\rangle ){(}^d{k}_e^{\gamma }-\langle {}^d{{\displaystyle k}}^{\gamma }\rangle )}}{\sqrt{{\displaystyle {\sum }_{e=1}^{N_e}{{(}^u{k}_e^a-\langle {}^u{{\displaystyle k}}^{\alpha }\rangle )}^2}}\sqrt{{\displaystyle {\sum }_{e=1}^{N_e}{{(}^d{k}_e^{\gamma }-\langle {}^d{{\displaystyle k}}^{\gamma }\rangle )}^2}}}& (31)\end{array}$

where

$\begin{array}{ll}\langle {}^{u}k{}^{\alpha }\rangle =\frac{1}{N_e}{\displaystyle \sum _{e=1}^{N_e}{}^u}{k}_e^{\alpha }\text{ and }\langle {}^{d}k{}^{\gamma }\rangle =\frac{1}{N_e}{\displaystyle \sum _{e=1}^{N_e}{}^d}{k}_e^{\gamma },& (32)\end{array}$

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 (Bläsche et al., 2020). For example the upstream node's in-degree of the second edge would be ${}^u{{\displaystyle k}}_2^{in}$. 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)

$\begin{array}{ll}{c}_{\parallel }=(k_j^a-\langle {}^u{{\displaystyle k}}^a\rangle )(k_i^{\gamma }-\langle {}^d{{\displaystyle k}}^{\gamma }\rangle )+(k_l^{\alpha }-\langle {}^u{{\displaystyle k}}^a\rangle )(k_h^{\gamma }-\langle {}^d{{\displaystyle k}}^{\gamma }\rangle )& (33)\end{array}$

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 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 centers ${\widehat{k}}_{in,1},{\widehat{k}}_{in,2},\dots ,{\widehat{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({\widehat{k}}_{out,i},{\widehat{k}}_{in,j},t)$. We can easily convert this to a vector-valued ODE by stacking the columns of $b({\widehat{k}}_{out},{\widehat{k}}_{in},t)$, from left to right, into a vector, $\widehat{b}(t)$, where the sth entry is ${\widehat{b}}_s(t)=b({\widehat{k}}_{out,i},{\widehat{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

$\begin{array}{l}\frac{d{b}_s(t)}{dt}=\frac{ϵe^{-i\beta }{R}_s(t)}{2}+\left[i{\omega }_0-\Delta +iϵsin\beta R_s(t)\right]b_s(t)\\ -\frac{ϵe^{i\beta }{R}_s(t)}{2}{\left[b_s(t)\right]}^2& (35)\end{array}$

for s ∈ {1, 2, …225} where we define

$\begin{array}{ll}{G}_s(t)=a_n\left[C_0+{\displaystyle \sum _{j=1}^n}{C}_j\left\{{\left[b_s(t)\right]}^j+{\left[{\bar{b}}_s(t)\right]}^j\right\}\right]& (36)\end{array}$

We need to calculate R_s(t) from G_s(t) using the equivalent of (8). We can write the analog of (8) as

$\begin{array}{ll}{R}_s(t)=\frac{1}{\langle k\rangle }{\displaystyle \sum _{s^{\prime }=1}^{225}}E(s,s^{\prime })G_{s^{\prime }}(t)& (37)\end{array}$

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

$\begin{array}{ll}R(t)=\frac{1}{\langle k\rangle }EG(t)& (38)\end{array}$

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

$\begin{array}{ll}E(r)=B{r}^2+Cr+D& (39)\end{array}$

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.

4.1. 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).

FIGURE 7

Figure 7. (A–D) Hopf bifurcation curves as both Δ and one of the types of assortativity are varied. A stable periodic orbit exists below the curve. Other parameters: ϵ = 0.2, ω₀ = 1, m = 100, M = 400.

FIGURE 8

Figure 8. (A–D) SNIC bifurcation curve (β = 0, blue circles joined by line) and Hopf bifurcation curves (β = 0.5, 0.7 and 1) as both ϵ and one of the types of assortativity are varied. For each value of β a stable periodic orbit exists below the curve. Other parameters: Δ = 0.05, ω₀ = 1, m = 100, M = 400. In the top-left panel, for β = 0.5 the curve terminates as r is decreased.

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.

5. Correlated Heterogeneity

We have so far assumed that the parameters ω₀ 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 in-degree 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), and Boccaletti 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.

5.1. 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

$\begin{array}{l}\frac{\partial b(k_{in})}{\partial t}=\frac{ϵe^{-i\beta }R(k_{in})}{2}+\left[i{\omega }_0-\Delta +iϵsin\beta R(k_{in})\right]b(k_{in})\\ -\frac{ϵe^{i\beta }R(k_{in})}{2}{\left[b(k_{in})\right]}^2\end{array}$

for k_in = m, m + 1, …M where

$\begin{array}{l}R(k_{in},t)=\frac{k_{in}}{{\langle k\rangle }^2}{\displaystyle \sum _{k_{in}^{\prime }}}{\displaystyle \sum _{k_{out}^{\prime }}}p(k_{in}^{\prime })p(k_{out}^{\prime })k_{out}^{\prime }G(k_{in}^{\prime },t)\\ =\frac{k_{in}}{\langle k\rangle }{\displaystyle \sum _{k_{in}^{\prime }}}p(k_{in}^{\prime })G(k_{in}^{\prime },t)& (41)\end{array}$

where G is given by (12) but with the degree dependence being on only k_in.

We define a scaled in-degree

$\begin{array}{ll}{\widehat{k}}_{in}=2\left(\frac{k_{in}-m}{M-m}\right)-1& (42)\end{array}$

which varies linearly from −1 to 1 as k_in goes from m to M, respectively. We first consider the case where ω₀ is a function of k_in. We write ${\omega }_0(k_{in})=1+\sigma {\widehat{k}}_{in}$ where σ controls the strength of dependence between k_in and ω₀. (Recall that we previously set ω₀ = 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

$\begin{array}{ll}Z(t)=\frac{1}{M-m+1}{\displaystyle \sum _{k_{in}=m}^M}b(k_{in},t).& (43)\end{array}$

This is an appropriate definition since the distribution of in-degrees 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).

FIGURE 9

Figure 9. Intrinsic frequency dependence ${\omega }_0(k_{in})=1+\sigma {\widehat{k}}_{in}$. (A) Solid line: stable fixed point of Equation (40); dashed line: unstable fixed point. Circles: maximum and minimum over one period of Im(Z). Im(Z) calculated using (30) for (B) σ = −0.1, (C) σ = 0, and (D) σ = 0.2. Other parameters: Δ = 0.05, β = 0, ϵ = 0.2, m = 100, M = 400, N = 2, 000.

One might think that having ω₀ 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 ω₀ 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 $\beta =\sigma ({\widehat{k}}_{in}+1)$. We choose ω₀ = 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.

FIGURE 10

Figure 10. Phase dependency $\beta (k_{in})=\sigma ({\widehat{k}}_{in}+1)$. (A) Solid line: stable fixed point of Equation (40); dashed line: unstable fixed point. Circles: maximum and minimum over one period of Re(Z). Re(Z) calculated using (30) for (B) σ = 0, (C) σ = 0.4, and (D) σ = 0.8. Other parameters: Δ = 0.05, ω₀ = 1, ϵ = 0.8, m = 100, M = 400, N = 2, 000.

As a third possibility we let Δ depend on k_in. Δ (the width of the distribution of intrinsic frequencies) cannot be negative so we set $\Delta =0.09+\sigma {\widehat{k}}_{in}$ and consider only −0.09 ≤ σ ≤ 0.09. We set other parameters ω₀ = 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.

FIGURE 11

Figure 11. Heterogeneity dependency $\Delta =0.09+\sigma {\widehat{k}}_{in}$. (A) Solid line: stable fixed point of (40); dashed line: unstable fixed point. Circles: maximum and minimum over one period of Im(Z). Im(Z) calculated using (30) for (B) σ = 0 and (C) σ = 0.09. Other parameters: β = 1, ω₀ = 1, ϵ = 0.6, m = 100, M = 400, N = 2, 000.

5.2. Out-Degree

Now consider the possibility that one of ω₀, Δ, 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

$\begin{array}{ll}\frac{\partial b}{\partial t}=\frac{ϵe^{-i\beta }R}{2}+[i{\omega }_0-\Delta +iϵsin\beta R]b-\frac{ϵe^{i\beta }R}{2}{b}^2& (44)\end{array}$

where b is a function of both k_in and k_out, but R is a function of k_in only:

$\begin{array}{ll}R(k_{in},t)=\frac{k_{in}}{{\langle k\rangle }^2}{\displaystyle \sum _{k_{in}^{\prime }}}{\displaystyle \sum _{k_{out}^{\prime }}}p(k_{in}^{\prime })p(k_{out}^{\prime })k_{out}^{\prime }G(k_{in}^{\prime },k_{out}^{\prime },t)& (45)\end{array}$

where G is given by Equation (12).

5.2.1. Computational Approach

If J = M − m + 1 is the number of distinct in-degrees (and out-degrees) then b can be thought of as a J × J matrix with J² 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

$\begin{array}{ll}(f,g)={\displaystyle \sum _{k=m}^M}f(k)g(k)& (46)\end{array}$

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

$\begin{array}{ll}{\displaystyle \sum _{k_{in}^{\prime }}}{\displaystyle \sum _{k_{out}^{\prime }}}p(k_{in}^{\prime })p(k_{out}^{\prime })k_{out}^{\prime }G(k_{in}^{\prime },k_{out}^{\prime },t)\approx {\displaystyle \sum _{i=1}^s}{\displaystyle \sum _{j=1}^s}{w}_i{w}_j{k}_{j}G(k_i,k_j,t)& (47)\end{array}$

Note that the k_i are not integer-valued. We thus solve (44) on the non-uniform 2D grid of s² “virtual” degree. An example for s = 10 is shown in Figure 12.

FIGURE 12

Figure 12. Non-uniform 2D grid of degrees as explained in section 5.2.1 with m = 100, M = 400, s = 10. The color shows the weight, w_iw_j, associated with each point.

The convergence with s is observed to be geometric (not shown) and we use s = 20 to calculate the results below.

5.2.2. Results

We write

$\begin{array}{ll}{\widehat{k}}_{out}=2\left(\frac{k_{out}-m}{M-m}\right)-1.& (48)\end{array}$

Setting ${\omega }_0(k_{out})=1+\sigma {\widehat{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 $\beta (k_{out})=\sigma ({\widehat{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 $\Delta =0.09+\sigma {\widehat{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).

FIGURE 13

Figure 13. Intrinsic frequency dependence ${\omega }_0(k_{out})=1+\sigma {\widehat{k}}_{out}$. Solid line: stable fixed point of (44); dashed line: unstable fixed point. Circles: maximum and minimum over one period of Im(Z). Other parameters: Δ = 0.05, β = 0, ϵ = 0.2, m = 100, M = 400.

FIGURE 14

Figure 14. Phase shift dependence $\beta (k_{out})=\sigma ({\widehat{k}}_{out}+1)$. Solid line: stable fixed point of (44); dashed line: unstable fixed point. Circles: maximum and minimum over one period of Re(Z). Other parameters: Δ = 0.05, ω₀ = 1, ϵ = 0.8, m = 100, M = 400.

6. Parameter Assortativity

We now consider assortativity by a parameter other than degree, in this case ω₀ 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 ω₀ from a distribution p(ω₀), which is non-zero only if ${\omega }_0\in [{\underset{\_}{\omega }}_0,{\overline{\omega }}_0]$, i.e., ω₀ is the minimum value of ω₀ and ${\overline{\omega }}_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

$\begin{array}{ll}r=\frac{{\displaystyle {\sum }_{e=1}^{N_e}({\omega }_{0,e}^{\prime }-\langle {\omega }_0^{\prime }\rangle )({\omega }_{0,e}-\langle {\omega }_0\rangle )}}{\sqrt{{\displaystyle {\sum }_{e=1}^{N_e}{({\omega }_{0,e}^{\prime }-\langle {\omega }_0^{\prime }\rangle )}^2}}\sqrt{{\displaystyle {\sum }_{e=1}^{N_e}{({\omega }_{0,e}-\langle {\omega }_0\rangle )}^2}}}& (49)\end{array}$

where ${\omega }_{0,e}^{\prime }$ is the value of the target ω₀ associated with the oscillator at the start of edge e and ω_{0, e} is the value of the target ω₀ associated with the oscillator at the end of edge e, and N_e is the number of edges. The means are

$\begin{array}{ll}\langle {\omega }_0^{\prime }\rangle =\frac{1}{N_e}{\displaystyle \sum _{e=1}^{N_e}}{\omega }_{0,e}^{\prime }\langle {\omega }_0\rangle =\frac{1}{N_e}{\displaystyle \sum _{e=1}^{N_e}}{\omega }_{0,e}& (50)\end{array}$

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 ω₀ for that oscillator and with half-width-at-half-maximum Δ. This will result in the creation of a network in which all oscillators have the same in- and out-degree, but those with high ω₀ are more likely to connect to those also having high ω₀ and vice versa.

To derive the continuum equations we see that the state of an oscillator can only depend on its ω₀ value. We discretize the range of ω₀ values, $[{\underset{\_}{\omega }}_0,{\overline{\omega }}_0]$, into m equal-sized bins, and thus we have

$\begin{array}{l}\frac{d{b}_s(t)}{dt}=\frac{ϵe^{-i\beta }{R}_s(t)}{2}+\left[i{\omega }_s-\Delta +iϵsin\beta R_s(t)\right]b_s(t)\\ -\frac{ϵe^{i\beta }{R}_s(t)}{2}{[b_s(t)]}^2& (51)\end{array}$

for s = 1, 2, …m, where ω_s is the value of ω₀ in the center of the sth bin. The analog of Equation (8) is

$\begin{array}{ll}{R}_s(t)=\frac{1}{\langle k\rangle }{\displaystyle \sum _{u=1}^m}{E}_{su}{G}_u(t)& (52)\end{array}$

where 〈k〉 is the degree of each oscillator,

$\begin{array}{ll}{G}_s(t)=a_q\left[C_0+{\displaystyle \sum _{j=1}^q}{C}_j\left\{{[b_s(t)]}^j+{[{\bar{b}}_s(t)]}^j\right\}\right]& (53)\end{array}$

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

$\begin{array}{ll}E(r)=B{r}^2+Cr+D& (54)\end{array}$

where B, C, and D are m × m constant matrices.

As an example we choose β = 0, Δ = 0.01, and p(ω₀) to be the uniform distribution on [0, 2]. (p(ω₀) 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.

FIGURE 15

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(ω₀) 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.

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 ω₀ 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).

7. 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 in-degree, 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 culture—e.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.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This work was partially supported by the Marsden Fund Council from Government funding, managed by Royal Society Te Aparangi, grant number 17-MAU-054.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Anstee, R. (1982). Properties of a class of (0, 1)-matrices covering a given matrix. Can. J. Math. 34, 438–453. doi: 10.4153/CJM-1982-029-3

CrossRef Full Text | Google Scholar

Arenas, A., Díaz-Guilera, A., Kurths, J., Moreno, Y., and Zhou, C. (2008). Synchronization in complex networks. Phys. Rep. 469, 93–153. doi: 10.1016/j.physrep.2008.09.002

CrossRef Full Text | Google Scholar

Ariaratnam, J. T., and Strogatz, S. H. (2001). Phase diagram for the winfree model of coupled nonlinear oscillators. Phys. Rev. Lett. 86:4278. doi: 10.1103/PhysRevLett.86.4278

PubMed Abstract | CrossRef Full Text | Google Scholar

Avalos-Gaytan, V., Almendral, J. A., Papo, D., Schaeffer, S. E., and Boccaletti, S. (2012). Assortative and modular networks are shaped by adaptive synchronization processes. Phys. Rev. E 86:015101. doi: 10.1103/PhysRevE.86.015101

PubMed Abstract | CrossRef Full Text | Google Scholar

Bläsche, C., Means, S., and Laing, C. R. (2020). Degree assortativity in networks of spiking neurons. J. Comput. Dyn. 7, 401–423. doi: 10.3934/jcd.2020016

CrossRef Full Text | Google Scholar

Boccaletti, S., Almendral, J., Guan, S., Leyva, I., Liu, Z., Sendi na-Nadal, I., et al. (2016). Explosive transitions in complex networks' structure and dynamics: percolation and synchronization. Phys. Rep. 660, 1–94. doi: 10.1016/j.physrep.2016.10.004

CrossRef Full Text | Google Scholar

Chandra, S., Hathcock, D., Crain, K., Antonsen, T. M., Girvan, M., and Ott, E. (2017). Modeling the network dynamics of pulse-coupled neurons. Chaos 27:033102. doi: 10.1063/1.4977514

PubMed Abstract | CrossRef Full Text | Google Scholar

De Franciscis, S., Johnson, S., and Torres, J. J. (2011). Enhancing neural-network performance via assortativity. Phys. Rev. E 83:036114. doi: 10.1103/PhysRevE.83.036114

PubMed Abstract | CrossRef Full Text | Google Scholar

de Santos-Sierra, D., Sendi na-Nadal, I., Leyva, I., Almendral, J. A., Anava, S., Ayali, A., et al. (2014). Emergence of small-world anatomical networks in self-organizing clustered neuronal cultures. PLoS ONE 9:e85828. doi: 10.1371/journal.pone.0085828

PubMed Abstract | CrossRef Full Text | Google Scholar

Engblom, S. (2006). Gaussian Quadratures With Respect to Discrete Measures. Technical report, Uppsala University.

Google Scholar

Ermentrout, G. B., and Kopell, N. (1986). Parabolic bursting in an excitable system coupled with a slow oscillation. SIAM J. Appl. Math. 46, 233–253. doi: 10.1137/0146017

CrossRef Full Text | Google Scholar

Fenton, F. H., Cherry, E. M., Hastings, H. M., and Evans, S. J. (2002). Multiple mechanisms of spiral wave breakup in a model of cardiac electrical activity. Chaos 12, 852–892. doi: 10.1063/1.1504242

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, N. I., and Lee, A. (1983). A correlation coefficient for circular data. Biometrika 70, 327–332. doi: 10.1093/biomet/70.2.327

PubMed Abstract | CrossRef Full Text | Google Scholar

Foster, J. G., Foster, D. V., Grassberger, P., and Paczuski, M. (2010). Edge direction and the structure of networks. Proc. Natl. Acad. Sci. U.S.A. 107, 10815–10820. doi: 10.1073/pnas.0912671107

CrossRef Full Text | Google Scholar

Gallego, R., Montbrió, E., and Pazó, D. (2017). Synchronization scenarios in the winfree model of coupled oscillators. Phys. Rev. E 96:042208. doi: 10.1103/PhysRevE.96.042208

PubMed Abstract | CrossRef Full Text | Google Scholar

Gómez-Garde nes, J., Gómez, S., Arenas, A., and Moreno, Y. (2011). Explosive synchronization transitions in scale-free networks. Phys. Rev. Lett. 106:128701. doi: 10.1103/PhysRevLett.106.128701

PubMed Abstract | CrossRef Full Text | Google Scholar

Ha, S.-Y., Park, J., and Ryoo, S. W. (2015). Emergence of phase-locked states for the winfree model in a large coupling regime. Discrete Contin. Dyn. Syst. A 35:3417. doi: 10.3934/dcds.2015.35.3417

CrossRef Full Text | Google Scholar

Hagmann, P., Cammoun, L., Gigandet, X., Meuli, R., Honey, C., Van Wedeen, J., et al. (2008). Mapping the structural core of human cerebral cortex. PLoS Biol. 6:e052306. doi: 10.1371/journal.pbio.0060159

PubMed Abstract | CrossRef Full Text | Google Scholar

Kähne, M., Sokolov, I., and Rüdiger, S. (2017). Population equations for degree-heterogenous neural networks. Phys. Rev. E 96:052306. doi: 10.1103/PhysRevE.96.052306

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2014). Numerical bifurcation theory for high-dimensional neural models. J. Math. Neurosci. 4:1. doi: 10.1186/2190-8567-4-13

PubMed Abstract | CrossRef Full Text | Google Scholar

Laing, C. R. (2017). “Phase oscillator network models of brain dynamics,” in Computational Models of Brain and Behavior, ed A. A. Moustafa (Hoboken, NJ: Wiley Online Library), 505–517. doi: 10.1002/9781119159193.ch37

CrossRef Full Text | Google Scholar

Laing, C. R., and Bläsche, C. (2020). The effects of within-neuron degree correlations in networks of spiking neurons. Biol. Cybern. 114, 337–347. doi: 10.1007/s00422-020-00822-0

CrossRef Full Text | Google Scholar

LaMar, M. D., and Smith, G. D. (2010). Effect of node-degree correlation on synchronization of identical pulse-coupled oscillators. Phys. Rev. E 81:046206. doi: 10.1103/PhysRevE.81.046206

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, W., Wu, Y., Xiao, J., and Zhan, M. (2013). Effects of frequency-degree correlation on synchronization transition in scale-free networks. Europhys. Lett. 101:38002. doi: 10.1209/0295-5075/101/38002

CrossRef Full Text | Google Scholar

Martens, M. B., Houweling, A. R., and Tiesinga, P. H. (2017). Anti-correlations in the degree distribution increase stimulus detection performance in noisy spiking neural networks. J. Comput. Neurosci. 42, 87–106. doi: 10.1007/s10827-016-0629-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Milton, J., and Jung, P. (2013). Epilepsy as a Dynamic Disease. Berlin; Heidelberg: Springer Science & Business Media.

Google Scholar

Nelsen, R. B. (2007). An Introduction to Copulas. New York, NY: Springer Science & Business Media.

Google Scholar

Netoff, T. I., Banks, M. I., Dorval, A. D., Acker, C. D., Haas, J. S., Kopell, N., et al. (2005). Synchronization in hybrid neuronal networks of the hippocampal formation. J. Neurophysiol. 93, 1197–1208. doi: 10.1152/jn.00982.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Newman, M. (2003). The structure and function of complex networks. SIAM Rev. 45, 167–256. doi: 10.1137/S003614450342480

CrossRef Full Text | Google Scholar

Newman, M. (2018). Networks. New York, NY: Oxford University Press. doi: 10.1093/oso/9780198805090.001.0001

CrossRef Full Text | Google Scholar

Nykamp, D. Q., Friedman, D., Shaker, S., Shinn, M., Vella, M., Compte, A., et al. (2017). Mean-field equations for neuronal networks with arbitrary degree distributions. Phys. Rev. E 95:042323. doi: 10.1103/PhysRevE.95.042323

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. (2008). Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18:037113. doi: 10.1063/1.2930766

PubMed Abstract | CrossRef Full Text | Google Scholar

Ott, E., and Antonsen, T. (2009). Long time evolution of phase oscillator systems. Chaos 19:023117. doi: 10.1063/1.3136851

PubMed Abstract | CrossRef Full Text | Google Scholar

Pazó, D., and Gallego, R. (2020). The winfree model with non-infinitesimal phase-response curve: Ott-antonsen theory. Chaos 30:073139. doi: 10.1063/5.0015131

PubMed Abstract | CrossRef Full Text | Google Scholar

Pazó, D., and Montbrió, E. (2014). Low-dimensional dynamics of populations of pulse-coupled oscillators. Phys. Rev. X 4:011009. doi: 10.1103/PhysRevX.4.011009

CrossRef Full Text | Google Scholar

Pazó, D., Montbrió, E., and Gallego, R. (2019). The winfree model with heterogeneous phase-response curves: analytical results. J. Phys. A 52:154001. doi: 10.1088/1751-8121/ab0b4c

CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Kurths, J. (2001). Synchronization. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511755743

CrossRef Full Text | Google Scholar

Restrepo, J. G., and Ott, E. (2014). Mean-field theory of assortative networks of phase oscillators. Europhys. Lett. 107:60006. doi: 10.1209/0295-5075/107/60006

CrossRef Full Text | Google Scholar

Schmeltzer, C., Kihara, A. H., Sokolov, I. M., and Rüdiger, S. (2015). Degree correlations optimize neuronal network sensitivity to sub-threshold stimuli. PLoS ONE 10:e0121794. doi: 10.1371/journal.pone.0121794

PubMed Abstract | CrossRef Full Text | Google Scholar

Schultheiss, N. W., Prinz, A. A., and Butera, R. J. (2011). Phase Response Curves in Neuroscience: Theory, Experiment, and Analysis. New York, NY: Springer. doi: 10.1007/978-1-4614-0739-3

CrossRef Full Text | Google Scholar

Skardal, P. S., Restrepo, J. G., and Ott, E. (2015). Frequency assortativity can induce chaos in oscillator networks. Phys. Rev. E 91:060902. doi: 10.1103/PhysRevE.91.060902

PubMed Abstract | CrossRef Full Text | Google Scholar

Skardal, P. S., Sun, J., Taylor, D., and Restrepo, J. G. (2013). Effects of degree-frequency correlations on network synchronization: universality and full phase-locking. Europhys. Lett. 101:20001. doi: 10.1209/0295-5075/101/20001

CrossRef Full Text | Google Scholar

Strogatz, S. (2000). From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Phys. D 143, 1–20. doi: 10.1016/S0167-2789(00)00094-4

CrossRef Full Text | Google Scholar

Strogatz, S. (2003). Sync: The Emerging Science of Spontaneous Order. Hyperion.

Teller, S., Granell, C., De Domenico, M., Soriano, J., Gómez, S., and Arenas, A. (2014). Emergence of assortative mixing between clusters of cultured neurons. PLoS Comput. Biol. 10:e1003796. doi: 10.1371/journal.pcbi.1003796

CrossRef Full Text | Google Scholar

Tsubo, Y., Teramae, J.-N., and Fukai, T. (2007). Synchronization of excitatory neurons with strongly heterogeneous phase responses. Phys. Rev. Lett. 99:228101. doi: 10.1103/PhysRevLett.99.228101

CrossRef Full Text | Google Scholar

Vasquez, J., Houweling, A., and Tiesinga, P. (2013). Simultaneous stability and sensitivity in model cortical networks is achieved through anti-correlations between the in- and out-degree of connectivity. Front. Comput. Neurosci. 7:156. doi: 10.3389/fncom.2013.00156

PubMed Abstract | CrossRef Full Text | Google Scholar

Vegué, M., Perin, R., and Roxin, A. (2017). On the structure of cortical microcircuits inferred from small sample sizes. J. Neurosci. 37, 8498–8510. doi: 10.1523/JNEUROSCI.0984-17.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

Vegué, M., and Roxin, A. (2019). Firing rate distributions in spiking networks with heterogeneous connectivity. Phys. Rev. E 100:022208. doi: 10.1103/PhysRevE.100.022208

PubMed Abstract | CrossRef Full Text | Google Scholar

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of 'small-world' networks. Nature 393, 440–442. doi: 10.1038/30918

CrossRef Full Text | Google Scholar

Winfree, A. (2001). The Geometry of Biological Time. New York, NY: Springer. doi: 10.1007/978-1-4757-3484-3

CrossRef Full Text | Google Scholar

Winfree, A. T. (1967). Biological rhythms and the behavior of populations of coupled oscillators. J. Theor. Biol. 16, 15–42. doi: 10.1016/0022-5193(67)90051-3

PubMed Abstract | CrossRef Full Text | Google Scholar