# Emergence of global synchronization in directed excitatory networks of type I neurons

Abolfazl Ziaeemehr, Mina Zarei, Aida Sheshbolouki

► To cite this version:

Abolfazl Ziaeemehr, Mina Zarei, Aida Sheshbolouki. Emergence of global synchronization in directed excitatory networks of type I neurons. Scientific Reports, 2020, 10 (1), pp.3306. 10.1038/s41598-020-60205-0. hal-03599122

HAL Id: hal-03599122

https://hal.science/hal-03599122v1

Submitted on 5 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN

# Emergence of global synchronization in directed excitatory networks of type I neurons

Abolfazl Ziaeemehr<fnref n="1" />, Mina Zarei<fnref n="1" />,<fnref n="2" />* & Aida Sheshbolouki<fnref n="3" />

The collective behaviour of neural networks depends on the cellular and synaptic properties of the neurons. The phase-response curve (PRC) is an experimentally obtainable measure of cellular properties that quantifies the shift in the next spike time of a neuron as a function of the phase at which stimulus is delivered to that neuron. The neuronal PRCs can be classified as having either purely positive values (type I) or distinct positive and negative regions (type II). Networks of type 1 PRCs tend not to synchronize via mutual excitatory synaptic connections. We study the synchronization properties of identical type I and type II neurons, assuming unidirectional synapses. Performing the linear stability analysis and the numerical simulation of the extended Kuramoto model, we show that feedforward loop motifs favour synchronization of type I excitatory and inhibitory neurons, while feedback loop motifs destroy their synchronization tendency. Moreover, large directed networks, either without feedback motifs or with many of them, have been constructed from the same undirected backbones, and a high synchronization level is observed for directed acyclic graphs with type I neurons. It has been shown that, the synchronizability of type I neurons depends on both the directionality of the network connectivity and the topology of its undirected backbone. The abundance of feedforward motifs enhances the synchronizability of the directed acyclic graphs.

For several decades, there has been a continuing research interest in the synchronization phenomenon due to its widespread application in natural and artificial systems ranging from neural dynamics<fnref n="1" />⁻<fnref n="3" />, cardiac pacemaker cells⁴, and power grid networks⁵ to social networks⁶. Particularly, synchronization plays a key role in the proper functionality of brain neurons. The synchronization of neuronal networks has been associated with many cognitive processes including memory formation⁷, directed attention⁸,⁹, and the processing of sensory stimuli¹⁰.

Synchronization patterns mainly depend on the dynamical properties of individual oscillators as well as the underlying structural connectivity. It has been found that brain neurons have different types of intrinsic dynamics that close to threshold, they may be grouped into two excitability classes: Type I and Type II. These two types of neurons are different in terms of the bifurcation structure observed during their transition to firing. Type I oscillations arise via a saddle-node-onto-limit-cycle bifurcation, whereas type II oscillations are initiated by Andronov-Hopf bifurcation. In addition, type I neurons exhibit a continuous frequency-current curve, whereas type II neurons show a discontinuous frequency-current curve¹¹,¹².

The aforementioned excitability type of a single neuron can be quantified by the phase response curve (PRC) which is an experimentally obtainable measure based on the transient change in the cycle period of the neuron in response to an external stimulus¹³. Different profiles of the neuronal phase response curve arise for different types of neurons in that in type I neurons, any excitatory perturbation causes an acceleration of the next spike, while in type II neurons, perturbations cause acceleration or delay of the next spike depending on the phase at which the perturbation is delivered to that neuron. These qualitatively different responses to stimulation lead to dramatically different synchronization patterns in neural networks. Previous studies have shown that, type I cells exhibit relatively poor propensity for synchronization under excitatory couplings, while type II cells are synchronized

<footnote n="1">¹Institute of Advanced Studies in Basic Sciences (IASBS), Department of Physics, Zanjan, 45137-66731, Iran.</footnote>

<footnote n="2">²Institute for Research in Fundamental Sciences (IPM), School of Computer Science, Tehran, 19395-5531, Iran.</footnote>

<footnote n="3">³University of Waterloo, Cheriton School of Computer Science, Waterloo, N2L3G1, Canada. *email: mina.zarei@iasbs.ac.ir</footnote>
better<sup>11,14-16</sup>. Experimental results verify the neurons' ability to switch between types as well as the coexistence of both cell types in the brain<sup>17-19</sup>. For instance, changes in the neocortical synchronization during sleep-wake cycle are often associated with excitability changes in the cortical neurons<sup>20,21</sup>.

Regarding the underlying network connectivity, Master Stability Function is a well-known formalism providing estimation of the synchronization stability in the networks of coupled identical oscillators<sup>22</sup>. According to this formalism, the spread of the eigenvalues of the Laplacian matrix is a synchronizability index. That is, the more compact the eigenvalues, the more likely the synchronous state will be stable. This index works well when comparing synchronizability of chaotic oscillators and identical Kuramoto oscillators. Using the index, it can be shown that complete graphs and directed acyclic graphs with identical node in-degrees are optimal networks provided that they embed an oriented spanning tree<sup>23,24</sup>.

An entirely different approach to explore the impact of network structure on synchronization is investigating important and dominant topological features of networks known as motifs<sup>25-28</sup>. Motifs are significantly over represented subgraphs that have been recognized as the building blocks of many real-world networks in various domains<sup>29</sup>. Therefore, it's important to understand how the global dynamics are affected by the network motifs. For example, previous studies have shown that feedback loop motifs contribute to two phenomena: responses to noise and dynamical stability<sup>30,31</sup>.

