Bifurcation of the neuronal
population dynamics of the modified
theta model: transition tomacroscopic gamma oscillation
Kiyoshi Kotani1
Research Center for Advanced Science and Technology, The University ofTokyo, 4–6–1 Komaba, Meguro-ku, Tokyo 153-8904, Japan
Akihiko Akao2
Department of Mathematics, University of Pittsburgh, Pittsburgh,
Pennsylvania 15213, USA
and Graduate School of Engineering, The University of Tokyo, 7–3–1
Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
Hayato Chiba3
Advanced Institute for Materials Research, Tohoku University,
2-1-1 Katahira, Aoba-ku, Sendai, 980-8577, Japan

Abstract
Interactions of inhibitory neurons produce gamma oscillations (30–80 Hz)
in the local field potential, which is known to be involved in functions suchas cognition and attention. In this study, the modified theta model is considered to investigate the theoretical relationship between the microscopic
structure of inhibitory neurons and their gamma oscillations under a wideclass of distribution functions of tonic currents on individual neurons. The
stability and bifurcation of gamma oscillations for the Vlasov equation of themodel is investigated by the generalized spectral theory. It is shown thatas a connection probability of neurons increases, a pair of generalized eigenvalues crosses the imaginary axis twice, which implies that a stable gamma
oscillation exists only when the connection probability has a value within asuitable range. On the other hand, when the distribution of tonic currents
on individual neurons is the Lorentzian distribution, the Vlasov equation is
reduced to a finite dimensional dynamical system. The bifurcation analyses of the reduced equation exhibit equivalent results with the generalized
1

kotani@neuron.t.u-tokyo.ac.jp
aka64@pitt.edu
3hchiba@tohoku.ac.jp
21

spectral theory. It is also demonstrated that the numerical computations ofneuronal population follow the analyses of the generalized spectral theory as
well as the bifurcation analysis of the reduced equation.

1

Introduction

The local field potential known as collective macroscopic oscillations in thebrain is organized by interacting neurons. Among them, gamma oscillations(30–80 Hz) in the local field potential are known to be involved in functions
such as cognition and attention[1, 2, 3, 4]. In a physiological point of view,
inhibitory synaptic interactions introduced by GABAergic neurons play important role to generate gamma oscillations[1, 5]. However, there are only
limited approaches on how the properties of neurons at the microscopic levelaﬀect macroscopic gamma oscillations.
In this study, the modified theta model (MT model) [6, 7], which is a dynamical system of membrane potentials of interacting neurons, is considered
to investigate the theoretical relationship between the microscopic structureof inhibitory neurons and their gamma oscillations. The Vlasov equationof the model is introdueced under the assumption that the network structure of neurons is the Erdös-Rényi random graph with the probability ofconnection p. For the Vlasov equation, bifurcations from the steady state(de-synchronized state) to the collective gamma oscillations (synchronizedstate) will be studied as the parameter p varies under a wide class of distribution functions of tonic currents on individual neurons.
The stability and bifurcation analysis is not straightforward because of
the continuous spectrum. Let T be a linear operator obtained by the linearization of the Vlasov equation around the steady state. Since the operator T has the continuous spectrum on the imaginary axis, the standard
stability theory of dynamical systems is not applicable. To handle such a
diﬃculty caused by the continuous spectrum, the generalized spectral theory will be employed which was developed to treat a similar problem for theKuramoto model [8, 9]. Although there are no eigenvalues of T when p issmall enough, it will be shown that there exist generalized eigenvalues on
the left half plane. We further assume that the variance of tonic currents on
individual neurons is suﬃciently small. Then, as the parameter p increases,
the generalized eigenvalues cross the continuous spectrum on the imaginary
axis and it implies that the gamma oscillation occurs. Further as p increases,
the eigenvalues again cross the continuous spectrum from the right half planeto the left half plane. As a result, the gamma oscillation is destroyed. These
results imply that if the edge density is too large, the gamma oscillation is
not realized.
2

In Sec.2, the modified theta model and its Vlasov equation are introduced.
A few properties of the steady state will be shown. In Sec.3, the Vlasov equation is reduced to a certain system of evolution equations by Ott-Antonsen
reduction. The linearized system around the steady state with the linear operator T will be derived. In Sec.4, the distribution function for tonic currents
on individual neurons is assumed to be the delta function. In this case, theeigenvalue problem of the operator T is completely resolved. If the distribution function is not the delta function, the operator T has the continuous
spectrum in general, which is studied in Sec.5. In Sec.6, the generalized
spectral theory is applied to investigate the stability and bifurcations of thesteady state. In Sec.7 we consider the Lorentzian distribution of tonic currents where the Vlasov equation is reduced to a finite dimensional dynamicalsystem. It is demonstrated that the numerical computations of neuronal
population follow the analyses of the generalized spectral theory as well as
the bifurcation analysis of the reduced equation.

2

The modified theta model and its steadystate

For the dynamics of the membrane potential V (t) of a Type 1 single neuron,
the following type of quadratic integrate and fire (QIF) model with resting
potential and threshold [11] is sometimes usedcmdV(V − VR )(V − VT )
= gL
+ I,
dt
VT − VR(2.1)

where cm , VR , VT , gL and I are membrane capacitance, resting potential,
threshold potential, leak conductance and input current, respectively. Valuesof these numerical constants used for numerical simulations in this paper are
listed in Table 1, although for mathematical statement, these specific values
are not used. When I = 0, the system has two fixed points V = VR , VT , andmost orbits approach to the stable one VR (if VR < VT ). On the other hand,
if I is larger than some value, two fixed points disappear by a saddle-nodebifurcation and any solution diverges to +∞ in a finite time. We assume
that a neuron fires when V (t) reaches +∞. Then, V (t) is reset to −∞ andagain goes to +∞ periodically.
Next, we consider the large population of neurons with interactionscmdVi(Vi − VR )(Vi − VT )i
= gL
+ Ii − gsyn(t)(Vi − Vsyn ),dt
VT − VRi = 1, · · · , N,(2.2)

