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

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 1-3 , cardiac pacemaker cells 4 , and power grid networks 5 to social networks 6 . 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 7 , directed attention 8,9 , and the processing of sensory stimuli 10 .
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 11,12 .
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 13 . 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

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 bellow.
It has been shown that, the number of feedback loops increases by enhancing the correlation between the inand out-degrees of the nodes 34 . 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 35 .
Directed acyclic graphs (DAGs) are constructed by extending the residual degree gradient method 36,37 . 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 1 ij , if nodes i and j are connected, and a ij = 0, otherwise. Whereas, the adjacency matrix of a directed graph is defined as, = a 2 ij , if a link is directed from node j to node i, and = a 0 ij , 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 15 . It can be defined as, Where T 0 is the unperturbed cycle period of the neuron, and θ T is the cycle period when the neuron is perturbed at phase θ. 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 www.nature.com/scientificreports www.nature.com/scientificreports/ versus θ diagram, contains both positive and negative regions. Studies verify that the Hodgkin-Huxley neurons have PRC II excitability type, while the Wang-Buzsáki and Traub neurons have PRC I excitability type 38 . 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-Buzsá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 39,40 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: where θ 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 0 i and = u 1 i , respectively. The function θ G( ) corresponds to the phase response curve of the oscillator (see the Fig. S1). In our simulations, the initial values of θ i are randomly drawn from a uniform distribution in the interval π [0, 2 ], 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 Here, … represents averaging over different network realizations and initial conditions. The magnitude is ≤ ≤ r 0 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-Buzsáki model to analyse the dynamics of the networks with type I inhibitory interneurons. The Wang-Buzsáki neuron model includes the following differential equation 41 : ), 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 (α h , β h , α n , β n ). s, α and β are activation variable and rate functions for synapses. φ 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 Here, θ 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: 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 The parameters of the Wang-Buzsáki model 41 are summarized in the  Table S1. The rate functions are given by the following functions:  I  I  I  I  I   m   dv  dt  Na  K  L  syn a pp , where gating variables m n , and h obey equations of the type . The rate functions are: , and synaptic current to the ith neuron is evaluated from 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 : where, N is the number of neurons in the network and τ = − denotes the interspike interval of the merged set of network spikes. Here, the inter-spike intervals are calculated between different neurons. … τ 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 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 Scientific RepoRtS | (2020) 10:3306 | https://doi.org/10.1038/s41598-020-60205-0 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 49 . According to the plot of the Lyapunov exponents versus the coupling strengths, all the Lyapunov exponets 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 ( ) with different scaling exponents (γ). 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 synchronizabilty 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 30 . www.nature.com/scientificreports www.nature.com/scientificreports/ 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 50 . 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 41 . 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 , which is observed during behavioural arousal, can emerge in an undirected random network of interconnected inhibitory neurons 41 . 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 Φ − I app parameter 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 τ 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 τ 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 41 . Here, the decay time constant of 10 ms that we used is consistent with GABA A receptor 51 . The parameter φ also affects the synchronization of the networks. In fact, decreasing φ reduces the firing rate by produc-  www.nature.com/scientificreports www.nature.com/scientificreports/ 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, synchonizability of type II neurons is higher than type I neurons in excitatory networks, while inhibition is more stabilizing the synchrony of type I neurons [52][53][54][55][56] . Adaptation can stabilize the synchrony of excitatory type I neurons by changing their PRC 14 . But without adaptation, excitation is desynchronizing even for fully connected identical type I neurons beyond the weak coupling regime 57 . 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  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 50 . 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 20 . 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.