In this paper, we investigate the role of network motifs in synchronization of directed neural networks with type I or type II neurons. Among various network motifs, we focus on two distinct categories: (I) feedback loops (FBL) and (II) feedforward loop motifs (FFL). These motifs are commonly highlighted features observed in real-world networks and have attracted a great deal of attention in the literature<sup>32,33</sup>. In order to study the dynamics of a large collection of such motifs, we generate two different directed graphs, with the same undirected skeleton, constructed only from FFLs or FBLs. These directed graphs are called directed acyclic graph (DAG) and balanced directed graph (BDG), respectively. Both analytical and simulation results show that type I excitatory and inhibitory neurons are synchronized in directed acyclic networks (i.e. networks without any feedback loops). However, these neurons fail to synchronize in undirected networks. In fact, feedback loops destroy synchronization of the networks with identical type I neurons.

# Material and Methods

## Structure of the networks.
An orientation of an undirected graph is an assignment of a single direction to each of its edges, turning the initial graph into a directed graph. In this paper, using two different algorithms, we built different orientations of a graph such that the number of feedback loops (FBLs) is maximized in one orientation and minimized in the other one. These algorithms are described below.

It has been shown that, the number of feedback loops increases by enhancing the correlation between the in- and out-degrees of the nodes<sup>34</sup>. Therefore, we start with an undirected graph, and try to assign each edge a single direction in such a way that the in-degree and the out-degree of any given node are almost the same. The Eulerian cycle in an undirected graph is a cycle that visits every link exactly once. An undirected graph has an Eulerian cycle if and only if every vertex has an even degree. Eulerian orientation of a graph is an orientation that directs the edges along the Eulerian path such that every vertex has the same number of incoming and outgoing edges. The resulting directed graph is called a balanced directed graph (BDG) and has maximum number of feedback loops. In a network in which some nodes are of odd degree, we can still build an orientation in a way that the resulting directed graph is almost balanced. To this end, a new fake node is added to the undirected graph and any existing node with an odd degree is connected to this new node. Since every undirected graph has an even number of nodes with odd degree, all the nodes in the new undirected graph have even degree and we can build the Eulerian orientation of the new graph. After orienting the graph, the fake node is removed with all its adjacent links. In this manner, the difference between in- and out-degree of each node is at most one and the oriented graph is almost balanced<sup>35</sup>.

Directed acyclic graphs (DAGs) are constructed by extending the residual degree gradient method<sup>36,37</sup>. At first, each node is labeled by a residual degree which determines the number of adjacent undirected links. The residual values drop as more adjacent links are assigned directions. Thus, initial and final residual values for each node equal to its degree and zero, respectively. In each iteration of assigning the directions, the node with the smallest residual degree is selected and all of its adjacent undirected links are assigned incoming directions. The residual values are updated afterwards and iterations are repeated until there are no undirected edges left in the network. Since each node will have incoming links only from the previously chosen nodes, the constructed graph will be acyclic. It is probable that this method gives rise to a disconnected DAG which in turn leads to an incomplete synchronization. Therefore, to avoid drawing incorrect conclusions when comparing synchronizability of different graphs, we have considered just connected DAGs in our simulations.

It should be noted that, directed and undirected networks should have the same number of nodes and arcs (i.e. unidirectional links) to be comparable in terms of their synchronizability. To this end, we double the number of arcs while giving directions to the links of an undirected graph. Therefore, the adjacency matrix of an undirected graph is defined as $a_{ij} = 1$, if nodes i and j are connected, and $a_{ij} = 0$, otherwise. Whereas, the adjacency matrix of a directed graph is defined as, $a_{ij} = 2$, if a link is directed from node j to node i, and $a_{ij} = 0$, otherwise.

## Phase response curve (PRC).
The phase response curve (PRC) is an illustrative tool to determine the phase shift of an oscillating neuron in response to a brief current pulse delivered at various phases of the cell cycle<sup>15</sup>. It can be defined as, $PRC = 1 - \frac{T_\theta}{T_0}$. Where $T_0$ is the unperturbed cycle period of the neuron, and $T_\theta$ is the cycle period when the neuron is perturbed at phase $\theta$. Therefore, positive and negative values of PRC indicate phase advances and delays, respectively. As mentioned before, neurons are classified into two excitability classes: type I and type II. Type I neurons always respond to a small excitatory stimulus by advancing the next spike. Therefore, they have positive values of PRC for all phases. In contrast, type II neuron, can advance or delay the next spike depending on the phase at which the perturbation is delivered to them. Therefore, for type II neurons, the PRC
versus $\theta$ diagram, contains both positive and negative regions. Studies verify that the Hodgkin-Huxley neurons have PRC II excitability type, while the Wang-Buszáki and Traub neurons have PRC I excitability type³⁸.

**Dynamics of the neurons.** Neurons can be modeled using both phase oscillators and conductance-based models. Since phase models are simple enough to be mathematically tractable, the networks of coupled phase oscillators have been widely studied to model biological and physical systems including neuronal networks. On the other hand, the detailed investigation of the neuronal dynamics and systematic parameter variation is possible through the use of conductance-based models. Here, we choose the extended Kuramoto model (simulating dynamics of coupled phase oscillators) and Wang-Buszáki and Traub neuron model (simulating the dynamics of coupled spiking neurons) to consider different aspects in the synchronization of neuronal networks. In the following, each model is explained with its corresponding synchronization order parameter.

**Phase model.** We use the extended Kuramoto model³⁹,⁴⁰ to analyse the dynamics of directed networks with type I or type II oscillators. This model assumes a network as a collection of N coupled phase oscillators such that the evolution of the ith oscillator is given by:

$$ \dot{\theta}_i = \omega_i + \frac{K}{N} \sum_{j=1}^{N} a_{ij} G(\theta_i, \theta_j), \quad G(\theta_i, \theta_j) = u_i \sin(\theta_j - \theta_i) + (1 - u_i) \frac{(1 - \cos(\theta_j - \theta_i))}{2} \quad (1) $$