iwhere Vsyn is a reversal potential and gsyn(t) is an i-th synaptic conductance3

governed by the equation
i
∑∑
dgsyn
1 i
+ gpeakδ(t − tjk ).
= − gsyndtτconn f ire(2.3)

Here, gpeak and τ denote a peak conductanceand a decay time, respectively∑(see also Table 1). The summation conn runs over all neurons connectedwith the i-th neuron. Thus,of a network of neurons.
∑ it reflects the structurej
The second summation f ire runs over all times tk at which a j’s connected
neuron gets firing. Hence, a neuron aﬀects the dynamics of neighboring
neurons when it fires. Eq.(2.2) is known as the Riccati equation which isextended to a smooth vector field on a compactified phase space S 1 (circle),on which Vi = +∞ is topologically connected with −∞. Specifically, we
introduce the coordinate transformation Vi 7→ θi of the form
Vi =VR + VTVT − VR
θi
+
tan .
2
2
2

(2.4)

Then, Eq.(2.2) is rewritten as the MT model introduced in Kotani et al.
[6, 12]cmdθii= −gL cos θi + c1 (1 + cos θi )Ii + gsyn(t)(c2 (1 + cos θi ) − sin θi ), (2.5)
dt

where c1 = 2/(VT − VR ) and c2 = (2Vsyn − VR − VT )/(VT − VR ) are constants.
By this transformation, V = VT , VR , (VT + VR )/2 and ±∞ are mapped topoints θ = π/2, −π/2, 0 and π on a circle, respectively.
Without loss of generality, we can drop cm by rescaling time (Indeed,
typical physiological value of cm is known to be about 1 µF/cm2 ). For mathematical analysis of the dynamics, it is convenient to consider the continuous
limit (Vlasov equation) of (2.5) given by ∂P + ∂ (vP ) = 0, P = P (t, θ, I)
∂t
∂θ
(2.6)

v(t, θ, I) = −gL cos θ + c1 (1 + cos θ)I + gsyn (t)(c2 (1 + cos θ) − sin θ),
where P (t, θ, I) is a probability density of neurons on S 1 parameterized by
the current I and time t. Further, we assume that a network is the ErdösRényi random graph with the probability of connection p, and consider the
homogenization (averaging with respect to the realization of random graphs)of Eq.(2.3) given by
∫
1
dgsyn
= − gsyn + gpeak · p · N v(t, π, I)P (t, π, I)G(I)dI.
(2.7)
dt
τ
R
4

∑
N denotes the number of neurons, thus the first summationconn is replaced∑by the number of edges pN . The second summation f ire is replaced by theabove integral that denotes the flux of neurons at θ = π at which they fire.
Here, we suppose that the current I is a random variable drawn from a given
distribution function G(I). We introduce the following notation

1
1
1
i

f (I) = − gL + c1 I + c2 gsyn (t) + gsyn (t)



2
2
2
2


 h(I) = c I + c g (t)
1
2 syn
(2.8)
∫
∫



A(t) = v(t, π, I)P (t, π, I)G(I)dI = gL P (t, π, I)G(I)dI



R
R

µ = gpeak · p · N,
where A(t) is a firing rate of the population and the relation v(t, π, I) = gLis used. Then, our system is written as

∂
∂P


+ (vP ) = 0, P = P (t, θ, I)


∂θ

 ∂t
(2.9)
v(t, θ, I) = f (I)eiθ + h(I) + f (I)e−iθ





 dgsyn = − 1 gsyn + µA(t).
dtτIn particular, µ = gpeak · p · N controls the edge density of the network. Our
purpose is to investigate bifurcations of the system from the steady state (desynchronous state) to the synchronous collective firing (gamma oscillation)
as the bifurcation parameter µ varies.
In what follows, we assume that(1) gL , gpeak , c1 are positive and c2 is a negative number (for numerical simulations, we use the specific values shown in Table 1. They are physiologicallyplausible for GABAergic neurons [6]), and(2) G(I) is an even and unimodal function around the average η, which implies G(I + η) = G(−I + η) and G′ (I) < 0 (I > η).
Let us derive the steady state (de-synchronous state) (P, gsyn ) = (P0 , g0 )of the system (2.9). We need the following lemma.
Lemma 2.1. For the quadratic equationf ∈ C, h ∈ R,

