Structural constraints on the emergence of oscillations in
  multi-population neural networks
Oscillations arise in many real-world systems and are associated with both
functional and dysfunctional states. Whether a network can oscillate can be
estimated if we know the strength of interaction between nodes. But in
real-world networks (in particular in biological networks) it is usually not
possible to know the exact connection weights. Therefore, it is important to
determine the structural properties of a network necessary to generate
oscillations. Here, we provide a proof that uses dynamical system theory to
prove that an odd number of inhibitory nodes and strong enough connections are
necessary to generate oscillations in a single cycle threshold-linear network.
We illustrate these analytical results in a biologically plausible network with
either firing-rate based or spiking neurons. Our work provides structural
properties necessary to generate oscillations in a network. We use this
knowledge to reconcile recent experimental findings about oscillations in basal
ganglia with classical findings.
Jie Zang, Shenquan Liu, Pascal Helson, Arvind Kumar
# Structural constraints on the emergence of oscillations in multi-population neural networks

###### Abstract

Oscillations arise in many real-world systems and are associated with both functional and dysfunctional states. Therefore, it is important to determine the causes of oscillations in a network. Whether a network can oscillate can be estimated if we know the strength of interaction between nodes. But in real-world networks (in particular in biological networks) it is usually not possible to know the exact connection weights. Therefore, it is important to determine the structural properties of a network necessary to generate oscillations. Here, we provide a proof that uses dynamical system theory to prove that an odd number of inhibitory nodes and strong enough connections are necessary to generate oscillations in a single cycle threshold-linear network. We illustrate these analytical results in a biologically plausible network with either firing-rate based or spiking neurons. Our work provides structural properties necessary to generate oscillations in a network. We use this knowledge to reconcile recent experimental findings about oscillations in basal ganglia with classical findings.
Introduction

Oscillations are ubiquitous in dynamical systems [1; 2]. They have important functional consequences but can also cause system malfunction. In the brain for instance, oscillations take part in information transfer [3; 4]. However, persistent beta band (13-30 Hz) oscillations are associated with the pathological symptoms of Parkinson's disease [5]. Therefore, it is important to determine when and how a system of many interacting nodes (network) oscillates.

This question is usually very difficult to answer analytically. The main tool that can be used is the Poincare-Bendixson theorem [6; 7] which is only valid in 2 dimensions, which drastically reduce its applicability. In some cases when we know the model parameters it is possible to calculate whether the system will oscillate or not. However, often such parameters cannot be measured experimentally. For example, in most physical, chemical, and biological networks, it is usually not possible to get the correct value of connectivity strength. By contrast, it is much easier to know whether two nodes in a system are physically connected and what is the sign (positive or negative) of their interactions. Therefore, it is much more useful to identify necessary structural conditions for the emergence of oscillations. A good example is the conjecture postulated by Thomas [8]: when considering a coupled dynamical system (\(\dot{x}=f(x)\) and \(x(0)\in\mathbb{R}^{n}\)) with a Jacobian matrix that has elements of fixed sign, it can exhibit oscillations only if the directed graph obtained from the nodes' connectivity (Jacobian matrix) admits a negative loop of two or more nodes (loop with an odd number of inhibition). This conjecture has been proven using graph theory for smooth functions \(f\)[9; 10].

Thomas also conjectured that the assumption on the constant sign of the Jacobian matrix may not be necessary [11], i.e. having a negative loop in some domain of the phase space should be sufficient to generate oscillations. This condition is more realistic due to the ubiquity of the non-linearity in biological systems. For example, in the brain, even though neurons are (usually) either excitatory or inhibitory, the transfer function linking the neurons is non-linear and can thus lead to elements of the Jacobian matrix with non-constant sign. To the best of our knowledge, this last conjecture has not been proved yet but there are many examples of it. For instance, oscillation can emerge from a simple EI network Wilson-Cowan model [12]. Our study is an example of this conjecture on the threshold-linear
network (TLN) model [13] which can closely capture the neural population dynamics.

Here, we study the long term behaviour of the TLN model in the case of a single cycle interaction containing all nodes. We show analytically that regardless of the sign of this loop, the system cannot oscillate when connections are too weak as the system possesses a unique globally asymptotically stable fixed point. However, when connections are strong enough, the system either possesses two asymptotically stable fixed points (positive loop) or a unique unstable fixed point (negative loop). In addition, the system can be shown to be bounded and thus, it has one of the following long term behavior: limit cycle, quasi-periodic or chaotic behavior. Interestingly, we can show that such dynamics can be shut down by introducing positive external input to excited nodes.

Based on our analytical results, we used simulations of basal ganglia (BG) network models with either firing rate-based or spiking neurons to explain recent experimental findings about the origin of oscillations in Parkinson's disease (PD). Traditionally, the subthalamic nucleus and globus pallidus (STN-GPe) subnetwork is considered to be the key network underlying the emergence of oscillations in PD [14; 15; 16]. However, recent experiments have shown that near complete inhibition of GPe but not of STN is sufficient to quench oscillations [17]. This observation contradicts several previous models and even clinical observations in which surgical removal of STN is used to alleviate PD symptoms. Our theory suggests that there are at least 6 possible cycles in the Cortex-BG network that have the potential to oscillate based on the connectivity structure. We show that even if STN is inhibited, other 'cycles' can sustain pathological oscillations. Interestingly, we found that GPe features in 5 out of 6 oscillatory cycles and therefore GPe inhibition is likely to affect PD-related oscillations in most cases.

## II Results

We study how the emergence of oscillations in a network of excitatory and inhibitory populations depends on the connectivity structure. We first consider a network of nodes with dynamics representing the average firing rate of a population. We derive structural conditions for the emergence of oscillations when the dynamics of individual nodes are described according to the threshold-linear network (TLN) model. Next, we use numerical simulations to test whether such results might still hold on two other models: the Wilson-Cowan popu
lation rate-based model [18] and a network model of the basal ganglia with spiking neurons (see Methods).