where $\theta_i$ is the phase of the ith oscillator and K is the overall coupling strength. The adjacency matrix of the graph is given by $A = (a_{ij})$. Type I and type II oscillators are distinguished by $u_i = 0$ and $u_i = 1$, respectively. The function $G(\theta)$ corresponds to the phase response curve of the oscillator (see the Fig. S1). In our simulations, the initial values of $\theta_i$ are randomly drawn from a uniform distribution in the interval [0, $2\pi$], and natural frequencies are identical. The degree of synchrony of the phase oscillators is quantified by the Kuramoto order parameters R, which is defined as $r(t) = \frac{1}{N} \langle |\sum_{i=1}^{N} e^{i\theta_i(t)}| \rangle$. Here, $\langle \dots \rangle$ represents averaging over different network realizations and initial conditions. The magnitude is $0 \le r \le 1$. The extreme cases are $r = 1$ (coherent state) and $r = 0$ (incoherent state). The time average of r after achieving a steady state is symbolized by R.

**Neuronal model.** We used Wang-Buszáki model to analyse the dynamics of the networks with type I inhibitory interneurons. The Wang-Buszáki neuron model includes the following differential equation⁴¹:

$$ \begin{aligned} C_m \frac{dv}{dt} &= -I_{Na} - I_K - I_L - I_{syn} + I_{app}, \\ \frac{dh}{dt} &= \phi(\alpha_h(1 - h) - \beta_h h), \\ \frac{dn}{dt} &= \phi(\alpha_n(1 - n) - \beta_n n), \\ \frac{ds}{dt} &= \alpha F(v)(1 - s) - \beta s. \end{aligned} $$

where $C_m$ is the capacitance of the neuron and $v$ is the neuron's membrane potential. $h$ and $n$ are time-varying activation variables that depend on voltage-dependent rate functions ($\alpha_h, \beta_h, \alpha_n, \beta_n$), $s$, $\alpha$ and $\beta$ are activation variable and rate functions for synapses. $\phi$ is the control parameter which directly alters the time constants of sodium inactivation and potassium activation. $F$ is a scaling factor which is given by $F(v) = 1/(1 + \exp(-(v - \theta_{syn})/2))$. Here, $\theta_{syn}$ is a parameter which must be set high enough to have neurotransmitter release only when a spike has been emitted in the presynaptic neuron. $I_{Na}$, $I_K$, $I_L$, $I_{syn}$, and $I_{app}$ represent sodium ($Na^+$), delayed rectifier potassium ($K^+$), leakage, synaptic and external currents, respectively. The currents are modeled by the following equations:

$$ I_L = g_L(v - E_L), \quad I_{Na} = g_{Na} m_{\infty}^3 h(v - E_{Na}), \quad I_K = g_K n^4 (v - E_K), \quad I_{syn,i} = \frac{K}{N} \sum_{j=0}^{N_i} g_{syn} s(v_j - E_{syn}). $$

where $g_L$, $g_{Na}$, $g_K$, and $g_{syn}$ are the maximal values of the conductances for leak, sodium, potassium and synaptic currents, respectively, and $E_L$, $E_{Na}$, $E_K$, and $E_{syn}$ are their respective reversal potentials. K is the coupling strength and the summation in $I_{syn,i}$ occurs over all neighbours of the ith neuron. The steady state activation parameter of the sodium current is $m_{\infty} = \alpha_m/(\alpha_m + \beta_m)$. The parameters of the Wang-Buszáki model⁴¹ are summarized in the Table S1. The rate functions are given by the following functions:

$$ \begin{aligned} \alpha_m &= -0.1(v + 35)/(\exp(-0.1(v + 35)) - 1), \\ \alpha_h &= 0.07\exp(-(v + 58)/20), \\ \alpha_n &= -0.01(v + 34)/(\exp(-0.1(v + 34)) - 1), \\ \beta_m &= 4\exp(-(v + 60)/18), \\ \beta_h &= 1/(\exp(-0.1(v + 28)) + 1), \\ \beta_n &= 0.125\exp(-(v + 44)/80). \end{aligned} $$
We have also used the the Traub pyramidal excitatory type I neuron model$^{42}$. The equations describing the model is $C_m \frac{dv}{dt} = -I_{Na} - I_K - I_L - I_{syn} + I_{app}$, where gating variables m, n and h obey equations of the type $\dot{x} = \alpha_x(v)(1-x) - \beta_x(v)x$. The rate functions are:

$$\alpha_m = 0.32(v + 54)/(1 - \exp(-(v + 25)/4)),$$

$$\alpha_h = 0.128\exp(-(v + 50)/18),$$

$$\alpha_n = 0.032(v + 52)/(1 - \exp(-(v + 52)/5))$$

$$\beta_m = 0.28(v + 27)/(1 - \exp((v + 27)/5) - 1),$$

$$\beta_h = 4/(1 + \exp(-(v + 27)/5)),$$

$$\beta_n = 0.5\exp(-(v + 57)/40).$$

Synaptic gates satisfy $\dot{s} = \alpha(1-s)(1 + \exp(-(v+10)/10)) - \beta s$, and synaptic current to the ith neuron is evaluated from $I_{syn,i} = K \sum_{j \in N_i} g_{syn} s(V_{res})$, where $N_i$ are neighbors of the ith neuron and $V_{res}$ is resetting membrane potential. The parameters of the Traub neuron model, which corresponds to AMPA-receptor-mediated synapse, are summarized in the Table S2.

There are various measures to quantify the level of synchrony in a large population of neurons within a network. Among the various metrics, we choose the interspike distance synchrony measure, denoted as B. As one would expect, the minimum interval between spikes of different neurons in synchronized state is less than that of asynchronous state. Based on this idea, the spike synchrony measure is defined as follows$^{43}$:

$$B = \left( \frac{\sqrt{\langle \tau_x^2 \rangle_{\tau} - \langle \tau_x \rangle_{\tau}^2}}{\langle \tau_x \rangle} - 1 \right) \frac{1}{\sqrt{N}} \quad (2)$$

where, N is the number of neurons in the network and $\tau_x = t_{x+1} - t_x$ denotes the interspike interval of the merged set of network spikes. Here, the inter-spike intervals are calculated between different neurons. $\langle \dots \rangle_{\tau}$ denotes the averaging over all intervals. It can be shown that B is bounded between 0 and 1, where $B = 1$ and $B = 0$ represent fully synchronized and asynchronous states, respectively. The voltage synchrony$^{44}$ is another commonly used measure that has been investigated in the supplementary material.

## Results

First, we studied the synchronization stability of the FBLs and FFLs connecting identical type I or type II phase oscillators analytically. For different motifs, we considered the reduced two-dimensional systems of the phase differences (the equations are provided in the supplementary material). As expected, we found stable (unstable) synchronous states for type II excitatory (inhibitory) FFLs and FBLs. On the other hand, we observed stable asynchronous state for type I excitatory and inhibitory FBLs, and stable synchronous state for type I excitatory and inhibitory FFLs. In fact, both eigenvalues of the Jacobian matrix, for the synchronous state of the type I FFL, are zero. Consequently, this fixed point is globally stable, but not asymptotically. It means that, if one perturbs the system and then leaves the system alone, the time dependence of the perturbation is milder than exponential. The phase portraits of the type I and type II excitatory and inhibitory oscillators situated in different motifs are plotted in the Fig. 1.

We also numerically investigated the synchronization stability of the type I and II identical excitatory phase oscillators including different motifs (detailed explanation of the method is found in the supplementary material and the article written by Wolf et al.$^{45}$). The evolutions of the three Lyapunov exponents for the type I and the type II FBLs and FFLs are demonstrated in the first and the third rows of the Fig. 2. The steady Lyapunov exponents versus coupling strengths are plotted at the second and the fourth rows of the Fig. 2. We can see that, as we expected from the analytical results, all the three Lyapunov exponents for the type I FFLs are zero. The bottom panels of the Fig. 2 illustrate the evolution of the stationary order parameter for the type I and II excitatory phase oscillators connected by FBLs and FFLs with various coupling strengths. As anticipated, the results display the stable synchronous state for the type II FBLs and FFLs, and the type I FFLs. Therefore, the feedforward loop motif is an optimal structure for synchronization of the identical excitatory type I phase oscillators. This is also correct for the identical inhibitory type I phase oscillators (see Fig. 1).

Figure 3 shows the structural properties of the BDGs and DAGs constructed from the same undirected Watts-Strogatz small-world and Barabási-Albert scale-free networks$^{46,47}$. The adjacency matrices have been exhibited in the first row. The adjacency matrix of a directed acyclic graph with a single source node (i.e. a node with zero in-degree) can be transformed into strictly lower triangular form by a permutation matrix. Therefore, as we can see, our orientation algorithm constructs DAGs correctly from different undirected graphs. The Laplacian matrix of directed network is defined as $L = D_{in} - A$, where $D_{in}$ is a diagonal matrix of in-degrees and A is the adjacency matrix. Real and imaginary parts of the eigenvalues of the Laplacian matrices are shown in the bottom panels. The real parts and the imaginary parts of the complex eigenvalues have their own information as the real parts reflect the undirected topology of the network and the imaginary parts reflect the effects of directed links and loops$^{35,48}$. For example, since the eigenvalues of a lower triangular matrix are the same as its diagonal entries, the Laplacian eigenvalues of a DAG are real and equal to the in-degrees of its nodes. Therefore, as we expected DAGs do not have any FBL and their eigenvalues are real. On the other hand, we can see that the variance of the imaginary part of the small-world BDGs is higher than that of the scale-free BDGs. It is because of the fact that clustering coefficient of a small-world network is higher than clustering coefficient of scale-free networks
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Phase plane portraits for reduced two-dimensional systems of type I and type II motifs. The phase plane portraits for reduced two-dimensional systems of (top row) type I and (bottom row) type II phase oscillators situated in a FBL, and FFL motifs. The plots on the left-hand side (right-hand side) of the figure correspond to the excitatory (inhibitory) oscillators. The stable and unstable (or saddle) fixed points have been distinguished by red and green colours, respectively. The Jacobian matrices of the fixed points coloured in blue have two zero eigenvalues.</figcaption></figure>

with the same number of nodes and arcs. Therefore, small-world BDGs have higher number of FBLs than scale-free BDGs. In addition, because of the hubs (i.e. large-degree nodes in the undirected graphs), the real parts of the eigenvalues of the scale-free BDGs are more spread out than small-world BDGs. On the other hand, the real parts of the eigenvalues of the scale-free DAGs are compacter than small-world DAGs, because they have narrower in-degree distributions. Note that based on our orientation algorithm the source node would be a hub in a scale-free DAG. Now that we understood the structural differences between our directed networks, in the following we will compare the synchronizability of type I and type II oscillators on these structures.

As mentioned, unlike FBLs, FFLs favour synchronization in networks of type I oscillators. However, in real world networks motifs are merged together to form more complex structures. In order to see the effects of FFLs and FBLs in the synchronization of the large networks, we constructed two different directed networks from the same undirected backbone. Balanced Directed Graphs (BDGs) with many FBLs, and Directed Acyclic Graphs (DAGs) without any FBL. The flow of information in the BDGs is similar to their undirected backbone, but the information flow in the DAGs is directed from the source node to the rest of the graph. Therefore, BDGs and DAGs are very different, even if they have the same undirected skeleton.

