{'fraction_non_alphanumeric': 0.04590086915668311, 'fraction_numerical': 0.009302325581395349, 'mean_word_length': 4.283196823033011, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
In-phase oscillations from the cooperation of cellular and network positive feedback in synaptically-coupled oscillators

We study the emergent dynamics of a network of synaptically coupled slow-fast oscillators. Synaptic coupling provides a network-level positive feedback mechanism that cooperates with cellular-level positive feedback to ignite in-phase network oscillations. Using analytical bifurcation analysis, we prove that the Perron-Frobenius eigenvector of the network adjacency matrix fully controls the oscillation pattern locally in a neighborhood of a Hopf bifurcation. Besides shifting the focus from the spectral properties of the network Laplacian matrix to the network adjacency matrix, we discuss other key differences between synaptic and diffusive coupling.


I. INTRODUCTION
Synchronization is usually studied in the context of diffusive coupling [1], [2], i.e., when the interaction between the oscillators is proportional to the difference in their states. Focusing on diffusive coupling has various limitations. First, because diffusive coupling is passive, oscillators must be intrinsic, that is, they must exhibit limit cycle oscillations in the absence of network interactions. Second, the only type of emergent network activity is a practically synchronous one, where, for large enough diffusive coupling strength, the oscillators converge to the same state modulo a synchronization error. Motivated by understanding the emergence of sustained in-phase oscillations in the suprachiasmatic nucleus (SCN) in the mammal master circadian clock [3], [4], we introduce a model of slow-fast oscillators with synaptic-like coupling and study the emergence of in-phase oscillations in it. A fundamental experimental observation, reproduced in our model but impossible to reproduce in diffusively coupĺed models, is that in the SCN many clock neurons behave as sustained oscillators only in the presence of network interactions, whereas they behave as damped oscillators when isolated [5]. SCN dynamics are therefore emergent, in the sense that the collective behavior (sustained oscillations) relies on network interactions and it is distinctively different from the isolated node behavior (damped oscillations) [6]. 1 The intrinsic dynamics of our oscillators include a saturated fast cellular positive feedback loop and a linear slow negative feedback loop. It is a simplified version of excitable neural dynamics [7], [8]. The interaction of the two loops leads to relaxation (slow-fast) neural-like oscillations for strong enough positive feedback through a Hopf bifurcation.
The network synaptic-like couplings are approximated as saturated inputs to the receiving oscillator of the state of the sending oscillator. They provide network positive feedback, which cooperates with cellular positive feedback to ignite and shape emergent network oscillations.
The contributions of our analysis are the following. First, we prove a general lemma for the spectral properties of a class of block-defined matrices with the structure of the Jacobian matrix of our model. Second, we show that diffusive coupling cannot induce synchronous oscillations in a network of damped oscillators. Third, we prove that under a strongly directed network topology synaptic coupling can lead to inphase oscillations even when the uncoupled oscillators are damped. If the coupling is in-regular, in-phase oscillations become synchronous, i.e., all oscillators converge to the same state. In this work, we rely on (local) bifurcation analysis at the model equilibrium and show that the (dominant) Perron-Frobenius eigenvector of the network adjacency matrix fully determines the in-phase oscillation pattern. Our results are in line with existing ones on automata synchronization [9], [10] and, together with [11], they stress the importance of considering non-diffusive coupling in synchronization studies. In future works, we will couple our local results with a global analysis, using, for instance, dominance analysis [12]. Also, we only consider here homogeneous (identical) intrinsic dynamics. In future works we will relax this assumption as well by exploiting the power of synaptic coupling of being naturally apt to cope with non-synchronous in-phase oscillations, as those that are expected in heterogeneous populations.

II. NOTATION AND DEFINITIONS
N denotes the set of positive natural numbers, and R the set of real numbers. In general, N ∈ N will be a positive integer. As usual, Re(z) = x denotes the real part of a complex number z = x + iy ∈ C. R N denotes the set of real N -tuples, and x ∈ R N denotes an arbitrary N -tuple. Because of the specific models used, it will be convenient to denote R 2N = R N ×R N and its elements as (x, y) ∈ R 2N . The zero and one vectors 0 N ∈ R N , 1 N ∈ N, denote tuples which have all their entries equal to zeroes and ones, respectively. Finally, a vector is said to be positive if all its entries are strictly positive, denoted by x > 0.
The set M N ×N contains all real N × N matrices represented as M = (M ij ). I N denotes the identity matrix in dimension N , and O N denotes the zero matrix in dimension N . The determinant of a matrix A ∈ M N ×N is denoted by |A|, and its characteristic polynomial is denoted by p(λ) = |A − λI N |.
A matrix is said to be simple if all of its diagonal entries are equal to zero.
A weighted digraph G = (V, A) is a 2-tuple consisting of a set of vertices or nodes V = {1, . . . , N } and an adjacency matrix A ∈ M N ×N with the convention that there exists a directed edge from vertex j to vertex i if and only if A ij = 0, in which case A ij is the weight of the edge. We will always assume that A ii = 0, i.e., there are no self-loops in the digraph, so that every adjacency matrix considered is simple as defined before. Given a node i ∈ V of a weighted digraph G , its weighted in-degree is denoted by Throughout this paper we don't require the graphs to be undirected, that is, we don't assume that A is symmetric. Therefore, the eigenvalues µ 1 , . . . , µ N of the Laplacian matrix L − may be complex and they all satisfy Under such a condition, d − will denote the global in-degree of the weighted digraph. G is strongly connected if, for any two nodes i and j, there exists a directed path which connects i to j; in this case, its adjacency matrix is said to be irreducible.

III. A NETWORK OF SLOW-FAST OSCILLATORS
We present a single, general model that includes diffusive and excitatory synaptic coupling between slow-fast damped or sustained oscillatorṡ is the time constant of the slow variables y i , α i > 0 are cellular positive feedback gains, and S is a sigmoid function modeling intrinsic and synaptic nonlinearities. A d is the diffusive coupĺing adjacency matrix, and A e is the excitatory coupling adjacency matrix. Clearly, (x 0 , y 0 ) = (0 N , 0 N ) constitutes an equilibrium point for system (1). When A e = O N the coupling between the oscillators is purely diffusive, and when A d = O N the coupling is purely excitatory. The two matrices A d and A e define the diffusive G d (V, A d ) and excitatory G e (V, A e ) digraphs, respectively.

A. Transitions from damped to sustained oscillations ruled by a Hopf bifurcation
We now show that the parameter α rules the transition from damped to sustained slow-fast oscillations for uncoupled oscillators. Consider model (1) for N = 1, which reduces to the single-oscillator model The Jacobian matrix evaluated at equilibrium is readily For both eigenvalues are purely imaginary. Moreover, continuity of the discriminant function ∆ guarantees λ 1,2 ∈ C\R for α sufficiently close to α H . Furthermore, observe that B. In-phase oscillations from network positive feedback between two coupled oscillators It is easy to show that the model undergoes a network Hopf bifurcation along the parametric curve √ β 1 β 2 = 1 + ε, at which point the oscillators start to oscillate in phase, as shown in Fig. 1. Observe that the uncoupled oscillators are damped in this case. The positive feedback brought by network interactions has the double role of both igniting and synchronizing the emergent oscillations.
We will show that the behavior observed in this lowdimensional example is impossible in general if the coupling is diffusive, whereas it generalizes to arbitrary strongly connected excitatory coupling topologies.

IV. A USEFUL LEMMA
During our discussion, we will find several block-wise defined matrices of the form where M ∈ M N ×N is any real matrix, ε > 0 is a (small) real constant, and a, c, d ∈ R. The following general lemma will turn out very useful in our analysis.  (3). Then its characteristic polynomial p(λ), for λ = −ε, is obtained as Moreover, any eigenvector (x, y) ∈ R 2N of J, corresponding to an eigenvalue λ ∈ C\{−ε}, must satisfy Proof. Obtaining the characteristic polynomial only requires us to apply the determinant formula for block-wise defined matrices [14] to This will require −(ε+λ)I N to be invertible, which imposes λ = −ε. Thus, the formula calculates p(λ) as , whence (4) follows. As for the eigenvector condition (5), consider λ ∈ C\{−ε} a root of p(λ) as given above, and This yields the linear system Given ε + λ = 0, we may solve for y in the second equation as y = 1 ε+λ (εI N + dM )x. We then substitute this into our first equation, getting (−(a + λ)I N + cM )x − 1 ε+λ (εI N + dM )x = 0 N . From here the second eigenvector condition follows, thus concluding this proof.
The relevance of this lemma lies in that the matrix M fully characterizes the spectral properties of the higherdimensional matrix J. More precisely, equation (4) establishes a one-to-two correspondence between the eigenvalues of J and those of M . Equation (5)  Proof. Observe that Thus, the model Jacobian computed at equilibrium is given by where L − is the in-Laplacian matrix associated to G d (V, A d ), which is exactly in the form of Lemma 1, with matrices J = J d , M = L − , and parameters a = 1 − α, c = −1, d = 0. The associated characteristic polynomial reads Let µ 1 , . . . , µ N be the eigenvalues of L − and recall that for k ∈ {1, . . . , N }. Setting ε = 0 in (6) yields λ − k = 0 with multiplicity m = N , and λ + k = α − 1 − µ k , which satisfies Re(λ + k ) = α − 1 − Re(µ k ) < 0. By continuity, the real parts of eigenvalues λ + k remains negative for sufficiently small ε > 0. To guarantee a similar result for λ − k , we can split (6) into its real and imaginary parts. Letting λ = σ + iτ and µ = u + iv, we get which can be interpreted as zero-level sets of some functions F , G, respectively. Using the Implicit Function Theorem, variables σ and τ can be expressed as functions S, T of the remaining variables ε, u, v whenever is satisfied. The derivative ∂S ∂ε is readily obtained by implicit differentiation as which is negative for ε = σ = τ = 0, given that u = Re(µ) 0 and α < 1. recall that λ − k = 0 is obtained as a zero of (6) when ε = 0. Thus, for positive and sufficiently small values of ε, every J d -eigenvalue has negative real part, and the equilibrium at the origin is locally exponentially stable.
Theorem 1 shows that diffusive coupling requires intrinsic oscillators to lead to synchronous network oscillations. One could provide a global proof by means of Lyapunov functions and convergent systems analysis [15]. We omit this proof due to space constraints.

VI. NETWORK AND CELLULAR POSITIVE FEEDBACK COOPERATE IN TRIGGERING SYNCHRONOUS OSCILLATIONS IN NETWORKS OF SLOW-FAST DAMPED NODES
We now turn to the network positive feedback present in model (1), for A d = O N and non-negative A e . Throughout this section we will make the standing homogeneity assumption α i = α for every i ∈ {1, . . . , N }, i.e., we assume that the uncoupled oscillators are identical. We will relax this homogeneity assumption in future works. We also let A e = βA, where A is a simple matrix. The two parameters α > 0, β > 0 govern cellular and network positive feedback, respectively. Then, the Jacobian of model (1) evaluated at its equilibrium at the origin reads So we may apply Lemma 1 to matrix J = J e , considering M = A, a = 1 − α, c = β, d = 0, to arrive at the following result.
Lemma 2. Let A ∈ M N ×N be a simple, non-negative matrix, and consider model (1) with A d = O N , A e = βA and α i = α, where α 0 and β > 0 are non-negative parameters. Let J e be the Jacobian matrix of this system evaluated at equilibrium (0 N , 0 N ). Then any eigenvector (x, y) of J e , associated to eigenvalue λ, must satisfy the conditions Moreover, by letting µ 1 , . . . , µ N be the eigenvalues of matrix βA, we obtain for each µ k two J e -eigenvalues λ ± k , where Proof. Apply Lemma 1 with J = J e , M = βA, a = 1 − α, c = 1, d = 0. Then equation (8) holds. This implies that the eigenvalues of βA and J e are linked by the expression From here we obtain two J e -eigenvalues λ ± k , which are indeed given by expression (9), thus ending the proof.
The relationship established at the end of Lemma 1 has become clearer, in that equation (9) explicitly determines J eeigenvalues as functions of βA-eigenvalues (and, therefore, of A-eigenvalues). Thus, spectral analysis of matrix J e , i.e. local analysis of purely excitatory system (1) under homogeneity hypotheses, reduces to spectral analysis of adjacency matrix A.

A. In-regular homogeneous network
We start by showing that if the coupling topology is in-regular, then model (1) undergoes a Hopf-bifurcation for strong enough cellular and network positive feedback. Furthermore, because 1 N is the dominant eigenvector of the adjacency matrix, the Hopf bifurcation happens along the synchronization space where each oscillator has the same state. That is, the network Hopf bifurcation leads to synchronous network oscillations.
Theorem 2. Let A ∈ M N ×N be an irreducible, simple, non-negative matrix associated to a strongly connected inregular digraph of global in-degree d − > 0, and consider model (1) with A d = O N , A e = βA and α i = α, where α ∈ [0, 1) and β > 0. Then, for sufficiently small ε > 0 the system undergoes a Hopf bifurcation along the parametric curve β = 1+ε−α d − . Moreover, the center manifold associated to the bifurcation is tangent to the synchronization subspace and is locally exponentially stable.
Proof. Given A e = βA, where A is in-regular, we conclude that µ 1 = βd − > 0 is an eigenvalue with corresponding eigenvector x 1 = 1 N . Applying formula (9) to this eigenvalue yields two J e -eigenvalues, namely λ ± 1 , given by the expression Thus for β = 1+ε−α part. Indeed, irreducibility of matrix A makes it possible to apply the Perron-Frobenius Theorem which guarantees that the dominant eigenvalue d − > 0 has algebraic and geometric multiplicity one. Transversality is also easily verified (we omit details here due to space constraint), which yields the Hopf bifurcation. Conditions (8) give us the eigenvectors ε). Now, identifying the real and imaginary parts u = (1 N , ε1 N ), v = (0 N , ε(1 − ε)1 N ) of the spanning vectors, it follows that J e u ∓ iJ e v = J e (u ∓ iv) = ±iλ(u ∓ iv) = λv ± iλu. J e , u, v and λ are real (matrices, vectors and values), so this last equation implies J e (u) = λv and J e v = −λu. Therefore the center manifold is tangent to the span of vectors u, v, whence the form of subspace E is obtained. Because all other eigenvalues have negative real part, the associated center is locally exponentially attractive.
Theorem 2 shows that when β − 1+ε−α d − > 0 the model exhibits synchronous oscillations. This condition can be fulfilled both by increasing the cellular positive feedback α for fixed network positive feedback β or vice-versa. Figure  2 numerically illustrates the predictions of our theorem.