### Structural conditions to generate oscillations

#### Intuition behind the analytical results

There exist many ways to generate oscillations in a network. Oscillation can arise from individual nodes due to their intrinsic dynamics (a spiking neuron can have a periodic behaviour given the ionic channel composition [19]) or from the weights' dynamics when considering synaptic plasticity [20]. Here we assume that the system's ability to oscillate only depends on the connectivity structure: the presence of positive or negative loops and the connections strength (Jacobian matrix) within them. That is, neither plasticity nor biophysics of neurons is considered.

Consider a small network of two nodes. If we connect them mutually with excitatory synapses, intuitively we can say that the two-population network will not oscillate. Instead, the two populations will synchronize. The degree of synchrony will, of course, depend on the external input and the strength of mutual connections. If both these nodes are inhibitory, one of the nodes will emerge as a winner and the other will be suppressed [21]. Hence, a network of two mutually connected inhibitory populations cannot oscillate either. We can extend this argument to three population networks with three connections that form a closed loop or 'cycle' (Fig. 1a, top). When all three connections in the cycle are excitatory, the three populations will synchronize. Essentially, we will have a single population. Thus, these two and three population motifs are not capable of oscillations.

The simplest network motif which is capable of oscillating consists of two mutually connected nodes: one excitatory and one inhibitory (EI motif: Fig. 1a, bottom) [12]. When there are three populations connected with three connections to form a cycle, the potential to oscillate depends on the number of inhibitory connections. A cycle with one inhibitory connection (EEI motif) can be effectively reduced to an EI motif and therefore, can oscillate. However, when there are two inhibitory connections (EII motif, Fig. 1a, top)), the two inhibitory neurons engage in a winner-take-all type dynamics and the network is not capable of oscillations. Finally, if there are three inhibitory connections (i.e. all three nodes are
inhibitory, III motif) the network enters in a winner-less-competition [22] and can exhibit oscillations (Fig. 1a, bottom).

These examples of two or three nodes suggest that a network can generate oscillations if there are one or three inhibitory connections in the network. These observations form the basis for the conjecture of Thomas [8] that gives a necessary condition for oscillations to emerge. This condition is of course not sufficient. In the following we find additional constraints (input and minimum connection strength) needed to determine the emergence of oscillations in a network. To this end we use the TLN model which captures the neural population dynamics to a great extent. After proving the key theorems, we test with simulation whether similar results hold on a more realistic Wilson-Cowan model and a model of basal ganglia with spiking neurons.

#### Threshold linear network model

We consider the TLN\((W,b)\) in which individual nodes follow the dynamics

\[\frac{dx_{i}}{dt}=-x_{i}+\left[\sum_{j=1}^{n}W_{ij}x_{j}+b_{i}\right]_{+},\quad i =1,\ldots,n \tag{1}\]

where \(n\) is the number of nodes, \(x_{i}(t)\) is the activity level of the \(i\)th node at time \(t\geq 0\), \(W_{ij}\) is the connection strength from node \(j\) to node \(i\) and \([\,\cdot\,]_{+}\stackrel{{\rm def}}{{=}}\max\{\,\cdot\,,0\}\) is the threshold non-linearity. For all \(i\in[n]\stackrel{{\rm def}}{{=}}\{1,\ldots,n\}\), the external inputs \(b_{i}\in\mathbb{R}\) are assumed to be constant in time. We refer to a \(n\) neurons network with dynamics given by eq. 1 as TLN\((W,b)\).

In order to help the definition of cycle connectivity matrices, we define

\[C_{n}\stackrel{{\rm def}}{{=}}\{(i,j)\in[n]^{2}|i-j=1\}\cup\{(1,n )\}.\]

We denote by \(\delta_{i,I}\) the Kronecker delta which equals 1 when node \(i\) is inhibitory and 0 otherwise (node \(i\) is excitatory). In the following, we use the convention that node 0 is node \(n\) and node \(n+1\) is node 1. For a given set of elements \(\{y_{k}\}_{k\in\mathbb{N}}\) in \(\mathbb{R}\), we will use the convention:

\[\prod_{k=i}^{j}y_{k}=1\mbox{ when }j<i. \tag{2}\]

We define \(\mathcal{A}=\{a_{1},\ldots,a_{n_{I}}\}\) as the ensemble of inhibited nodes (\(\{k\in[n]|\delta_{k-1,I}=1\}\)) put in order such that \(a_{1}<\cdots<a_{n_{I}}\). Denoting by card\((\cdot)\) the cardinal function, we have that card\((a)=n_{I}\). We also use the cycle convention for \(\mathcal{A}\): \(a_{n_{I}+1}=a_{1}\).
#### Analytical results

**Theorem 1**.: _Let a network of inhibitory and excitatory nodes be connected through a graph \(G\) which does not contain any directed cycle. Assume that its nodes follow TLN(\(W,b\)) dynamics

Figure 1: **Structural condition for oscillations: odd inhibitory cycle rule and its illustrations.****a**: Examples of oscillating motifs and non-oscillating motifs in Wilson-Cowan model. Motifs that cannot oscillate show features of Winner-take-all: the winner will inhibit other nodes with a high activity level. Inversely, the oscillatory ones all show features of winner-less competition, which may contribute to oscillation. **b**: The odd inhibitory cycle rule for oscillation prediction with the sign condition of a network. **c**: Illustrations of oscillation in complex networks. Based on the odd inhibitory cycle rule, Network I can’t oscillate, while Network II could oscillate by calculating the sums of their motifs. The red or black arrows indicate inhibition or excitation, respectively. Hollow nodes and solid nodes represent excitatory and inhibitory nodes, respectively.
(eq. 1) with_