The synchronizability of the scale-free networks with different link-directionalities has been investigated in the Fig. 4. We can see that type II oscillators are synchronized on all the three network structures. However, type I excitatory oscillators can be synchronized only on DAGs. This proves that FBLs are harmful for synchronization of type I oscillators, while FFLs enhance their synchronization. In addition, we found similar synchronization level for BDGs and undirected graphs. This result is in agreement with previous finding that the synchronization condition for the symmetrized network guarantees synchronization in the asymmetrical network with node balance⁴⁹. According to the plot of the Lyapunov exponents versus the coupling strengths, all the Lyapunov exponents of the type I DAG are small. Therefore, we considered very long time simulation to estimate the Lyapunov exponents of type I DAG. The results confirm that the absolute values of all the exponents are smaller than $10^{-5}$ (see Fig. S2). Thus, the synchronous state for type I scale-free DAGs is globally stable but not asymptotically. This is similar to the synchronous stability of their building blocks, i.e. the FFL motifs (see Fig. 1). The type I inhibitory phase oscillators are also synchronized well on the DAG structures (see Fig. S3).

In order to better investigate the role of undirected topology of the DAGs on their synchronization, we construct the DAGs from scale-free networks ($p(d) \propto d^{-\gamma}$) with different scaling exponents ($\gamma$). As we know, increasing the scaling exponents leads to more degree-homogeneous networks. Figure 5 shows the synchronization of the DAGs constructed from the mentioned scale-free networks. The result shows that increasing the scaling exponents, decreases the clustering coefficient and synchronizability of scale-free DAGs. The underlying reason is that, by reducing the scaling exponent in networks with power-law degree distributions, the source node (the node with zero in-degree where the information flow starts from) would be able to reach every other node in the network easier. This characteristic has been previously shown for synchronization of directed networks of type II oscillators³⁰.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Synchronization stability of excitatory FBL and FFL motifs. (Top four rows) the Lyapunov exponents versus time (log-log scale) and the coupling for identical type I and type II excitatory phase oscillators connected by FBL (left column) and FFL loops (right column). The sign of Lyapunov exponents in the log-log plots are colour coded, where blue and red indicate negative and positive values, respectively. (Bottom row) The stationary order parameter versus coupling for identical type I and type II phase oscillators connected by FBL (left column) and FFL loops (right column).</figcaption></figure>

The effect of gradually adding feedback loops to the DAGs is also investigated. The results show that the synchronizability of the directed network decreases nonlinearly with increasing the number of the feedback loops (see Fig. S4).

In reality inhibitory and excitatory neurons work together to perform complex tasks. It has been suggested that excitatory and inhibitory inputs of a neuron are balanced, and this balance is important for the highly irregular firing observed in the cortex⁵⁰. Our results show that disconnected hybrid networks that are formed by feedforward motifs that their effective nodes (i.e. nodes which have outgoing edges) are purely inhibitory or excitatory, can provide balance networks with incoherent dynamics (see Fig. S5).

In the preceding parts, we investigated the synchronization of type I and type II phase oscillators, and we showed that type I excitatory and inhibitory oscillators which are connected by DAGs are synchronized. However, previous studies have shown that type I neurons are only synchronized with inhibitory synaptic connections⁴¹. To address this concern, we study the synchronizability of our oriented networks using conductance-based neuron models. First, we consider the type I Wang-Buzsáki inhibitory neurons connected by scale-free directed networks. The characteristics of the Wang-Buzsáki model neurons have been presented in the Supplementary Fig. S6. This model has been used to show that gamma rhythm (20–80 Hz), which is observed during behavioural arousal, can emerge in an undirected random network of interconnected inhibitory neurons⁴¹. However this rhythmic behaviour is only observed in the dense networks. Top panels of the Fig. 6 represent interspike distance synchrony level of the type I Wang-Buzsáki neurons connected via scale-free undirected (U), BDG, and DAG by inhibitory synapses in the $\Phi - I_{app}$ parameter
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Structural properties of two different directed networks. Adjacency matrices (top) and the location of Laplacian eigenvalues in the complex plane (bottom), for small-world and scale-free BDGs and DAGs with 200 nodes and average degree of 6. Schematic illustrations of BDGs and DAGs constructed from the same undirected backbone are shown at the right. The source node is highlighted with red colour.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Differences in the synchronization stability of excitatory scale-free networks with various directionalities and neuronal types. Lyapunov exponents versus coupling for (top row) identical type I and (middle row) type II phase oscillators situated on (left column) scale-free undirected graphs, (middle column) BDGs and (right column) DAGs, with $N = 200$ nodes are illustrated. The sign of the Lyapunov exponents is colour coded, where blue and red indicate negative and positive values, respectively. (Bottom row) The stationary order parameter versus the coupling strength of the scale-free networks with different link-directionalities (i.e. U(Undirected), BDG and DAG) are illustrated.</figcaption></figure>