B. Strongly-connected homogeneous network
In-regular networks are too restrictive to accurately model biological networks like the SCN. In this section we relax the in-regularity assumption. Irreducibility of the adjacency matrix for strongly connected coupling topologies implies uniqueness of a Perron-Frobenius eigenvector that fully determines the pattern of in-phase oscillations emerging at the network Hopf bifurcation.
Theorem 3. Let A ∈ M N ×N be an irreducible, simple, nonnegative matrix associated to a strongly connected digraph, and consider model (1) associated to A d = O N , A e = βA and α i = α, α ∈ [0, 1) and β > 0. Let ρ > 0 be the leading eigenvalue of A, and x 0 > 0 the Perron eigenvector associated to ρ. Then, for ε > 0 and sufficiently small, the system undergoes a Hopf bifurcation along the parametric curve β = 1+ε−α ρ . Moreover, the center manifold associated to the bifurcation is tangent to the real subspace and is locally exponentially stable.
Proof. Given A e = βA, where A is irreducible and nonnegative, we conclude that µ 1 = βρ > 0 is its leading real eigenvalue with corresponding eigenvector x 0 > 0. Applying formula (9) to this eigenvalue yields two J e -eigenvalues λ ± 1 , where Thus, for β = 1+ε−α ρ > 0, λ ± 1 are purely imaginary complex eigenvalues while all other eigenvalues have negative real part. Indeed, irreducibility of matrix A makes it possible to apply the Perron-Frobenius Theorem which guarantees that leading eigenvalue ρ > 0 has algebraic and geometric multiplicity one. Transversality is once again easily verified, which yields the Hopf bifurcation. Setting λ = ε(1 − ε), conditions (8) give us the J e -eigenvector z 1,2 = (x 0 , (ε ∓ i ε(1 − ε))x 0 ) associated to λ ± 1 = ±iλ at bifurcation. As in the previous Theorem, writing z 0 = u 0 ∓iv 0 in its real and imaginary parts, one again sees that u 0 and v 0 span the real tangent subspace to the center manifold, whence we conclude the form of subspace E. Because all other eigenvalues have negative real part at the bifurcation, the associated center manifold is locally exponentially stable.
Theorem 3 shows that when β = 1+ε−α ρ > 0 the model exhibits in-phase oscillations. This condition can be fulfilled both by increasing the cellular positive feedback α for fixed network positive feedback β or vice-versa. Furthermore, it shows that the Perron-Frobenius eigenvector of the adjacency matrix fully controls the oscillation pattern, at least close to the Hopf bifurcation. Figure 3 numerically illustrates the predictions of our theorem.