f z 2 + hz + f = 0,
the roots
z± =

)
√
1 (
−h ± h2 − 4|f |22f5

Vi (t)membrane potential
i
gsyn
(t)
synaptic conductanceIi
input current
cm
membrane capacitanceVR
resting potential
VT
threshold potential
Vsynreversal potential
gL
leak conductancegpeak
peak conductanceτ
decay timec1
2/(VT − VR )c2
(2Vsyn −VR −VT )/(VT − VR )N
the number of neurons

unknown variableunknown variable
i.i.d. drawn from a distribution G(I)
1 (µF/cm2 )−62 (mV)−55 (mV)−70 (mV)
0.1 (mS/cm2 )0.0214 (mS/cm2 )5 (ms)
2/7
−23/7
800

Table 1: Names of variables and values used for numerical simulations.

satisfy the following properties.
(i) When h2 − 4|f |2 ≤ 0, |z± | = 1.
(ii) When h2 − 4|f |2 > 0, |z+ | < 1 and |z− | > 1.
Let (f0 (I), h0 (I)) denotes (f (I), h(I)) for which the value of the steadystate gsyn (t) = g0 is substituted. Define two sets
Ω1 = {I ∈ R | h0 (I)2 − 4|f0 (I)|2 > 0}
Ω2 = {I ∈ R | h0 (I)2 − 4|f0 (I)|2 ≤ 0}.
Due to Lemma 2.1, the equation f0 (I)eiθ + h0 (I) + f0 (I)e−iθ = 0 has real
solution only when I ∈ Ω2 . Hence, using unknown g0 , the steady state of thefirst equation of (2.9) is given by

C(I)


(I ∈ Ω1 )iθ
f0 (I)e + h0 (I) + f0 (I)e−iθ
P0 (θ, I) =(2.10)δ(θ − θ0 (I))(I ∈ Ω2 ),where C(I) is a normalization constant and δ(θ − θ0 (I)) is the delta function
supported at θ0 (I). θ0 (I) satisfies f0 (I)eiθ0 (I) + h0 (I) + f0 (I)e−iθ0 (I) = 0 and−π < θ0 (I) ≤ 0. The normalization condition provides
∫ 2π
C(I)
1=
dθ.
f0 (I)eiθ + h0 (I) + f0 (I)e−iθ06Changing variables eiθ = z yields
I
I−iC(I)−iC(I)
1=
dz =
dz,
2
f0 (I)(z − z+ )(z − z− )f0 (I)z + h0 (I)z + f0 (I)where the integral path is the unit circle. Lemma 2.1 shows z+ is inside theunit circle while z− is not. Thus, the residue theorem gives
1=1
2πC(I)
f0 (I) z+ − z−⇒

C(I) =1√h0 (I)2 − 4|f0 (I)|2 .
2π

The next task is to calculate g0 . When P is in the steady state P0 (θ, I), we
have
∫
∫
∫
A(t) = gL P0 (π, I)G(I)dI =C(I)G(I)dI + gL
δ(π − θ0 (I))G(I)dI.
R

Ω1

Ω2

Since θ0 (I) ̸= π, the second term in the right hand side disappears. The
steady state gsyn (t) = g0 is given by a solution of the following equation∫ √τµh0 (I)2 − 4|f0 (I)|2 G(I)dI.
(2.11)
g0 = τ µA(t) =2π Ω1Changing the variables from I to x by the relation
h0 (I)2 − 4|f0 (I)|2 = 2gL (c1 I + c2 g0 ) − g02 − gL2 =: x(2.12)

yieldsτµg0 =4πgL c1

∫ ∞
0

√

(
x·G

x + g02 + gL2 − 2c2 gL g0
2gL c1)Further change of variables by x = ω 2 gives∫ ∞τµe (ω) dω,g0 =ω2 · G
4πgL c1 −∞e is an even function defined bywhere G)( 2
ω + g02 + gL2 − 2c2 gL g0e.
G(ω) = G
2gL c1dx.

(2.13)(2.14)(2.15)The steady state (de-synchronous state) (P0 , g0 ) is given by (2.10) and asolution of (2.13) or (2.14), though g0 can be calculated only numerically.

7

Proposition 2.2. Suppose that G(I) is even and unimodal around theaverage η. For any µ ≥ 0, a nonnegative solution g0 = g0 (µ) of (2.13)
uniquely exists. g0 (µ) is monotonically increasing in µ and g0 (µ) → ∞ as
µ → ∞.
Proof. PutτF (z) =
4πgL c1

∫ ∞

√

(
x·G

0

x + z 2 − 2c2 gL z + gL2
2gL c1)dx.

Then, g0 is given by the intersection of the line y = z/µ and the curvey = F (z). Since F (z) > 0 for any z ≥ 0, it is suﬃcient to show that F ′ (z) < 0
for z > 0. Then, F (z) is monotonically decreasing and Proposition is proved.
()∫
τ 2z − 2c2 gL ∞ √x + z 2 − 2c2 gL z + gL2′′
F (z) =x·G
dx.
4πgL c1 2c1 gL
2gL c10
Because of our assumption for the signs of parameters, the number in thefront of the integral is positive (recall that c2 < 0).
By the assumption for G, G′ (I) is an odd function around I = η satisfying
G′ (I) < 0 for I > η and G′ (η) = 0. Thus, there is a unique number x0 suchthat)(x + z 2 − 2c2 gL z + gL2′
H(x) = G
2gL c1is an odd function around x = x0 satisfying H(x) < 0 for x > x0 andH(x0 ) = 0. If x0 ≤ 0, then H(x) < 0 for x > 0 and F ′ (z) is negative. If
x0 > 0, we decompose the above integral as
∫ ∞
∫ x0
∫ 2x0
∫ ∞
√
√
√
√
xH(x)dx =xH(x)dx +xH(x)dx +xH(x)dx.
0

0

x0

2x0

By the assumption, the sum of the first two terms is negative, and the third
term is also negative, which proves that F ′ (z) < 0.

3

Ott-Antonsen reduction

For the system (2.9), we expand a solution in a Fourier series. Putting∫ 2π
P (t, θ, I)eikθ dθ,zk (t, I) =0

8

we rewrite (2.9) as the system of equations of zk (t, I)’s given byż1 = i( f (I)z2 + h(I)z1 + f (I) )






żk = ik( f (I)zk+1 + h(I)zk + f (I)zk−1 )



1
(3.1)
ġsyn = − gsyn + µA(t)

τ

)
(

∫

∞
∞

∑
∑

gL


z
)
−
(z2k−1 + z 2k−1 ) G(I)dI.
A(t)
=
1
+
(z
+
2k
2k

2π R
k=1
k=1
It is easy to verify that the set in the phase space defined by {zk = z1k | k =1, 2 · · · } is an invariant set of this system, on which it is reduced to thesystem (Ott-Antonsen reduction [13, 14, 15]). Thus, the reduced equations,which has been introduced in Ref.[12], are given byż1 = i( f (I)z12 + h(I)z1 + f (I) )





2


 ż 1 = −i( f (I) z 1 + h(I)z 1 + f (I) )1(3.2)
ġsyn = − gsyn + µA(t)



τ
)
∫ (



z1
gL
z1

 A(t) =1−−
G(I)dI.
2π R
1 + z1 1 + z 1
We regard this system as a dynamical system on the Hilbert space H :=L2 (R, G(I)dI) ⊕ L2 (R, G(I)dI) ⊕ R (the equation for the complex conjugate
z 1 is added because L2 (R, G(I)dI) is considered as an R-vector space).
Lemma 3.1. The steady state (P0 , g0 ) shown in Section 2 lies on theinvariant set.
Proof. We show that the steady state satisfies zk = z1k . When I ∈ Ω2 ,∫ 2π
δ(θ − θ0 (I))eikθ dθ = eikθ0 (I) = z1k .
zk =0When I ∈ Ω1 ,

∫ 2π
zk =
0

C(I)
f0 (I)eiθ + h0 (I) + f0 (I)e−iθ

eikθ dθ.

As in the previous section, this integral is calculated by using the residue
theorem with the definition of C(I) aszk =k
2πC(I) z+k= z+= z1k .
f0 (I) z+ − z−9This proves the lemma.
In what follows, we investigate the linear stability of the steady state(z1 , z 1 , g0 ) = (z+ (I), z+ (I), g0 ) of Eq.(3.2) given by
z+ (I) =)√1 (−h0 (I) + h0 (I)2 − 4|f0 (I)|2
2f0 (I)and a positive solution of (2.13). The linearized equation of (3.2) around thissteady state is
 
  
 
Q(I)
0
P (I)
u
u
u
d   u =0Q(I) P (I)   u  := T  u  ,
(3.3)
dt
v
v
v
R(I) R(I) −1/τ
where
√

Q(I)
=
i
h0 (I)2 − 4|f0 (I)|2 ,




i
P (I) = (1 + z+ )((c2 + i)z+ + c2 − i),
∫
2

u
µg


 R(I)u = − L
G(I)dI.
2π R (1 + z+ )2
Q(I) and P (I) are multiplication operators on L2 (R, G(I)dI) and R, respectively, and R(I) is the integral operator on L2 (R, G(I)dI). Thus, the above
matrix defines a linear operator T on H. For the linear stability of the steadystate, we will investigate the spectrum of T . The following relations will beoften used.
h0 (I)2 − 4|f0 (I)|2 = 2gL (c1 I + c2 g0 ) − g02 − gL2 ,
)
1 (√
P (I)
2 − 4|f (I)|2 + i(c g − g ) .
=
h
(I)
0
0
2 L
0
(1 + z+ )2
2gL

4

(3.4)
(3.5)

Delta distribution

First, we consider the simplest case: we suppose that the distribution G(I) =δ(I − η) is the delta function supported at η. In this case, the integral
operator is not involved and the system (3.2) is reduced to a 3-dim system.
The equation (2.11) to obtain g0 is reduced to
√
µτ
µτ √
2
2
h(η) − 4|f (η)| =2gL (c1 η + c2 g0 ) − g02 − gL2 .
g0 =
2π
2π

10

The eigen-equation of the linearized system around the steady state (z+ (η), z+ (η), g0 )is given by
λ − Q(η)
0−P (η)
0λ − Q(η)−P (η)  = 0. (4.1)
f (λ) := det 

 µgL
1
µgL
1
λ + 1/τ2
2π (1 + z+ (η))
2π (1 + z+ (η))2When µ√= 0, g0 = 0. Thus, three eigenvalues are given by λ = −1/τ and
λ = ±i 2gL c1 η − gL2 . Since we are interested in the macroscopic oscillation
which will appear through a Hopf bifurcation, we assume that 2c1 η − gL > 0;i.e. there exist two eigenvalues on the imaginary axis when µ = 0.
Proposition 4.1. Suppose η > gL /(2c1 ) and 1/τ < −c2 gL . For anyµ > 0, f (λ) = 0 has a negative real root α and a pair of complex rootsβ, γ = β whose real parts are positive. The pair transversely crosses theimaginary axis at µ = 0 from the left half plane to the right half plane.
Proof. As in Eq.(2.12), we put(

x = 2gL (c1 η + c2 g0 ) − g02 − gL2 =

2πg0
τµ

)2
.

With the aid of Eqs.(3.4) and (3.5), we can show that the eigen-equation(4.1) is rearranged as
1
1
1
f (λ) = λ3 + λ2 + (x + g0 )λ + (x + g02 − c2 gL g0 ) = 0τττThe relation between roots and coeﬃcients implies111
α + β + γ = − , αβ + βγ + γα = x + g0 , αβγ = − (x + g02 − c2 gL g0 ).
τττSince all coeﬃcients are positive, it turns out that there are no positive realroots. By the assumption, we have
g0
1
f (−1/τ ) = (g0 − − c2 gL ) > 0,ττwhich implies that there exists a negative real root α satisfying α < −1/τ .
This shows β + γ = −1/τ − α > 0. Hence, if β and γ are real, at least oneof them is positive. This is a contradiction and they should be a complex
pair with positive real parts. The last statement follows from the implicit
function theorem to f (λ) = 0.
Fig.1 represents the motions of eigenvalues as µ increases from zero. The
proposition implies that a Hopf bifurcation occurs at µ = 0 and the gamma
oscillation is stable for any µ > 0.
11

Figure 1: (Color online)The motions of eigenvalues as µ increases from zero.

5

Stability analysis 1

In general, the system (3.2) is an infinite dimensional dynamical system on
H and the linear operator T has the continuous spectrum. Let us derive theeigen-equation of the operator T defined by (see Eq.(3.3)) 
λ − Q(I)−P (I)0u0λ − Q(I) −P (I)   u  = 0.
v
−R(I)−R(I) λ + 1/τ
The first two equations give u = (λ − Q(I))−1 P (I)v. By substituting it into
the third equation, we obtain
)
∫ (
µgL
1
P (I)
1
P (I)
1
λ+ =−
+
G(I)dI.(5.1)τ
2π R (1 + z+ )2 λ − Q(I) (1 + z+ )2 λ − Q(I)The factors (λ − Q(I))−1 and (λ − Q(I))−1 suggest that if the support of G is
the whole real axis, which will be assumed in Sec.6, the continuous spectrumof T is given by the set{Q(I) | I ∈ R} ∪ {Q(I) | I ∈ R} = iR ∪ R<0 ,(5.2)

that is, the whole imaginary axis and the negative real axis (Fig.2). Becauseof the spectrum on the imaginary axis, the steady state looks neutrally stablein H-topology. In the next section, we will employ the generalized spectral
theory to treat the continuous spectrum.

12

(a)

(b)

Figure 2: (Color online)(a) The continuous spectrum of the operator T .
(b) The generalized spectrum of T . The dotted lines denote the motion ofgeneralized eigenvalues and solid lines denote that of usual eigenvalues as µincreases. See Section 6.

We further simplify Eq.(5.1). When I ∈ Ω2 , Q(I) = Q(I) is a real
number. We also have
P (I)
P (I)
+
2
(1 + z+ )(1 + z+ )2()
i (c2 + i)z+ + c2 − i (c2 − i)z+ + c2 + i
1 − |z+ |2
.
=
−
=
2
1 + z+
1 + z+|1 + z+ |2
Due to Lemma 2.1, |z+ | = 1 and the above quantity becomes zero. Thus, we
can assume I ∈ Ω1 in Eq.(5.1):
)
∫ (
1
µgL
P (I)
1
1
P (I)
λ+ =−
+
G(I)dI.
τ2π Ω1 (1 + z+ )2 λ − Q(I) (1 + z+ )2 λ − Q(I)
Changing the variables by (2.12) and using the formula (3.5) yield
1
=
τ ∫
(√
)
)
(
√
∞
1
µ
x + i(c2 gL − g0 )x − i(c2 gL − g0 )x + g02 + gL2 − 2c2 gL g0
√
√
dx.
−
·G
+
4πc1 0 2gL
2gL c1
λ−i x
λ+i x
λ+

Further, we put x = ω 2 to obtain
)
∫ ∞(
1
µ
ω + i(c2 gL − g0 ) ω − i(c2 gL − g0 ) e
λ+ =−+G(ω)ωdω,τ
4πgL c1 0
λ − iω
λ + iω
13eewhere G(ω)
is defined by Eq.(2.15). Note that G(ω)
is an even function.
Thus, putting ω 7→ −ω only for the second term gives∫ ∞1µω + i(c2 gL − g0 ) eω G(ω)dω,
λ+ =−
(5.3)
τ
4πgL c1 −∞λ − iω
Finally, by using the relation
ω2iωλ= iω −,
λ − iω
λ − iω
we obtain1µλ+ =(λ + g0 − c2 gL )
τ
4πgL c1∫ ∞iω eG(ω)dω.
−∞ λ − iω(5.4)

When µ = 0, there are no interactions of neurons and de-synchronous stateshould be asymptotically stable. Thus, one may expect the existence of
eigenvalues on the left half plane. Unfortunately, this is not true.
Lemma 5.1. When µ = 0, there are no eigenvalues with nonzero realparts.
Proof. The right hand side of Eq.(5.4) is a bounded holomorphic functionin λ on the right and left half planes. Thus, a possible eigenvalue for µ = 0 is
only λ = −1/τ . However, it is not a true eigenvalue because it is embedded in
the continuous spectrum and the corresponding eigenvector does not belong
to H.
Lemma 5.2. Suppose µ > 0.
(i) If λ is an eigenvalue, so is its complex conjugate.
(ii) There are no eigenvalues on the real axis.
This lemma implies that if a bifurcation occurs at some µ > 0, it will bea Hopf bifurcation.
Proof. Putting λ = x + iy in (5.3), we obtain
∫ ∞

1
µ
xω + (c2 gL − g0 )(y − ω) e
ω G(ω)dω x + τ = − 4πg cx2 + (y − ω)2
L 1 −∞∫ ∞ 2ω − yω + (c2 gL − g0 )x e
µ

 y=−
ω G(ω)dω.
4πgL c1 −∞x2 + (y − ω)2e is an even function, it is easy to verify that if (x, y) is a solution, so
Since Gis (x, −y). Further, y = 0 always satisfies the second equation. Then, the14

first equation gives
1
µ
x+ =−
τ
4πgL c1∫ ∞x − c2 gL + g0 2 e
ω G(ω)dω
x2 + ω 2−∞

If x > 0, the left hand side is positive, while the right hand side is negative. If
x ≤ 0, it is embedded in the continuous spectrum and not a true eigenvalue.

6

Stability analysis 2 (generalized spectral theory)

In Section 5, the spectrum of the operator T on H = L2 (R, G(I)dI) ⊕L2 (R, G(I)dI) ⊕ R is conisdered: The operator T has the continuous spectrum on the imaginary axis and there are no eigenvalues on the left half plane
when µ = 0. This implies that the de-synchronous state is neutrally stable
but not asymptotically stable in the topology of H. However, the unknown
function P (t, θ, I) of Eq.(2.9) is not a L2 function but a distribution. Thus,
it is better to study the dynamics in a weaker topology rather than that of
H. For this purpose, we employ the generalized spectral theory in Section 6.
The generalized spectral theory was developed by Chiba [8, 9] to treat problems related to continuous spectra. He defined generalized eigenvalues andproved that it plays a similar role to usual eigenvalues. In this section, a brief
review of the generalized spectral theory is given. All proofs are included in[9]. Then, we will show that the operator T has the generalized eigenvalues
and they induce the stability and bifurcations of the de-synchronous state.
For a given function f (ω) on R, let us consider the function∫
1
A(λ) =f (ω)dωR λ − iωof λ ∈ C. It is holomorphic on the right half plane, though singular on theimaginary axis in general because of the factor (λ − iω)−1 . However, if f (ω)
is a “nice” function, A(λ) may be well-defined on the imaginary axis. The
following lemma known as Sokhotskii formula [16] is fundamental.
Lemma 6.1. If f (ω) is a holomorphic function around the imaginary
axis, A(λ) is well-defined on the imaginary axis and it has an analytic con-15

tinuation from the right half plane to the left half plane given by ∫1
f (ω)dω(Re(λ) > 0)
R λ − iω


∫

1
lim
f (ω)dω(Re(λ) = 0)A(λ) =
Re(λ)→0+ R λ − iω



∫


1


f (ω)dω + 2πf (−iλ) (Re(λ) < 0).
R λ − iω(6.1)

In [9], this idea is applied not to functions but to linear operators. Let T be a
linear operator defined on a Hilbert space H having the continuous spectrum
on the imaginary axis. By the definition of the spectrum set, the resolvent
operator (λ−T )−1 is singular on the imaginary axis. To consider the analytic
continuation of (λ−T )−1 , let X ⊂ H be a “nice” dense subspace of H and X ′be its dual space (the set of continuous linear functionals on X). If a topology
of X is stronger than that of H, the dual space is larger than H and we havea triplet X ⊂ H ⊂ X ′ , which is called the Gelfand triplet. We can provefor some class of operators that if the resolvent (λ − T )−1 is regarded as an
operator from X into X ′ , it has an analytic continuation from the right halfplane to the left half plane beyond the continuous spectrum on the imaginary
axis. If the analytic continuation has a singularity on the left half plane, it
is called a generalized eigenvalue. Namely, if the domain of the resolvent
is restricted to a “nice” space X (which depends on the problem at hand.
See [9, 10] for several examples), then the continuous spectrum disappears.
Instead of the continuous spectrum, generalized eigenvalues appear on theleft half plane. An associated eigenvector of a generalized eigenvalue is notincluded in H but an element of the dual space X ′ . Thus, we cannot find ageneralized eigenvalue in the Hilbert space theory. A generalized eigenvalue
plays a similar role to a usual eigenvalue. If it exists on the left half plane,
it induces an exponential decay of a solution of a linear system du/dt = T u
with respect to the topology of X ′ , which is weaker than that of H.
For our operator T given in (3.3), the resolvent has an analytic continuation to the left half plane if the domain is restricted to a certain space ofholomorphic functions. The formula to give generalized eigenvalues is easily
obtained from (5.4) by applying Lemma 6.1. For it, we assume that the function G(I) defined on R has a continuation to the complex plane and regard
it as a function on C.
Proposition 6.2. Suppose that G(I) has a meromorphic continuation
to the complex plane. The analytic continuation of the eigen-equation (5.4)

16from the right half plane to the left half plane is given by
λ+

1
=
τ

∫
µ
iω e
(λ + g0 − c2 gL )G(ω)dω
(Re(λ) > 0)4πgL c1
λ − iω
)
(∫
(6.2)

µ
iω e

e

G(ω)dω + 2πλG(−iλ)(Re(λ) < 0)
(λ + g0 − c2 gL )4πgL c1
λ − iω

A root of the first line is a usual eigenvalue, while that of the second line
gives a generalized eigenvalue of T . It is easy to confirm that if λ is a root,
then so is the complex conjugate λ.
Theorem 6.3. Suppose that G(I) has a meromorphic continuation to
the complex plane without poles on the real axis.
(i) When µ = 0, T has generalized eigenvalues λ = λ(µ) on the left halfplane. One of them is λ = −1/τ , and the others are determined by poles ofethe function G(ω)
on the upper half plane.
(ii) Suppose that η > gL /(2c1 ) and 1/τ < −c2 gL . If the variance of G(I)is suﬃciently small, generalized eigenvalues λ(µ) cross the imaginary axis
transversely from left to right, and become usual eigenvalues.
(iii) If µ is suﬃciently large, there are no (generalized) eigenvalues on theregion {λ | Re(λ) ≥ −1/τ }. In particular, there are no eigenvalues on theright half plane.
See Fig.2 (b) for a schematic picture of the motions of generalized eigenvalues. This result implies that a pair of (generalized) eigenvalues crosses
the imaginary axis at least twice, so that a Hopf bifurcation occurs twice(see also Lemma 5.2). At the first bifurcation, the stable synchronous firing
state (gamma oscillation) bifurcates from the de-synchronous state, and at
the second bifurcation, it disappears and the de-synchronous state becomesstable again. It is remarkable that if the edge density µ of the network is toolarge, synchronization does not occur.
Proof. (i) λ = −1/τ is a trivial solution of (6.2) when µ = 0. To findeanother solution for µ = 0, note that if G(ω)
has a pole ω∗ on the upper halfeplane, then iω∗ is a pole of G(−iλ)
on the left half plane. Let us express theLaurent series ase
G(−iλ)
=

dp
dp−1
+
+ ···
p
(λ − iω∗ )(λ − iω∗ )p−1

for some positive integer p. Substituting it into the second line of (6.2) and17

multiplying the factor (λ − iω∗ )p giveµ
(λ + g0 − c2 gL ) · 2πλdp + O(λ − iω∗ ) = 0,
4πgL c1
as λ → iω∗ . Here note that the integral in (6.2) is bounded at λ = iω∗because of Lemma 6.1. This proves that iω∗ is a solution for µ = 0.
(ii) Since G(I − η) is an even and unimodal function around the average
η, it converges to the delta function supported at η as the variance tends tozero with respect to the topology of Schwartz distributions. Thus, Eq.(5.4)
converges to the eigen-equation for the delta function (4.1). Fix two values
0 < µ1 << 1 and 1 << µ2 < ∞. On the closed interval [µ1 , µ2 ], the trajectoryof a (generalized) eigenvalue λ(µ) is uniformly approximated by that of theeigenvalue for the delta function case if the variance is suﬃciently small.
Then, the statement follows from Proposition 4.1.
(iii) Eq.(5.4) is also written as
)
∫ ∞(
−µ
λ
1
e
G(ω)dω.
(λ + g0 − c2 gL )
1−
λ+ =
τ
4πgL c1
λ − iω
−∞By Proposition 2.2, g0 (µ) → ∞ as µ → ∞. By changing the variable asω = g0 ω̂, the above integral is rewritten as) ( 2 2)∫ ∞(λg0 (ω̂ + 1) + gL2 − 2c2 gL g0
g0
1−
G
dω̂.
λ − ig0 ω̂
2gL c1−∞∫∞eThis shows that the second term −∞ λ/(λ − iω)G(ω)dωis infinitesimally∫∞esmaller than the first term −∞ G(ω)dω
as g0 → ∞. Hence, considering only
leading terms as µ → ∞, we obtain
∫ ∞
1
−µ
e
G(ω)dω
+ o(µ).
λ+ =(λ + g0 )τ
4πgL c1−∞
If Re(λ) > −1/τ , the real part of the left hand side is positive, while that of
the right hand side is negative. This proves that there are no eigenvalues on
the region Re(λ) > −1/τ . The same argument is also valid for the secondline of (6.2). Thus, there are no generalized eigenvalues on the same region.
eRemark. If G(ω)
does not have poles on the upper half plane, there areno generalized eigenvalues except −1/τ when µ = 0. Then, the statement(ii) implies that generalized eigenvalues pop out from the infinity ∞ to theleft half plane when µ = +0, and they approach to the imaginary axis as thevariance tends to zero for fixed µ > 0.
Although the Lorentzian distribution G(I) = (∆/π)/((I − η)2 + ∆2 ) does
not have the variance, the statement (ii) is applicable because it converges
to the delta function as ∆ → 0.
18

7

Numerical results

Here we consider the Lorentzian distribution G(I) = (∆/π)/((I − η)2 + ∆2 ).
Then, the reduced equation (Eq. (3.2)) can be further simplified by virtue ofthe residue theorem. We also introduce the order parameter α(t) ∈ C as
∫ ∫ 2πα(t) =
P (t, θ, I)eiθ G(I)dθdI.
(7.1)
R

0

Ott Antonsen ansatz and the residue theorem yield the following relations
α(t) = z1 (t, η + i∆),
[
(
)]
gL
−α(t)
A(t) =
1 + 2Re.
2π
1 + α(t)

(7.2)
(7.3)

By using these relations, the dynamics of the order parameter and the synaptic conductance are described asα̇ = i( f (η + i∆)α2 + h(η + i∆)α + f (η + i∆) )





1
ġsyn = − gsyn + µA(t)
(7.4)
τ[
(
)]



g
−α(t)

 A(t) = L 1 + 2Re,
2π
1 + α(t)which has been presented in Ref. [12]. We employ bifurcation analyses forthe equations. Here, we adopt (η, ∆) = (2, 0.05) and set µ as a bifurcation
parameter. The bifurcation diagram is shown in Fig. 3(a). We can see thatas µ increases, the populational gamma oscillation emerges at µ ∼ 0.18 via
Hopf bifurcation. The oscillatory state disappears by further increase of µ
at µ ∼ 4.7.
These dynamics are confirmed by direct numerical computations of the individual neuron (Eq. (2.5)) with homogenized synaptic dynamics (Eq. (2.7))
to avoid the eﬀects of random variation and finite size eﬀect that appear theindividual network structure. We show raster plots of individual firings with
µ = 8.6 · 10−2 , 3.2, and 15.0 in Figs. 3(b-1),(c-1), and (d-1), respectively. The
vertical axis is a neuron index that is sorted in ascending order of I. We cansee that the emergence and disappearance of the populational oscillation coincides with the bifurcation analysis of the macroscopic equation (Eq. (7.4)).
Additionally, we show time-course of g with µ = 8.6 · 10−2 , 3.2, and 15.0 in
Figs. 3(b-2),(c-2), and (d-2), respectively. The ferquency of the oscillation in
Figs. 3(c-2) is 34 Hz, therefore, it is gamma frequency oscillation. Numericalresults for the individual neuron model (Eq. (2.5) and (2.7)) well match tothose for the macroscopic equation (Eq. (7.4)).
19

Figure 3: (Color online)(a) Bifurcation diagram of Eq. (7.4). The green
line is a stable steady state and the red one is an unstable one. The bluelines indicate the upper and lower bounds of the gamma oscillation. (b-1)
Raster plot with µ = 8.6 · 10−2 , (b-2) Time-course of the conductance g with
µ = 8.6 · 10−2 . The blue line is a result of the individual neuron model(Eq. (2.5) and (2.7)) and the red line is that of the macroscopic equation(Eq. (7.4)). (c-1) Raster plot with µ = 3.2, (c-2) Time-course of g obtainedby the neuron model (blue) and the macroscopic equation (red) with µ = 3.2.
(d-1) Raster plot with µ = 15.0, (c-2) Time-course of g obtained by theneuron model (blue) and the macroscopic equation (red) with µ = 15.0.

208Discussion

For the analysis of populational neuronal firings with steady-states or slow
frequency oscillations, it is eﬀective to use the firing rate model, such as
Wilson-Cowan model [17], derived by the adiabatic approximation. On theother hand, for analysis of oscillations with high frequency such as gamma
oscillations, the firing rate model is not appropriate as an approximation [18].
Instead, Vlasov equation of the MT model posesses the dynamics of distribution itself. This is why we utilized the MT model for the analyses of gamma
oscillations. Although some studies have analyzed the coupled dynamics ofthe conventional theta neuron model [14, 15, 19, 20], the strength of using
the MT model lies in its ability to evaluate the collective dynamics undersynaptic interactions with physiologically appropriate conductance and reversal potential.
The inhibitory synaptic interactions by GABAergic neurons are known
to generate a macroscopic oscillation in the gamma frequency range, which isnamed ING (interneuron gamma oscillation) [1, 5]. In addition, synaptic connectivity in the cerebral cortex is diﬀerent from individual to individual anddepends on the developmental stage. It has been reported that the synaptic
connections increase in early phase of development followed by decrease by
activity-dependent synaptic pruning [21, 22, 23].
Therefore, we investigate the relationship between the synaptic connectivitiy and macroscopic gamma oscillations in the population of the inhibitory
neurons under meromorphic function G(I) that is even and unimodal aroundthe average and has no poles on the real axis. If the variance of G(I) is suﬃciently small, the generalized spectral theory indicate that the populational
gamma oscillation emerges only when the synaptic connectivity has a value
within an approariate range, which is confirmed by the numerical bifurcation
analysis under the Lorentzian distribution for G(I). As shown in Fig. 3(c-2),
the macroscopic oscillation with frequency of 33.6 Hz is observed, which iswithin the range of gamma oscillation (30–80 Hz). In addition, the dynamics of two asynchronous states by the numerical computations (Fig. 3(b-1)
and (d-1)) look very diﬀerent. In the asynchronous state with smaller µ,almost every neuron fires aperiodically (Fig. 3(b-1)). On the contrary, inthe asynchronous state with larger µ, only partial population fires and others keep silent (Fig. 3(b-1)). We could say that normal neuronal populations
avoid the latter state that has excess connections between inhibitory neurons.
Besides biological plausibility, the second bifurcation seems mathematicallymeaningful. The generalized spectral theory indicates that these transitionsof gamma oscillations are introduced by the nontrivial motion of the pair of
the generalized eigenvalues that cross the imaginary axis twice. Remarkably,both types of asynchronous state are well characterized by the generalized21

eigenvalues on the left half plane.

Acknowledgment
This study is supported in part by JST PRESTO (JPMJPR14E2) to KK,
and by JST PRESTO (JPMJPR16E7) to HC.

References[1] X. J. Wang. Neurophysiological and computational principles of corticalrhythms in cognition. Physiological reviews, 90(3), 1195-1268 (2010).
[2] T. Womelsdorf, J. M. Schoﬀelen, R. Oostenveld, W. Singer, R. Desimone, A. K. Engel, and P. Fries. Modulation of neuronal interactions
through neuronal synchronization. Science 316, 1609-1612, (2007).
[3] E. Rodriguez, N. George, J. P. Lachaux, J. Martinerie, , B. Renault,
and F. J. Varela. Perception’s shadow: long-distance synchronization ofhuman brain activity. Nature, 397(6718), 430-433 (1999).
[4] C. M. Gray, P. König, A. K. Engel, and W. Singer. Oscillatory responses
in cat visual cortex exhibit inter-columnar synchronization which reflects
global stimulus properties. Nature, 338(6213), 334 (1989).
[5] M. Bartos, I. Vida, and P. Jonas, Synaptic mechanisms of synchronizedgamma oscillations in inhibitory interneuron networks, Nat. Rev. Neurosci., vol. 8, no. 1, pp. 45 窶・6, (2007)[6] K. Kotani, I. Yamaguchi, L. Yoshida, Y. Jimbo, and B. Ermentrout,
Population dynamics of the modified theta model: Macroscopic phasereduction and bifurcation analysis link microscopic neuronal interactions
to macroscopic gamma oscillation, J R Soc Interface. 11(95): 20140058,(2014).
[7] A. Akao, Y. Ogawa, Y. Jimbo, B. Ermentrout, and K. Kotani. Relationship between the mechanisms of gamma rhythm generation and themagnitude of the macroscopic phase response function in a population of
excitatory and inhibitory modified quadratic integrate-and-fire neurons.
Physical Review E, 97(1), 012209. (2018)[8] H. Chiba, A proof of the Kuramoto conjecture for a bifurcation structureof the infinite-dimensional Kuramoto model, Ergodic Theory Dynam.
Systems 35, no. 3, 762-834, (2015).
[9] H. Chiba, A spectral theory of linear operators on rigged Hilbert spaces
under analyticity conditions, Adv. in Math. 273, 324-379, (2015).
22

[10] H. Chiba, A spectral theory of linear operators on rigged Hilbert spaces
under analyticity conditions II : applications to Schrodinger operators,
Kyushu Journal of Math. 72, 375-405 (2018).
[11] W. Gerstner and W. Kistler . Spiking Neuron Models. Single Neurons,
Populations, Plasticity. Cambridge University Press (2002)[12] A. Akao, Y. Ogawa, B. Ermentrout, Y. Jimbo, and K. Kotani. Macroscopic dynamics of heterogeneous theta neurons connected with chemical
synapses, SIAM conference on Life sciences 2016, PP1, (2016)[13] E. Ott and T. M. Antonsen, Low dimensional behavior of large systemsof globally coupled oscillators, Chaos 18 (2008), no. 3, 037113.
[14] T. B. Luke, E. Barreto, and P. So, Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons, NeuralComputation. 25(12), 3207-3234. 2013[15] P. So, T. B. Luke, and E. Barreto, Networks of theta neurons with timevarying excitability: Macroscopic chaos, multistability, and final-state
uncertainty, Phys. D Nonlinear Phenom., 267, 16 窶・6, 2014[16] F. D. Gakhov, Boundary value problems, Dover Publications, Inc., New
York, 1990, Translated from the Russian, Reprint of the 1966 translation.
[17] H. R. Wilson and J. D. Cowan. Excitatory and inhibitory interactions
in localized populations of model neurons. Biophys. J., 12:124, (1972)[18] B. Ermentrout and D. H. Terman, Mathematical foundations of neuroscience (Vol. 64). Springer (2010).
[19] B. Gutkin, T. Hely,and J. Jost. Noise delays onset of sustained firing ina minimal model of persistent activity. Neurocomputing, 58, 753-760.
(2004)[20] Z. P. Kilpatrick and B. Ermentrout. Sparse gamma rhythms arisingthrough clustering in adapting neuronal networks. PLoS computational
biology, 7(11), (2011).
[21] P. R. Huttenlocher. Morphometric study of human cerebral cortex development. Neuropsychologia, 28(6), 517-527, (1990)[22] N. Uesaka, M. Abe, K. Konno, M. Yamazaki, K. Sakoori, T. Watanabe,
T-H. Kao, T. Mikuni, M. Watanabe, K. Sakimura, and M. Kano. Retrograde signaling from progranulin to Sort1 counteracts synapse elimination in the developing cerebellum. Neuron, 97(4), 796-805, (2018).
[23] X. Wu, Y. Fu, G. Knott, J. Lu, G. Di Cristo, and Z. J. Huang. GABA
signaling promotes synapse elimination and axon pruning in developing
cortical inhibitory interneurons. J. Neurosci. 32(1), 331-343, (2012)23