\[W_{ij}=\begin{cases}w_{ij}(-1)^{\delta_{j,I}}&\text{when edge $i\gets j\in G$}\\ 0&\text{otherwise},\end{cases}\]

_where \(w_{ij}\in\mathbb{R}^{+}\ \forall\ i,j\in[n]\). Then, TLN(\(W,b\)) has a unique globally asymptotically stable fixed point._

**Theorem 2**.: _Let G be a cyclical graph with \(n_{I}\in\mathbb{N}^{+}\) inhibitory nodes and \(n_{E}\in\mathbb{N}\) excitatory nodes such that \(n_{I}+n_{E}\geq 2\) (\(\geq 3\) when \(n_{I}=1\)). Assume that the nodes follow the TLN(\(W,b\)) dynamics (eq. 1) with for all \(i,j\in[n]\), \(w_{j}\in\mathbb{R}^{+}\),_

\[W_{ij}=\begin{cases}w_{j}(-1)^{\delta_{j,I}}&\text{when $(i,j)\in C_{n}$}\\ 0&\text{otherwise},\end{cases}\]

_and \(b_{i}=0\) when the node \(i-1\) is excitatory and \(b_{i}>0\) otherwise. Moreover, using convention (eq. 2), assume that the initial state is bounded,_

\[\forall j\in\{x_{a_{k}},\ldots,x_{a_{k+1}-1}\},\quad x_{j}(0)\in[0,b_{a_{k}} \prod_{i=a_{k}}^{j-1}w_{i}]. \tag{3}\]

_Then, the long time behaviour of the network depends on the following conditions,_

\[\forall k\in[n_{I}],\quad\prod_{i=a_{k}}^{a_{k+1}-1}w_{i}<\frac{b _{a_{k+1}}}{b_{a_{k}}}, \tag{4}\] \[\prod_{i=a_{k}}^{a_{k+1}-1}w_{i}>\frac{b_{a_{k+1}}}{b_{a_{k}}},\] (5) \[\sqrt[n]{\prod_{i=1}^{n}w_{i}}<\frac{1}{cos(\pi/n)},\] (6) \[\sqrt[n]{\prod_{i=1}^{n}w_{i}}>\frac{1}{cos(\pi/n)}. \tag{7}\]

_If \(n_{I}\) is even and_

* _eq._ 4 _is satisfied, TLN(_\(W,b\)_) has a unique globally asymptotically stable fixed point with support_ \([n]\)
* _eq._ 5 _is satisfied, TLN_\((W,b)\) _has two asymptotically stable fixed points with strict complementary subsets of_ \([n]\) _as supports._

_If \(n_{I}\) is odd and_

* _eq._ 4 _is satisfied, TLN_\((W,b)\) _has a unique fixed point which is globally asymptotically stable and its support is_ \([n]\)_,_
* _eq._ 5 & eq._ 6 _is satisfied, TLN_\((W,b)\) _has a unique fixed point which is asymptotically stable (not globally) and its support is_ \([n]\)_,_
* _eq._ 5 & eq._ 7 _are satisfied, TLN_\((W,b)\) _has a unique fixed point which is unstable and has_ \([n]\) _as support._

_Remark 1_.: First, note that eq. 4 implies

\[\sqrt[n]{\prod_{i=1}^{n}w_{i}}<1, \tag{8}\]

and similarly, eq. 5 implies

\[\sqrt[n]{\prod_{i=1}^{n}w_{i}}>1. \tag{9}\]

In addition, the bound on the initial state eq. 3 can be easily removed. We use it because it eases the proof as we then don't need to introduce technical details that are not interesting for this study.

Then, Theorem 2 says that a possible condition for the one cycle TLN to oscillate is that the number of inhibitory nodes is odd when the connection strength are strong enough (i.e. eq. 5 & eq. 7). In that case, the system has no stable fixed point and from Lemma 1 it is bounded so it has one of limit cycle, quasi-periodic or chaotic behaviours. In particular, Theorem 2 states that the odd number of inhibitory nodes is not sufficient. Indeed, when eq. 4 holds and \(n_{I}\) is odd, no oscillations are possible as the fixed point is globally stable. It is also the case when \(n_{I}\) is even which corresponds to Thomas' conjecture.

Finally, there is a gap in between the conditions (between eq. 4 and eq. 5 for example) for which the long term behaviour is not determined.
_Remark 2_.: In particular, if for all \(i\in[n]\), \(w_{i}=w\in\mathbb{R}_{+}^{*}\) and for all \(k\in[n_{i}]\), \(b_{a_{k}}=b\in\mathbb{R}_{+}^{*}\), then the dynamics of the system only depends on \(w\). When \(n_{I}\) is even: \(w<1\) implies that \(\operatorname{TLN}(W,b)\) have a unique globally asymptotically stable fixed point; \(w>1\) implies that the fixed point for \(w<1\) becomes unstable and \(\operatorname{TLN}(W,b)\) has two more asymptotically stable fixed point. If \(n_{I}\) is odd, \(\operatorname{TLN}(W,b)\) only has a unique fixed point which is asymptotically stable when \(w<\frac{1}{cos(\pi/n)}\) (globally when \(w<1\)) and unstable when \(w>\frac{1}{cos(\pi/n)}\).

_Remark 3_.: In Theorem 2, we assume that the external inputs are absent for excited nodes. Assume that the external input to any excited node, say node \(a_{k}<i<a_{k+1}\), is strictly positive. Then, bounding its dynamics as in Lemma 1, we know that its activity will be more than \(b_{i}\). Hence, the next inhibited node \(a_{k+1}\) can be silenced forever if

\[b_{i}\prod_{j=i}^{a_{k+1}-1}w_{j}>b_{a_{k+1}},\]

thus destroying the cycle structure and thus preventing oscillation from emerging.

On the other hand, if the external inputs to excited nodes are strictly negative, Theorem 2 conclusion will be similar but now with condition described in eq. 4 replaced by