spaces. The results verify that type I Wang-Buzsáki neurons are not synchronized on the sparse undirected graph and BDG. In contrast, there exists a parameter region in which neurons can be synchronized in the scale-free DAGs. Moreover, it has been shown that synaptic time constants potentially play a role in the determination of network frequency, and our results confirm this. In fact larger time constant makes the overall decay slower. Therefore, increasing $\tau_{syn}$ decreases the oscillation frequency for different values of the applied currents. In addition, we can see that the synchrony index displays a peak at the small values of the $\tau_{syn}$. This would mean that synchrony is high only within a restricted range of frequencies. These findings are consistent with the previous results obtained for dense undirected networks⁴¹. Here, the decay time constant of 10 ms that we used is consistent with $GABA_A$ receptor⁵¹. The parameter $\phi$ also affects the synchronization of the networks. In fact, decreasing $\phi$ reduces the firing rate by produc-
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of synchronization of type I excitatory phase oscillators on excitatory scale-free DAGs with different exponents. (First Column) The degree distribution of the scale-free networks with different scaling exponents, 1024 nodes and average degree of 10. Dashed lines represent linear curve fittings to the probability distribution functions. (Second Column) The stationary order parameters of scale-free DAGs versus scaling exponents, for $K = 3$. (Third and fourth Columns) display the average clustering coefficient (C) of the undirected backbones and the average shortest path from the source node to other nodes (S), versus the scaling exponents. The results are averaged over 50 realizations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Synchronization of Wang-Buzsáki inhibitory neurons coupled via scale-free undirected and directed graphs. (First three columns) Top panels display the interspike distance synchrony of the networks on $I_{app} - \phi$ phase space. Bottom panels represent the raster plots and the sum of the synaptic gates $\sum_i s_i$ for three network instances using the parameters specified by the arrows at the top plots. The synapse parameters are $\alpha = 12 \text{ ms}^{-1}$, $\beta = 0.1 \text{ ms}^{-1}$, $g_{syn} = 0.1 \text{ mS/cm}^2$. (Right column) Top and bottom panels show the effects of synaptic time constant ($\tau_{syn} = 1/\beta$) on synchrony and frequency of oscillation for different values of applied currents, respectively at $K = 1.25$. The parameters of the networks are valued as $N = 200$, and $\gamma = 3$.</figcaption></figure>

ing a deep afterhyperpolarization (see Fig. S7). In our simulation, a spike is detected when the voltage of membrane reaches a fixed threshold (–55 mV) from a lower value.

A similar behaviour is seen for the synchronization of the DAGs constructed by excitatory type I neurons by using Traub neuron model (see Figs. 7 and S8). However, for excitatory neurons the oscillation frequency increases as the time constant increases. The results on neuron models are also confirmed by using voltage synchrony measure (see the Figs. S9 and S10).

## Discussion

The observed complex spatiotemporal neural activity patterns reflect both the connectivity and dynamic properties of individual neurons. Two types of neuronal activities are found based on their different responses to the small depolarizing current pulses. Type I neurons only speed up the oscillation when they receive the stimulus, whereas type II neurons experience both phase advance and delay, depending on the timing of the perturbation. Previous studies have shown that, synchronizability of type II neurons is higher than type I neurons in excitatory networks, while inhibition is more stabilizing the synchrony of type I neurons<sup>52–56</sup>. Adaptation can stabilize the synchrony of excitatory type I neurons by changing their PRC<sup>14</sup>. But without adaptation, excitation can desynchronizing even for fully connected identical type I neurons beyond the weak coupling regime<sup>57</sup>. However, the master stability formalism suggest that fully connected network is an optimal structure for the synchronization of identical type II neurons. Therefore, the master stability approach can not be applied to investigate the synchronizability of the networks of type I neurons. On the other hand, neurons are connected to each other through synapses which are mostly
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Synchronization of Traub neurons coupled via scale-free undirected and oriented graphs. (First three columns) Top panels show the spike synchrony of the networks on $I_{app}$ – K phase space. Bottom panels exhibit the raster plots and the sum of the synaptic gates $\sum_{i} \varsigma_{i}$ for three network instances using the parameters specified by the arrows at the top plots. The synapse parameters are $\alpha = 12 \text{ ms}^{-1}$, $\beta = 0.5 \text{ ms}^{-1}$, $g_{syn} = 0.01 \text{ mS/cm}^2$. (Right column) Top and bottom panels represent the effects of synaptic time constant ($\tau_{syn} = 1/\beta$) on spike synchrony and frequency of oscillation, respectively at $K = 0.3$. The parameters of the network are valued as $N = 200$, and $\gamma = 3$.</figcaption></figure>

unidirectional in passing signals. Based on the observations mentioned above, we studied synchronization properties of directed networks of identical type I or type II neurons, assuming excitatory and inhibitory interactions.

Motifs are significantly over-represented subgraphs of a complex network, and play important roles in shaping its emergent behaviours. Therefore, firstly we focused on the synchronization properties of two different motifs in directed networks: feedback loops and feedforward loop motifs. Performing the linear stability analysis and the numerical simulation, we showed that synchrony can emerge in the feedforward loops with identical type I oscillators and purely excitatory or inhibitory connections. In fact, this synchronous state is globally stable but not asymptotically. Therefore, a small perturbation to the synchronous state would cause the system to converge to the fixed point milder than exponentially. However, the asynchronous state turns out to be asymptotically stable for feedback loops constructed from identical type I phase oscillators with purely excitatory or inhibitory connections. The analytical results obtained are compared and found to comply with the results obtained by using the extended Kuramoto model simulation. Next, we studied the synchronization of various types of neurons coupled via large directed networks. To this end, we used different methods for assigning the link directions to construct two different directed networks from the same undirected backbone, referred to as balanced degree graph (BDG) and directed acyclic graph (DAG). Precisely, many feedback loops are merged together to construct a BDG and only feedforward loop motifs are combined to make a DAG. Both the linear stability analysis and the numerical simulation confirmed that identical type I neurons connected by DAGs via strictly excitatory or inhibitory synapses are fully synchronized, while they are not synchronized when they are connected by BDGs. We can see that the steady Lyapunov exponents of type I DAGs are zero. This is in agreement with the result we found for the synchronization of feedforward loops. Therefore, studying the dynamics of the motifs can provide an intermediate step to better understand the synchronization of the larger collections. Besides the directionality, we showed that the undirected topology of the graphs also affects the synchronization of directed networks. The DAGs constructed from undirected networks with higher clustering coefficients, have higher number of feedforward motifs and higher synchronizability. In fact, a DAG is more synchronizable when all of its nodes are easily reachable from the source node. The results are also verified using conductance-based Wang-Buzsáki and Traub neuron models. In reality inhibitory and excitatory neurons work together to perform complex tasks. In general, excitatory and inhibitory inputs of a neuron are said to be balanced, and this balance is important for the highly irregular firing observed in the cortex⁵⁰. The global dynamic of motifs with hybrid effective nodes is very complex and directly related to the initial phase values. However, Our investigation suggests that excitatory-inhibitory balance and incoherent dynamics can be provided by merging inhibitory and excitatory feedforward motifs. Investigating the synchronization of neurons with different excitability types is important to understand the function of the diffuse modulatory cholinergic systems in the brain²⁰. Our main result is that edge directionality significantly alter the synchronizability of the type I neurons with purely excitatory connections, and it can not be ignored. Further studies are clearly required to determine the interplay between neuronal excitability, frequency, edge directionality, and network synchronization.