A. Model extension
As it is the case for many physiological networks, sometimes outputs from one node affect the receiving node along multiple timescales. To incorporate such an effect in our model one could consider an extended case of the excitatory version, adjusting equations (1) to account for the effect of x-variables over y-variables. Under homogeneity conditions, computations of this extended case would be very similar to those made before, as seen by applying Lemma 1 when d = 0.

B. Extension to global results
The analysis made in this paper relies solely on local properties of dynamical systems near equilibrium. Therefore, a more complete and formal approach should also incorporate x 0 x(t) , where ·, · and · denote the standard scalar product and 2-norm, respectively, and x 0 is the adjacency matrix Perron-Frobenious eigenvector as defined in Theorem 3. Observe that l(t) → 1 at almost all time points, i.e., excluding time points where x(t) = 0. In other words, along the in-phase network oscillations, the x-component of the state vector is parallel to the Perron eigenvector x 0 , as predicted by Theorem 3. global tools before and after bifurcation to guarantee convergence to either a stable steady state or a stable limit cycle. For example, in the diffusive case one could propose a Lyapunov function [17] or use the theory of convergent systems [15]. Alternatively one can use dominance analysis [12], through which it might be possible to show the existence of a globally attractive and invariant 2-dimensional manifold corresponding to the center manifold of the Hopf bifurcation, which would effectively make our local bifurcation analysis global.

C. Heterogeneous populations
It is evident that real life networks won't maintain, in general, the homogeneous hypothesis which were used extensively in the proofs of the excitatory case. A more delicate analysis should be provided when considering heterogeneous networks, much more likely to be found in real life phenomena, through higher dimensional bifurcation theory [13].

D. Application to circadian rhythmogenesis
The model in this work was originally motivated by the synchronization phenomena observed in the suprachiasmatic nucleus (SCN) of the mammal hypothalamus. One may identify different subpopulations and connections inside the SCN (e.g. spatially [18], GABAergic [19], neuropeptidergic [4], [20]). Neuromodulation here not only affects the electrophysiological rhythms, but also gives input to the molecular clock through a much slower loop [21]. Other works have pointed out that neural appositions and density of connections vary according from one neuropeptidergic subpopulation to the other [22], [23]. Therefore, a multilayer digraph model may result useful in capturing the dynamic properties of this circadian phenomenon. Additional topologies A e j could be incorporated to account for different neuropeptide release (VIP, AVP, GRP being the main ones).