\[b_{a_{k}}\prod_{i=a_{k}}^{a_{k+1}-1}w_{i}-\sum_{j=a_{k}+1}^{a_{k+1}-1}b_{j} \prod_{i=j}^{a_{k+1}-1}w_{i}<b_{a_{k+1}}.\]

This means that cycles with even (odd) inhibitory nodes need strong enough connections to generate multi-stability (limit cycle). We now clarify that the latter condition relates to weights' strength. With \(w=(w_{1},\cdots,w_{n})\) and using the set of functions increasing functions \((f_{i})_{1\leq i\leq n}\) such that

\[f_{i}^{w,b}(x)=w_{i-1}x_{i-1}-b_{i}\]

one can write the last condition as

\[f_{a_{k+1}}^{w,b}\circ\cdots\circ f_{a_{k}}^{w,b}(b_{a_{k}})<0.\]

Hence, the left term is increasing with any weight strength.

One should also note that under this negative input assumption to excited nodes, when weights are weak, the support of the fixed point might be different from \([n]\). In particular, some excited nodes might not belong to the support.
_Remark 4_.: When the decay rates are not the same (here all of them are \(-1\)), similar results hold but then the conditions for stability are more difficult to state precisely. Finally, when considering the EI network (two nodes), the system always admits to a unique globally asymptotically stable fixed point with support \(\{1,2\}\). Indeed, it is easy to show that the system will always reach the domain where the inhibitory node is small enough so that one can remove the threshold function in eq. 1 and thus the eigenvalues of the Jacobian matrix are \(\pm i\sqrt{w_{1}w_{2}}-1\). No oscillations are then possible, which is an easy example to show that negative loops are not sufficient to generate oscillation in non smooth dynamical systems.

Similar results have been shown by Snoussi [9] and Gouze [10]. Considering dynamical systems of the form \(\dot{x}=f(x)\) where \(f\) is a continuously differentiable function on a given open convex set and \(f\) has a constant sign Jacobian matrix, they used graph theory methods to show that negative loop in this matrix is a necessary condition to generate oscillations. In our case, \(f\) is not continuously differentiable, the Jacobian matrix elements can change sign within the state space and we show that there is a need of additional constraints for oscillations to arise. A formal proof of the aforementioned theorems is provided in Appendix by using classical dynamical theory tools.

### Intuition behind the proof of the theorems

The idea behind our proof can be explained graphically. We assume that nodes cannot oscillate due to their intrinsic activity and a fixed external input only drives them to a non-zero activity which does not change over time. Therefore, they need input from their pre-synaptic (upstream) nodes to change their state in a periodic manner to generate oscillations. In such a network, if we perturb the node \(i\) with a pulse-like input, it is necessary that the perturbation travels through the network and returns to the node \(i\) with a \(180^{o}\) phase shift (i.e. with an inverted sign). Otherwise, the perturbation dies out and each node returns to a state imposed by its external input.

In a network without directed cycles, it is possible to sort the nodes into smaller groups where nodes do not connect to each other (Fig. 2a). That is, a network with no directed cycles, can be rendered as a feed-forward network in which the network response by definition does not return to the node (or group) that was perturbed. Such a network can only oscillate when the intrinsic dynamics of individual nodes allow for oscillatory dynamics.
However, having a directed cycle is no guarantee of oscillations because network activity must return to the starting node with a 180\({}^{\circ}\) phase shift. This requirement puts a constraint on the number of inhibitory connections in the cycles. When we assume that there are no delays (or the delay is constant) in the connections, excitatory connections do not introduce any phase shift, however, inhibitory connections shift the phase by 180\({}^{\circ}\) (in the simplest case invert the sign of the perturbation). Given this, when a cycle has an even number of inhibitory connections the cycle cannot exhibit oscillations (Fig. 2b,top). However, replacing an inhibitory connection by an excitatory one can render this cycle with an ability to oscillate (Fig. 2b, bottom). Therefore, odd number of inhibitory appears to be necessary for oscillation to emerge.

### The effect of network parameters on oscillations

To test the validity of our theorems in more realistic biological neuronal networks, we numerically simulated the dynamics of the Wilson-Cowan model. Specifically, we investigated the role of synaptic transmission delays, synaptic weights, external inputs and self-connection in shaping the oscillations when the network has directed cycles. In particular, we focused on two networks: the III motif with three inhibitory nodes (odd inhibitory links) and the EII motif with one excitatory and two inhibitory nodes (even inhibitory links).

Our numerical simulation showed that for a wide range of parameters (synaptic delays, synaptic weights, external input and self-inhibition), while III network showed oscillations, EII network did not show any oscillations (Fig. 3 and S1). The oscillation frequency however depended on the exact value of the synaptic delays, synaptic weights and external inputs. For instance, increasing the synaptic delay reduced the oscillation frequency (Fig. 3 b). Synaptic delays play a more important role in shaping the oscillations in an EI type network (see Supplementary Fig. S2). The effect of increasing the synaptic strength was contingent on the external inputs. In general, increasing the synaptic strength resulted in a reduction in the oscillation frequency (Fig. 3 c). Next, oscillation frequency changed in a non-monotonic fashion as a function of external input irrespective of the choice of other parameters (Fig. 3 d). Typically, a mid-range input strength resulted in maximum oscillation frequency. Finally, increasing the self-connection of nodes increased the oscillation frequency but beyond a certain self-connection the node was completely silenced and it changed the network topology
Figure 2: **The intuitive explanations of Theorem 1 and 2.****a**: A visual representation of why directed cycles are important in network oscillation. By rearranging all nodes, any network without directed cycles can be seen as a feed-forward network which will make the system reach a stable fixed point. **b**: An intuitive explanation of the odd inhibitory cycle rule by showing the activities of two 6-node-loops. Odd inhibitory connections (bottom) can help the system oscillate, while even inhibitory connections has the opposite effect.
and oscillations disappeared (Fig. 3 e).