Received: 19 August 2019; Accepted: 31 January 2020;
Published online: 24 February 2020
References

1. Belykh, I., de Lange, E. & Hasler, M. Synchronization of bursting neurons: What matters in the network topology. *Phys. review letters* **94**, 188101 (2005).

2. Gregoriou, G. G., Gotts, S. J., Zhou, H. & Desimone, R. High-frequency, long-range coupling between prefrontal and visual cortex during attention. *Science* **324**, 1207–1210 (2009).

3. Wang, Q., Chen, G. & Perc, M. Synchronous bursts on scale-free neuronal networks with attractive and repulsive coupling. *PLoS one* **6**, e15851 (2011).

4. Glass, L. & Mackey, M. C. *From clocks to chaos: The rhythms of life* (Princeton University Press, 1988).

5. Motter, A. E., Myers, S. A., Anghel, M. & Nishikawa, T. Spontaneous synchrony in power-grid networks. *Nat. Phys.* **9**, 191 (2013).

6. Chowdhury, S. N., Majhi, S., Ozer, M., Ghosh, D. & Perc, M. Synchronization to extreme events in moving agents. *New J. Phys.* (2019).

7. Axmacher, N., Mormann, F., Fernández, G., Elger, C. E. & Fell, J. Memory formation by neuronal synchronization. *Brain research reviews* **52**, 170–182 (2006).

8. Missonnier, P. et al. Frontal theta event-related synchronization: comparison of directed attention and working memory load effects. *J. Neural Transm.* **113**, 1477–1486 (2006).

9. Tiitinen, H. T. et al. Selective attention enhances the auditory 40-hz transient response in humans. *Nature* **364**, 59 (1993).

10. Bartlett, N. R. & Bartlett, S. C. Synchronization of a motor response with an anticipated sensory event. *Psychol. review* **66**, 203 (1959).

11. Ermentrout, B. Type I membranes, phase resetting curves, and synchrony. *Neural computation* **8**, 979–1001 (1996).

12. Izhikevich, E. M. *Dynamical systems in neuroscience* (MIT press, 2007).

13. Schultheiss, N. W., Prinz, A. A. & Butera, R. J. *Phase response curves in neuroscience: theory, experiment, and analysis* (Springer Science and Business Media, 2011).

14. Ermentrout, B., Pascal, M. & Gutkin, B. The effects of spike frequency adaptation and negative feedback on the synchronization of neural oscillators. *Neural computation* **13**, 1285–1310 (2001).

15. Smeal, R. M., Ermentrout, G. B. & White, J. A. Phase-response curves and synchronized neural networks. *Philos. Transactions Royal Soc. B: Biol. Sci.* **365**, 2407–2422 (2010).

16. Bolhasani, E. & Valizadeh, A. Stabilizing synchrony by inhomogeneity. *Sci. reports* **5**, 13854 (2015).

17. Reid, M. A., Flores-Otero, J. & Davis, R. L. Firing patterns of type ii spiral ganglion neurons in *in vitro*. *J. Neurosci.* **24**, 733–742 (2004).

18. Stiefel, K. M., Gutkin, B. S. & Sejnowski, T. J. The effects of cholinergic neuromodulation on neuronal phase-response curves of modeled cortical neurons. *J. computational neuroscience* **26**, 289–301 (2009).

19. Stiefel, K. M., Gutkin, B. S. & Sejnowski, T. J. Cholinergic neuromodulation changes phase response curve shape and type in cortical pyramidal neurons. *PLoS one* **3**, e3947 (2008).

20. Bear, M. F., Corners, B. W. & Paradiso, M. A. *Neuroscience: exploring the brain* (Wolters Kluwer, 2016).

21. Platt, B. & Riedel, G. The cholinergic system, eeg and sleep. *Behav. brain research* **221**, 499–504 (2011).

22. Pecora, L. M. & Carroll, T. L. Master stability functions for synchronized coupled systems. *Phys. review letters* **80**, 2109 (1998).

23. Nishikawa, T. & Motter, A. E. Synchronization is optimal in nondiagonalizable networks. *Phys. Rev. E.* **73**, 065106 (2006).

24. Nishikawa, T. & Motter, A. E. Maximum performance at minimum cost in network synchronization. *Phys. D: Nonlinear Phenom.* **224**, 77–89 (2006).

25. D'Huys, O., Vicente, R., Erneux, T., Danckaert, J. & Fischer, I. Synchronization properties of network motifs: Influence of coupling delay and symmetry. *Chaos: An Interdiscip. J. Nonlinear Sci.* **18**, 037116 (2008).

26. Lodato, I., Boccaletti, S. & Latora, V. Synchronization properties of network motifs. *EPL (Europhysics Lett.* **78**, 28001 (2007).

27. Soriano, M. C., Van der Sande, G., Fischer, I. & Mirasso, C. R. Synchronization in simple network motifs with negligible correlation and mutual information measures. *Phys. review letters* **108**, 134101 (2012).

28. Esfahani, Z. G., Gollo, L. L. & Valizadeh, A. Stimulus-dependent synchronization in delayed-coupled neuronal networks. *Sci. reports* **6**, 23471 (2016).

29. Sporns, O. & Kötter, R. Motifs in brain networks. *PLoS biology* **2**, e369 (2004).

30. Sheshbolouki, A., Zarei, M. & Sarbazi-Azad, H. Are feedback loops destructive to synchronization? *EPL (Europhysics Lett.* **111**, 40010 (2015).

31. Gérard, C., Gonze, D. & Goldbeter, A. Effect of positive feedback loops on the robustness of oscillations in the network of cyclin-dependent kinases driving the mammalian cell cycle. *The FEBS journal* **279**, 3411–3431 (2012).

32. Mangan, S. & Alon, U. Structure and function of the feed-forward loop network motif. *Proc. Natl. Acad. Sci.* **100**, 11980–11985 (2003).

33. Alon, U. Network motifs: theory and experimental approaches. *Nat. Rev. Genet.* **8**, 450 (2007).

34. Bianconi, G., Gulbahce, N. & Motter, A. E. Local structure of directed networks. *Phys. review letters* **100**, 118701 (2008).

35. Zarei, M., SamaniAghababaei, K. & Omidi, G. R. Complex eigenvectors of network matrices give better insight into the community structure. *J. Stat. Mech. Theory Exp.* **2009**, P10018 (2009).

36. Son, S. W., Kim, B. J., Hong, H. & Jeong, H. Dynamics and directionality in complex networks. *Phys. review letters* **103**, 228702 (2009).

37. Sheshbolouki, A., Zarei, M. & Sarbazi-Azad, H. The role of leadership in the synchronization of directed complex networks. *J. Stat. Mech. Theory Exp.* **2015**, P10022 (2015).

38. Mato, G. & Samengo, I. Type i and type ii neuron models are selectively driven by differential stimulus features. *Neural computation* **20**, 2418–2440 (2008).

39. Kuramoto, Y. *Chemical oscillations, waves, and turbulence*, vol. 83 (Springer-Verlag Berlin Heidelberg, 1984).

40. Mofakham, S., Fink, C. G., Booth, V. & Zochowski, M. R. Interplay between excitability type and distributions of neuronal connectivity determines neuronal network synchronization. *Phys. Rev. E.* **94**, 042427 (2016).

41. Wang, X. J. & Buzsáki, G. Gamma oscillation by synaptic inhibition in a hippocampal interneuronal network model. *J. neuroscience* **16**, 6402–6413 (1996).

42. Traub, R. Simulation of intrinsic bursting in ca3 hippocampal neurons. *Neuroscience* **7**, 1233–1242 (1982).

43. Tiesinga, P. H. & Sejnowski, T. J. Rapid temporal modulation of synchrony by competition in cortical interneuron networks. *Neural computation* **16**, 251–275 (2004).

44. Lim, W. & Kim, S. Y. Coupling-induced spiking coherence in coupled subthreshold neurons. *Int. J. Mod. Phys. B.* **23**, 2149–2157 (2009).

45. Wolf, A., Swift, J. B., Swinney, H. L. & Vastano, J. A. Determining lyapunov exponents from a time series. *Phys. D: Nonlinear Phenom.* **16**, 285–317 (1985).

46. Watts, D. J. & Strogatz, S. H. Collective dynamics of small-world networks. *Nature* **393**, 440 (1998).

47. Barabási, A. & Albert, R. Emergence of scaling in random networks. *Science* **286**, 509–512 (1999).

48. Caughman, J. S. & Veerman, J. Kernels of directed graph laplacians. *Electronic journal combinatorics* **13**, 39 (2006).

49. Belykh, I., Belykh, V. & Hasler, M. Synchronization in asymmetrically coupled networks with node balance. *Chaos: An Interdiscip. J. Nonlinear Sci.* **16**, 015102 (2006).

50. Haider, B., Duque, A., Hasenstaub, A. R. & McCormick, D. A. Neocortical network activity *in vivo* is generated through a dynamic balance of excitation and inhibition. *J. Neurosci.* **26**, 4535–4545 (2006).
51. Börgers, C. *An introduction to modeling neuronal dynamics*, vol. 66 (Springer, 2017).

52. Whittington, M. A., Traub, R. D., Kopell, N., Ermentrout, B. & Buhl, E. H. Inhibition-based rhythms: experimental and mathematical observations on network dynamics. *Int. journal psychophysiology* **38**, 315–336 (2000).

53. Bartos, M. et al. Fast synaptic inhibition promotes synchronized gamma oscillations in hippocampal interneuron networks. *Proc. Natl. Acad. Sci.* **99**, 13222–13227 (2002).

54. Chow, C. C., White, J. A., Ritt, J. & Kopell, N. Frequency control in synchronized networks of inhibitory neurons. *J. computational neuroscience* **5**, 407–420 (1998).

55. Brunel, N. & Hansel, D. How noise affects the synchronization properties of recurrent networks of inhibitory neurons. *Neural Comput.* **18**, 1066–1110 (2006).

56. Rich, S., Booth, V. & Zochowski, M. Intrinsic cellular properties and connectivity density determine variable clustering patterns in randomly connected inhibitory neural networks. *Front. neural circuits* **10**, 82 (2016).

57. Hansel, D., Mato, G. & Meunier, C. Synchrony in excitatory neural networks. *Neural computation* **7**, 307–337 (1995).

**Acknowledgements**

M.Z. wishes to acknowledge financial support of “The Abdus Salam International Centre for Theoretical Physics (ICTP)”.

**Author contributions**

M.Z. conceived the study, A.Z. and A.S. jointly designed and implemented the simulation model, M.Z. and A.Z. created the analytic results, M.Z. wrote the paper and A.S. edited the manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

Supplementary information is available for this paper at https://doi.org/10.1038/s41598-020-60205-0.

Correspondence and requests for materials should be addressed to M.Z.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020