Overall, these results are consistent with our rule that odd number of inhibitory nodes and strong enough connections are necessary to induce oscillations in a directed cycle. The actual frequency of oscillations depends on specific network parameters.

## Oscillators in the cortex-basal ganglia network

Next, we use our theorem to explain recent experimental observations about the mechanisms underlying the emergence of oscillations in the basal ganglia. Emergence of 15-30 Hz (beta band) oscillations in the cortico-basal ganglia (CBG) network is an ubiquitous feature of Parkinson's disease (PD) [23; 24; 25; 26]. Based on their connectivity and activity subthalamic nucleus (STN) and the globus pallidus externa (GPe) subnetwork has emerged as the most likely generator of beta oscillations [27; 28]. The STN-GPe subnetwork becomes oscillatory when their mutual connectivity is altered [15; 29] or neurons become bursty [30; 31] or striatal inputs to GPe increase [16; 32; 33; 34]. However, oscillations might also be generated by the striatum [35], by the interaction between the direct and hyperdirect pathways [36] and even by cortical networks that project to the BG [37]. Recently, de la Crompte et al. [17] used optogenetic manipulations to shed light on the mechanisms underlying oscillation generation in PD. They showed that GPe is essential to generate beta band oscillations while motor cortex and STN are not. These experiments force us to rethink the mechanisms by which beta band oscillations are generated in the CBG network.

To better understand when GPe and/or STN are essential for beta band oscillations, we identified the network motifs which fulfill the odd inhibitory cycle rule. For this analysis, we excluded D1 SPNs because they have a very low firing rate in the PD condition [33]. In addition, cortex is assumed as a single node in the CBG network.

The CBG network can be partitioned into 238 subnetworks with 2, 3, 4, 5 or 6 nodes (see Supplementary Fig. S3-S8). Among these partitions, there are five loops (or cycles) in the CBG network with one or three inhibitory projections: Proto-STN, STN-GPi-Th-cortex, Proto-Arky-D2, Proto-FSN-D2, and Proto-GPi-Th-Cortex-D2 (Fig. 4a). One or more of these 5 loops appeared in 88 (out of 238) subnetworks of CBG (see Fig. 4b, colors indicate different loops). Larger subnetworks consisting of 5 and 6 nodes have multiple smaller subnetworks (with 2 or 3 nodes) that can generate oscillations (boxes with multiple colors
Figure 3: **Influence of network properties on the oscillation frequency in motifs III and EII with Wilson-Cowan model.****a**: The changed network parameters are shown in the table. Red (green) connections are inhibitory (excitatory) and black arrows are the external inputs. **b-e**: We systematically varied the synaptic delay time **b**, synaptic weights **c**, external input **d**, and self-connection **e**. These parameters were varied simultaneously for all the synapses i.e. in each simulation all synapses were homogeneous. Green, orange, red and turquoise respectively show the effect of synaptic delay, synaptic strength, external input and self-inhibition. See the Supplementary Fig. S1 and Fig. S2 for more detailed results about III and EI network motifs.
in Fig. 4b).

Based on our odd inhibitory cycles in BG, we found three oscillatory subnetworks which do not involve the STN (Fig. 4a, cyan, green and purple subnetworks). However, each of these oscillatory subnetworks involves Prototypical neurons (from the GPe) which receive excitatory input from STN. Therefore, it is not clear whether inhibition of STN can affect oscillations or not. To address this question, we first simulated the dynamics of a four-node motif (Fig. 4c top) using the Wilson-Cowan type model (see Methods). In this subnetwork, we have three cycles: Proto-STN loop with one inhibitory connection, Proto-STN-Arky-D2 loop with three inhibitory connections and Proto-Arky-D2 with three inhibitory connections.

We systematically varied external inputs to the STN and D2-SPNs and measured the frequency of oscillations (see Methods). We found that for weak inputs to the D2-SPNs, the Proto-STN subnetwork generated oscillations for weak positive input (Fig. 4c bottom). However, as the input to D2-SPNs increased, the oscillation frequency decreased and oscillations were observed even for very strong drive to STN (Fig. 4c bottom). That is, in this model, both Proto-STN and Proto-D2-Arky subnetworks compete for oscillations, which subnetwork wins depending on their inputs. To disentangle the oscillations of each of these two subnetworks, we performed 'lesion' experiments in our model (see Methods). These experiments also mimicked lesions performed in non-human primates [30].

When we removed the D2-SPN to Proto projections, the network could oscillate but only because of the Proto-STN subnetwork (Fig. 4d). In this setting, we get relatively high-frequency beta band oscillations but only for a small range of excitatory inputs to the STN (Fig. 4d, bottom). In this setting, inhibition of STN would certainly abolish any oscillation. Next, we removed the STN output (equivalent to inhibition of STN), the Proto-D2-Arky subnetwork generated oscillations for weak positive inputs to the D2-SPNs (Fig. 4d, bottom). Note that unlike in Fig. 4c, here we injected additional input to Proto to compensate for the loss of excitatory input from STN and to ensure that it had sufficient baseline activity. The frequency of Proto-D2-Arky oscillations was smaller than that observed for the Proto-STN subnetwork because the former involves a three synapses loop. However, as we have shown earlier, frequency of oscillation can be changed by scaling the connection weights or external inputs (Fig. 3). Overall, these results suggest that, in principle, it is possible for CBG network to oscillate even when STN is removed from the network.
Figure 4: **Schematic of CBG network model with potential oscillators and the interaction between two oscillators in Wilson-Cowan model.****a**: CBG structure with red lines denoting inhibition and green lines denoting excitation, along with five potential oscillators based on the odd inhibitory cycle rule. **b**: Oscillation in all BG motifs from 2 nodes to 6 nodes based on the odd inhibitory cycle rule. Each grid represents a separate motif. We use different colors to mark motifs that can oscillate, and each color means an oscillator from panel **a**. **c**: The reaction of oscillation frequency to different external inputs to D2 and STN in a BG subnetwork. External inputs to Proto and Arky are 1 and 3, respectively. **d**: Same thing as **c** but ruining the connection from D2 to Proto. **e**: Same thing as **c** but destroying the connections from STN and increasing the input to Proto from 1 to 4.
#### Oscillations in model of basal ganglia with spiking neurons

Thus far we have only illustrated the validity of our theorems in a firing rate-based model. To be of any practical value to brain science, it is important to check whether our theorems can also help in a network with spiking neurons. To this end, we simulated the two subnetworks with 3 inhibitory connections: Proto-D2-FSN and Proto-D2-Arky (see Methods). These subnetworks were simulated using a previous model of BG with spiking neurons [34].

The subnetworks Proto, Arky, D2-SPN, FSN have very little recurrent connectivity to oscillate on their own. We provided Poisson type external input. All neurons in a subnetwork received the same input rate but a different realization of the Poisson process. Both Proto-D2-FSN (Fig. 5a) and Proto-D2-Arky (Fig. 5b) subnetworks showed \(\beta\)-band oscillations. In the loop Proto-D2-FSN, D2-SPN neurons have a relatively high firing rate. This could be a criterion to exclude this loop as a potential contributor to the beta oscillations.

Next, we mimicked the STN inhibition experiments performed by de la Crompe et al. [17] in our model. To this end, we simulated the dynamics of BG network excluding D1-SPNs (because of their low firing rate in PD condition) and FSN (because with FSNs in the oscillation loop, D2-SPNs may have non-physiological firing rates). In this reduced model of BG, we changed inputs to operate in a mode where either Proto-D2-Arky (Fig. 5c) or Proto-STN (Fig. 5d) loop was generating the oscillations. In both cases, we systematically increased the inhibition of STN neurons.

In the Proto-D2-Arky mode, as we inhibited STN neurons, firing rate of the Proto neurons decreased and oscillations in STN population diminished but Proto neurons showed clear beta band oscillations (Fig. 5c). By contrast, and as expected when STN-Proto loop was generating the oscillations, increasing the STN inhibition abolished the oscillations in both STN and Proto neurons (Fig. 5d). In STN-Proto loop, when STN is inhibited, there is no cycle left in the network and therefore oscillations diminished, whereas the Proto-D2-Arky loop remains unaffected by the STN inhibition (except for a change in the firing rate of the Proto neurons). As shown in Fig. 4c, whether oscillations are generated by the Proto-D2-Arky or STN-Proto loop depends on the relative input to the D2 or STN neurons. So it is possible that in rodents, D2-SPN have stronger input from the cortex than STN and therefore, oscillations survive despite near complete inhibition of STN.
## Discussion

Here we prove in a single cycle TLN model and illustrate with numerical simulations of biological networks that when the number of inhibitory nodes in a directed cycle is odd

Figure 5: **Oscillations in a leaky integrate-and-fire (LIF) spiking neuronal network model of specific BG motifs.****a-b**: Average peristimulus time histograms (PSTH) of all neurons in **a** Proto-FSN-D2 and **b** Proto-Arky-D2 motifs under Parkinson condition with power spectral density (PSD) at the top right. **c**: PSTH of Proto and STN in a BG subnetwork with motif Proto-Arky-D2 as the oscillator during different STN inhibition. **d**: Same thing as **c** but changing the oscillator from Proto-Arky-D2 to Proto-STN.
and connections are strong enough, then the system has the potential to oscillate. In 1981, Thomas [8] conjectured that at least one negative feedback loop (i.e., a loop with an odd number of repressors) is needed for gene regulatory networks to have periodic oscillating behavior. This conjecture was proven for Boolean dynamical systems by Snoussi [9] and Gouze [10]. But their proof required that node transfer-function is differentiable everywhere. We here prove a more complete theorem for a case where node transfer-function is threshold-linear as is the case for many network in the brain. Thus, together with previous results of Snoussi [9] and Gouze [10] we further expand the scope within which we can comment on the potential of a network to generate oscillation based only on the connectivity structure alone. In addition, we complete this condition by one on weights' strength stating that the latter needs to be strong enough for the system to possibly oscillate. Eventually, oscillations can be quenched by adding positive external input to excited nodes.

A key assumption of our analysis is that there are no delays in the network. Indeed, delays within and between subnetwork connections can have a big effect on the oscillations [38]. In the numerical simulations of basal ganglia network, we included biologically realistic synaptic delays (i.e. connection delays were shorter than the time constants of the neurons). Our results suggest that such delays do not influence our results and they only determine the oscillation frequency. But it is not possible to comment on how the results may change when delays become longer than the time constant of the node.

#### Interactions between input and network structure

Previous models suggest that when we excite the excitatory node or inhibit the inhibitory node oscillations can emerge and strengthen [12, 16]. By contrast, when we inhibit the excitatory node or excite the inhibitory node, oscillations are quenched. This can be summarised as the 'Oscillations Sign Rule'. Let us label the excitatory population as positive and inhibitory as negative. Let us also label excitatory inputs as positive and inhibitory inputs as negative. Now if we multiple the sign of the node and sign of the stimulation, we can comment on the fate of oscillations in a qualitative manner. For example, inhibition of inhibitory nodes would be \(-\times-=+\) i.e. oscillations should be increased and when we inhibit excitatory nodes, it would be \(-\times+=-\) i.e. oscillations should be decreased. The 'Oscillation sign rule' scales to larger network with more nodes. With the 'Odd Cycle Rule'
as we have shown we can comment on whether a directed cycle will oscillate or not from the count of inhibitory links. When we combine the 'Oscillations Sign Rule' with the 'Odd Cycle Rule' we can get a more complete qualitative picture of whether a stimulating a node in a network will generate oscillations or not.

#### Interaction between node properties and network structure

In our proof we have assumed that nodes follow rather simple dynamics and have a threshold-linear transfer-function. In reality nodes in physical, chemical and biological systems can have more complex dynamics. For instance, biological neurons have the property of spike frequency adaptation or rebound spiking. Similarly, synapses in the brain can increase or decrease their weights based on the recent history of inputs which is referred to as the short-term-facilitation or short-term-depression of synapses [39]. Such biological properties can be absorbed in the network structure in the form of an extra inhibitory or excitatory connection. When nodes can oscillate given their intrinsic dynamics then the question becomes more about whether a network structure can propagate oscillations to other nodes.

#### Oscillations in the basal ganglia

We applied our results to understand the mechanisms underlying the emergence of PD-related pathological oscillations in the basal ganglia. Given that there are 8 key neuron populations in the basal ganglia, we enumerated 238 possible directed cycles. From 2-node-motifs to 6-node-motifs, our odd cycle rule identified 88 potential directed cycles that can generate oscillations. Among these, 81 cycles feature GPe (either Proto or Arky type or both) and 66 feature STN. Which specific cycle underlies oscillations depends on the exact input structure. For instance, when input to STN is higher than the D2 neurons, the STN-GPe network generates oscillations. But when inputs to D2 neurons are stronger, the D2-Proto-Arky cycle can become the oscillator. That is, STN is not necessary to generate oscillations in the basal ganglia. Our results also suggest that besides focusing on the network connectivity, we should also estimate the inputs to different nodes in order to pinpoint the key nodes underlying the PD-related pathological oscillations - that would be the way to reconcile the recent findings of de la Crompe [17] with previous results.
#### Beyond neural networks

In this work we have used the odd cycle rule to study oscillations in basal ganglia. However, oscillatory dynamics and the odd cycle rule show up in many chemical, biological and even social systems such as neuronal networks [40], psychological networks [41], social and political networks [42, 43, 44, 45], resting-state networks in autism [46] and gene networks [47, 48]. In fact, originally Thomas' conjecture [8] about the structural conditions for oscillations was made for gene regulatory networks. Therefore, we think that insights obtained from our analytical work can be extended to many other chemical, biological and social networks. It would be interesting to check to what extent our prediction of quenching oscillation by exciting the excitatory nodes holds in other systems besides biological neuronal networks.

## Methods

To study the emergence of oscillations in the basal ganglia, we used three models: Threshold-Linear Network (TLN), Wilson-Cowan model and network with spiking neurons. TLN model (eq. 1) was used here to rigorously prove that simple conditions, such as the odd inhibitory cycle rule, can lead to oscillations (Theorems 1 and 2). Wilson-Cowan type firing rate-based model was used to find the structural constraints on oscillations and to determine the effect of network properties (such as delays, synaptic weights, external inputs, and self-inhibition) on the emergence of oscillations. Finally, to demonstrate the validity of the odd inhibitory cycle rule in a more realistic model, we use a network with spiking neurons.

### Wilson-Cowan dynamics

In the firing rate-based models, we reduced each cortex-basal ganglia (CBG) subnetwork to a single node. To describe firing rate dynamics of such a node, we used the classic Wilson-Cowan model [18]

\[\tau\frac{dr_{i}(t)}{dt}=-r_{i}(t)+F\left(\sum_{j=1}^{n}w_{ij}r_{j}+I_{i}^{ext }\right) \tag{10}\]

where \(r_{i}(t)\) is the firing rate of the \(i\)th node, \(\tau\) is the time constant of the population activity, \(n\) is the number of nodes (or subnetworks), \(w_{ij}\) is the strength of connection from node
to \(i\), and \(I_{i}^{ext}\) is the external input to the population. \(F\) is a nonlinear activation function relating output firing rate to input, given by

\[F(x)=\frac{1}{1+e^{-a(x-\theta)}}-\frac{1}{1+e^{a\theta}} \tag{11}\]

where the parameter \(\theta\) is the position of the inflection point of the sigmoid, and \(\frac{a}{4}\) is the slope at \(\theta\). Here, \(\tau\), \(\theta\), and \(a\) are set as 20, 1.5, and 3. Other parameters of the model varied with each simulation. The simulation specific parameters for Fig. 3 are shown in Tables 1, 2 and for Fig. 4 are shown in Table 3, respectively.

#### Network model with spiking neurons

The basal ganglia network with spiking neurons was taken from a previous model by Chakravarty et al. [34]. Here we describe the model briefly and for details we refer the reader to the paper by Chakravarty et al. [34].

#### Spiking neuron model

Here, we excluded D1-SPNs because they have a rather small firing rate in PD conditions. The striatal D2-type spiny neurons (D2-SPN), fast-spiking neurons (FSNs) and STN neurons were modelled as standard LIF neurons with conductance-based synapses. The membrane potential \(V^{x}(t)\) of these neurons was given by:

\[C_{m}^{x}\frac{dV^{x}(t)}{dt}=I_{e}(t)+I_{syn}(t)-g_{L}^{x}\left[V^{x}(t)-V_{ reset}^{x}\right] \tag{12}\]

where \(x\in\{\)D2-SPN, FSN, and STN\(\}\), \(I_{e}(t)\) is the external current induced by Poisson type spiking inputs (see below), \(I_{syn}(t)\) is the total synaptic input (including both excitatory and inhibitory inputs). When \(V^{x}\) reached the threshold potential \(V_{th}^{x}\), the neuron was clamped to \(V_{reset}^{x}\) for a refractory duration \(t_{ref}\) = 2 ms. All the parameter values and their meaning for D2-SPN, FSN and STN are summarized in Tables 2, 3, 4, respectively.

We used the LIF model with exponential adaptation (AdEx) to simulate Proto and Arky neurons of the globus pallidus externa (GPe), with their dynamics defined as

\[C^{x}\frac{V^{x}(t)}{dt} =-g_{L}^{x}[V^{x}(t)-V_{\rm reset}^{x}]-w^{x}+I_{\rm syn}^{x}(t)+ I_{e}+g_{L}^{x}\Delta_{T}\exp\left(\frac{V^{x}(t)-V_{T}^{x}}{\Delta_{T}}\right)\] \[\tau_{w}\dot{w}^{x} =a\left(V^{x}(t)-V_{\rm reset}^{x}\right)-w^{x}\]
where \(x\in\{\)Proto, Arky\(\}\). Here when \(V^{x}(t)\) reaches the threshold potential (\(V_{th}^{x}\)), a spike is generated and \(V^{x}(t)\) as well as \(w^{x}\) will be reset as \(V_{\text{reset}}^{x}\), \(w^{x}+b\), respectively, where b denotes the spike-triggered adaptation. The parameter values and their meaning for Proto and Arky are specified in Table S5. Neurons were connected by static conductance-based synapses. The transient of each incoming synaptic current is given by:

\[I_{\text{syn}}^{x}\left(t\right)=g_{\text{syn}}^{x}\left(t\right)\left[V^{x}( t)-E_{rev}^{x}\right]\]

where \(x\in\{\)D2-SPN, FSN, STN, Arky, and Proto\(\}\). \(E_{rev}^{x}\) is the synaptic reversal potential and \(g_{\text{syn}}^{x}(t)\) is the time course of the conductance transient, given as follows:

\[g_{\text{syn}}^{x}(t)=\left\{\begin{array}{l}J_{syn}^{x}\frac{t}{\tau_{syn }}\exp\left(\frac{-(t-\tau_{syn})}{\tau_{syn}}\right),\text{ for }t\geq 0\\ 0,\text{ for t}<0\end{array}\right.,\]

where \(syn\in\{\)exc, inh\(\}\), \(J_{syn}^{x}\) is the peak of the conductance transient and \(\tau_{syn}^{x}\) is synaptic time constant. The synaptic parameters are shown in Table S6.

Some of model parameters were changed to operate the BG model in specific modes dominated by a 2 or 3 nodes cycle. The Table S7 and Table S8 show the parameters of Fig. 5c and Fig. 5d, respectively.

### External input

Each neuron in each sub-network of the BG received external input in the form of excitatory Poisson-type spike trains. This input was provided to achieve a physiological level of spiking activity in the network. For more details please see Chakravarthy et al. [34]. Briefly, the external input was modelled as injection of Poisson spike-train for a brief period of time by using the inhomogeneous_poisson_generator device in NEST. The strength of input stimulation can be controlled by varying the amplitude of the EPSP from the injected spike train.

### STN inhibition experiment

We set a subnetwork of basal ganglia to study how STN inhibition affects oscillation when different motifs dominate the system. The connections and external inputs to each neuron in Fig. 5c and 5d are shown in Tables S7 and S8. To simulate the increasing inhibition to
STN, the external input to STN was reduced from 1 pA to -99 pA in Fig. 5c and from 30 pA to -50 pA in Fig. 5d.

### Data analysis

The estimate of oscillation frequency of the firing rate-based model was done using the power spectral density calculated by pwelch function of MATLAB. The spiking activity of all the neurons in a sub-population were pooled and binned (rectangular bins, bin width = 0.1 ms). The spectrum of spiking activity was then calculated for the binned activity using pwelch function of MATLAB.

### Simulation tools

Wilson-Cowan type firing rate-based model was simulated using Matlab. All the relevant differential equations were integrated using Euler method with a time step of 0.01 ms. The network of spiking neurons was simulated in Python 3.0 with the simulator NEST 2.20 [49]. During the simulation, differential equations of BG neurons were integrated using Runga-Kutta method with a time step of 0.1 ms.

### Code availability

The simulation code will be made available on GitHub upon publication of the manuscript.
\begin{table}
\begin{tabular}{l c c c c c}  & \multicolumn{3}{c}{Synaptic weights} & \multicolumn{2}{c}{population properties} \\ \cline{2-6} Populations & E1 & I1 & I2 & external input & delay \\ \hline E1 & 0 & 0 & -15 (-20 – 0) & 6 & 0 (0 – 10) \\ I1 & 15 (0 – 20) & 0 (-20 – 0) & 0 & 6 (0 – 20) & 0 (0 – 10) \\ I2 & 0 & -15 (-20 – 0) & 0 (-20 – 0) & 6 (0 – 20) & 0 (0 – 10) \\ \end{tabular}
\end{table}
Table 2: Parameters of EII network for Fig. 3

\begin{table}
\begin{tabular}{l c c c c c}  & \multicolumn{3}{c}{Synaptic weights} & \multicolumn{2}{c}{population properties} \\ \cline{2-6} Populations & I1 & I2 & I3 & external input & delay \\ \hline I1 & 0 (-20 – 0) & 0 & -15 (-20 – 0) & 6 (0 – 20) & 0 (0 – 10) \\ I2 & -15 (-20 – 0) & 0 (-20 – 0) & 0 & 6 (0 – 20) & 0 (0 – 10) \\ I3 & 0 & -15 (-20 – 0) & 0 (-20 – 0) & 6 (0 – 20) & 0 (0 – 10) \\ \end{tabular}
\end{table}
Table 1: Parameters of III network for Fig. 3 and S1
## Acknowledgments

We thank Kingsshuk Chakravarthy for sharing the code of the basal ganglia network with spiking neurons. We thank Dr. Henri Rhiimaki for helpful comments and suggestions. This work was funded in parts by Swedish Research Council (VR), StratNeuro (to AK), Digital Futures grants (to AK and PH), the Inst. of Advanced Studies, University of Strasbourg, France Fellowship (to AK), and the National Natural Science Foundation of China under Grant No.11572127 and 11872183 (to SL).
