Multistable attractors in a network of phase oscillators with three-body interaction
Takuma Tanaka∗
Department of Computational Intelligence and Systems Science, Tokyo Institute of Technology, Japan.

Toshio Aoyagi
Graduate School of Informatics, Kyoto University, Japan andJST, CREST(Dated: May 2, 2011)
Three-body interactions have been found in physics, biology, and sociology. To investigate theireﬀect on dynamical systems, as a ﬁrst step, we study numerically and theoretically a system of phaseoscillators with three-body interaction. As a result, an inﬁnite number of multistable synchronizedstates appear above a critical coupling strength, while a stable incoherent state always exists forany coupling strength. Owing to the inﬁnite multistability, the degree of synchrony in asymptotic
state can vary continuously within some range depending on the initial phase pattern.

To analyze the neuronal networks with three-body interaction, we exploit the fact that neurons exhibit periodic ﬁrings in many cases. Periodic activities are ubiquitous in not only neuronal networks, but also phenomena studied in other ﬁelds of biology, including gene expression in E. coli, synchronous ﬂashing of ﬁreﬂies, and(a) 500
400
Neuron

(c)

(b)

300
200
100
0
29985 29990 29995 30000
t

80
78
76
74
72
70
0 100 200 300 400 500
Neuron

1
0.8
0.6
R

Interaction among particles or elements in classicalmechanics, electromagnetism, and many other ﬁelds ofphysics is often modeled by two-body interaction. Description by the linear superposition of two-body interactions has allowed us to predict the future orbits of theplanets and to design drug molecules that tightly bindto the target protein. However, it has been revealed thatthe net interaction experienced by an element cannot be
written as the linear superposition of the two-body interaction in several systems, including physical systems [1],social and economic systems [2], and neuronal networks[3]. A typical example is signal transmission from one
neuron to another. The signals are mediated by the release of neurotransmitters from synapses, and some neurotransmitters modulate the response of neurons to inputs from other neurons (heterosynaptic plasticity) [4].
This modulation can be regarded as three-body interaction, although synaptic transmission is conventionallymodeled as a two-body interaction. To show what occurs in such neuronal networks with three-body interactions, we present a numerical simulation of a networkof Hodgkin-Huxley neurons [5] with short-term heterosynaptic plasticity (see [6] for details). In this model, theinput from neuron j to i is modulated by the relative
spike timing of neuron i and other neurons in this model.
Figure 1(a) shows that this neuronal network exhibits
multistability, in which the numbers of synchronized neurons at the steady state vary depending on the initialconditions [Fig. 1(b)]. This seems to be a novel behavior
not observed in systems with only two-body interactions.
However, this system is too complicated to show analytically why this multistability arises.

Firing rate (Hz)PACS numbers: 05.45.Xt,05.65.+b

0.4
0.2
0

0.1

1

10
t

100

1000

FIG. 1. (color online). Two examples of multistability arisingfrom three-body interaction. (a) Raster plot of the spikes in a
network of the Hodgkin-Huxley neurons with short-term heterosynaptic plasticity (N = 500). Black, red, and blue dots
represent the ﬁrings in the networks starting from diﬀerent
initial conditions. Neurons are sorted in ascending order ofthe ﬁring rate. Firing rate of each neuron is shown in (b).
(c) Time evolution of the order parameter R for three different initial conditions in the phase-oscillator systems with
non-uniform random coupling (N = 500). Phase distributions of oscillators at t = 0 and t = 1, 000 are shown on theleft and right, respectively. These two examples demonstratethat the same system can show diﬀerent degrees of synchrony
depending on the initial conditions.

pedestrians’ gait [7]. The behavior of these periodic activities is described by a form of phase oscillators in aquite general context [8]. However, three-body interaction among phase oscillators has not been studied yet.
Since phase oscillators are simple enough to be analytically tractable and structurally stable, theory of phaseoscillators is a powerful tool in interpreting and elucidating complicated experimental results in which three-body2interactions play an essential role. In this Letter, we thus
examine the eﬀect of three-body interaction on the dynamics of globally-coupled phase oscillators.
As a natural extension of the system of limit-cycle
oscillators with two-body interaction, N -oscillator system with two- and three-bodyinteraction is described∑by Ẋi = Fi (Xi ) + j,k Vijk (Xi , Xj , Xk ), where Fi
describes the dynamics of uncoupled oscillator i and
Vijk is the phase coupling function. Two-body interaction Vij (Xi , Xj ) is then included as a special case of
the three-body interaction Vijk (Xi , Xj , Xk ). Using thephase reduction technique, we can describe the dynamics
of oscillator i with one variable, phase φi . Thus, the dynamics of the system of phase oscillators with three-bodyinteraction is generally given by∑φ̇i = ωi +
Γijk (φji , φki ),(1)j,kwhere ωi is the natural frequency of oscillator i, φji =φj − φi , and Γijk is the coupling function.
We present one example in which typical novel features
arising from three-body interactions can be seen:
1 ∑[aij sin(φji + α1ij ) + bij sin(2φji + α2ij )]φ̇i = ωi +
N j
1 ∑
+ 2
cijk sin(φji + α3ijk ) cos(φki + α4ijk ),Nj,kwhere aij , bij ∼ N (0.3, 0.01), cijk ∼ N (6, 4), andα1ij , α2ij , α3ijk , α4ijk ∼ N (0, 0.09). Here N (µ, σ 2 ) denotes the normal distribution with mean µ and varianceσ 2 . In all simulations throughout this Letter, the natural frequencies are drawn from N (0, 1). Some typicaltime evolutions of the order parameter R representing
the degree of synchrony is shown in Fig. 1(c), in which
the above system starts from diﬀerent initial conditions.
The order parameter R is deﬁned by
R exp(iψ) =1 ∑exp(iφj ),
N j

(2)

where ψ is the average phase associated with the order parameter. As illustrated in Fig. 1(c), the system starting from a completely uniform initial distribution remains desynchronized, while the system with nonuniform initial distribution can go to various synchronized states in a similar way as in Fig. 1(a). Two numerical simulations shown in Fig. 1 suggest that the system
containing three-body interaction can exhibit multistablebehaviors in a structurally stable manner.
To investigate these behaviors analytically, we here impose three assumptions which do not spoil the essence of
the above dynamical behaviors. First we assume that thephase coupling functions are identical for all oscillators,
that is, Γijk (φji , φki ) = Γ(φji , φki )/N 2 . Second, without
loss of generality, we can assume that the phase coupling

function is symmetric, that is, Γ(φji , φki ) = Γ(φki , φji ),because replacing the asymmetric coupling Γasym (x, y)with symmetric coupling Γsym (x, y) = [Γasym (x, y) +Γasym (y, x)]/2 does not change the dynamics. The
last assumption is that inverting the phases of oscillators inverts the sign of forces among them, that is,Γ(φji , φki ) = −Γ(φij , φik ). Although this seems a rather
tight constraint, this antisymmetricity is a propertyof the classical two-body coupling function Γ(φji ) =sin φji . We conﬁrmed that the system under these constraints could exhibit the qualitatively same behavior as
in Fig. 1(c). Finally, we note that, owing to the 2πperiodicity, Γ can be approximated by the ﬁnite Fourier
series Γ(φji , φki ) = K2 (sin φji +sin φki )/2+K20 (sin 2φji +sin 2φki )/2 + K3 (sin φji cos φki + cos φji sin φki ), where
K2 , K20 and K3 are constants. Thus, the dynamics ofglobally-coupled phase oscillators with this type of threebody coupling is given byφ̇i = ωi ++1 ∑(K2 sin φji + K20 sin 2φji )
N j

2K3 ∑
sin φji cos φki .
N2

(3)

j,k

We further simplify this model equation to make it analytically tractable. Using order parameter R and setting
K2 = 0, K30 = 0, and K3 = K, we obtain the equationsof dynamics with pure three-body interaction,θ˙i = ωi − KR2 sin 2θi ,(4)where θi = φi − ψ is the relative phase of oscillator i
to the average phase ψ [Eq. (2)]. Because we are using
a co-rotating frame, we may here assume that the average phase ψ is constant. We assume that the frequencyof the average phase ψ equals the mean of the distribution g(ω) of the natural frequency, the standard normal
distribution. This assumption simpliﬁes the equations
to be derived, and, in addition, the solution of the derived self-consistent equation ﬁts substantially well with
the numerical results, although this assumption may nothold in some cases.
Numerical simulations of Eq. (4) with N = 10, 000
oscillators and K = 3 from three initial conditions areshown in Fig. 2(a). Order parameter R takes various
values depending on the initial conditions. Synchronizedand desynchronized states coexist in the same parameterregion. The relationship between the natural frequency
ωi and the phase φi is also shown in Fig. 2(b,c,d). Figure 2(c) indicates that oscillators can be phase locked totwo speciﬁc phases. Indeed, an oscillator with natural freωiquency ωi can be phase locked to θi = 12 arcsin KR
2, π+
ωi
1
2
2
arcsin
,
if
−KR
≤
ω
≤
KR
.
On
the
other
hand,
2
i
2
KR
Fig. 2(d) shows that the system with the same parametervalues can exhibit a completely desynchronized state.
If all of the phase-locked oscillators are locked to θi =3(b) 2π

1

φi

(a)

π
0

0.8

(c)

R

φi

0.6

2π

0.4

-4-3-2-1 0 1 2 3 4ωiπ0(d) 2πφi

0.2-4-3-2-1 0 1 2 3 4
ωisolutions r2 > r1 , there exists 0 < α < 1 with which
the largest non-zero solution of R = S[R, K, αq(θ)] issmaller than r2 because S[R, K, αq(θ)] = αS[R, K, q(θ)][Fig. 3(a)]. Hence, to obtain the lowest attainable R, wehave to ﬁnd q(θ) with which R = S[R, K, q(θ)] has only
one non-zero solution. In other words, we ﬁnd the smallest r satisfying S[r, K, q(θ)] = r and S 0 [r, K, q(θ)] = 1 invarying q(θ), where∂
S[R, K, q(θ)]
∂R
R=r
∫ π/4
q(θ)W (θ, r, K)C(θ, r, K) dθ,= 4Kr

S 0 [r, K, q(θ)] ≡

π

0
0.1

1

10

0

100 1000 10000-4-3-2-1 0 1 2 3 4
ωi

t

−π/4

FIG. 2. (a) Time evolution of the order parameter R fromthree diﬀerent initial conditions in the mean ﬁeld model with
N = 10, 000 and K = 3. (b,c,d) ωi -φi relationship for diﬀerent initial conditions at t = 10, 000.

ωi
1
2 arcsin KR2 , R is given by

(

∫ KR2

R=

cos
−KR2

∫ π/4

= 2KR2

1
ωarcsin2KR2)g(ω) dωcos θ cos 2θg(KR2 sin 2θ) dθ−π/4≡ S(R, K),(5)2

where we used dω/ dθ = 2KR cos 2θ, and assumed thatthe non-phase-locked oscillators do not contribute to thevalue of the order parameter because the distributiong(ω) of natural frequency is the standard normal distribution. The self-consistent equation R = S(R, K) has
a solution R = 0 for any K. In addition, S 0 (0, K) =∂S
∂R |R=0 = 0 suggests that this solution is stable. For
some K, the self-consistent equation has a non-zero solution R = r > 0 or two non-zero solutions r2 > r1 > 0[Fig. 3(a)]. Equation (5) gives the order parameter ofthe system in which all of the phase-locked oscillators
ωitake θi = 12 arcsin KR
2 . Oscillator i, however, can also
ωibe phase locked to θi = π + 21 arcsin KR
2 . Deﬁning n(θ)
as the number of oscillators phase locked to θ, we characterize the distribution of the phase-locked oscillators withthe function q(θ) = [n(θ) − n(θ + π)]/[n(θ) + n(θ + π)].
Note that |q(θ)| ≤ 1. Then, the order parameter of thesystem is given by∫ π/4q(θ)C(θ, R, K) dθ ≡ S[R, K, q(θ)],
R = 2KR2−π/4(6)where C(θ, R, K) = cos θ cos 2θg(KR2 sin 2θ).
The largest attainable R for coupling strength K isgiven by the largest solution r2 of Eq. (5), while the smallest attainable non-zero R for the coupling strength K isgiven by the minimum of the largest positive solution of
the self-consistent equation Eq. (6) over all possible realizations of q(θ). If R = S[R, K, q(θ)] has two non-zeroandW (θ, r, K) = 1 + Kr2 sin 2θ

g 0 (Kr2 sin 2θ)g(Kr2 sin 2θ)[Fig. 3(c)]. To this end, ﬁrst we ﬁx R to r and examine
whether there exists a solution −1 ≤ q(θ) ≤ 1 of theequations S[r, K, q(θ)] = r and S 0 [r, K, q(θ)] = 1. If it
exists, there is a solution −1 ≤ q1 (θ) ≤ 1 of equations
S[r, K, q(θ)] = r and S 0 [r, K, q(θ)] = s1 ≤ 1 [Fig. 3(b),blue line], and there is a solution −1 ≤ q2 (θ) ≤ 1 ofequations S[r, K, q(θ)] = r and S 0 [r, K, q(θ)] = s2 ≥ 1[Fig. 3(b), green line]. Conversely, if q1 (θ) and q2 (θ)
are given, −1 ≤ q(θ) = uq2 (θ) + (1 − u)q1 (θ) ≤ 1,
1
where 0 ≤ u = s1−s
≤ 1, is a solution of the equa2 −s1
tions S[r, K, q(θ)] = r and S 0 [r, K, q(θ)] = 1 [Fig. 3b, redline]. Thus, the existence of q1 (θ) and q2 (θ) which satisfy
S[r, K, q1 (θ)] = r, S 0 [r, K, q1 (θ)] ≤ 1, S[r, K, q2 (θ)] = r,
and S 0 [r, K, q2 (θ)] ≥ 1 is a necessary and suﬃcient condition of the existence of the solution of S[r, K, q(θ)] = r
and S 0 [r, K, q(θ)] = 1. It is suﬃcient for us to calculate
the maximum and the minimum of S 0 [r, K, q(θ)] under
the constraints S[r, K, q(θ)] = r and |q(θ)| ≤ 1.
S[r, K, q(θ)] and S 0 [r, K, q(θ)] have the same domain
of integration, and their integrands diﬀer by a factor ofW (θ, r, K). Hence, the maximum of S 0 [r, K, q(θ)] under
the conditions S[r, K, q(θ)] = r and |q(θ)| ≤ 1 is given
by S 0 [r, K, q2 (θ)] where q2 (θ) = 2Θ[W (θ, r, K) − w2 ] − 1.
Here Θ(x) is the Heaviside function, and w2 is set tosatisfy S[r, K, q2 (θ)] = r. In this case, the phase-locked
oscillators are distributed according to n(θ)/[n(θ)+n(θ +π)] = Θ[W (θ, r, K) − w2 ]. In other words, ﬁrst we adjust w2 to set S[r, K, q2 (θ)] = r [Fig. 3(d)], and next wecheck whether S 0 [r, K, q2 (θ)] is larger than 1 [Fig. 3(e)].
In the same way, we vary w1 to set S[r, K, q1 (θ)] = r,
where q1 (θ) = 2Θ[w1 −W (θ, r, K)]−1, and check whether
S 0 [r, K, q1 (θ)] is smaller than 1.
Thus, we have theoretically obtained the region of theorder parameter which can be achieved by choosing suitable initial conditions [Fig. 3(f), red line]. In this ﬁgure,the dots represent the data from numerical simulations(N = 10, 000). The theoretical results agree with the numerical ones, though several points with K < 3 lie outside of the theoretically derived region. This discrepancy

4

j,k

As Fig. 3(g) shows, as K3 increases, the system ﬁrst startsout exhibiting either a single synchronized or desynchronized state depending on K2 . Then brieﬂy, a small window in which these two states are bistable, appears. Finally, multistable synchronized states, or for smaller K2 ,a coexistence of desynchronized and multistable synchronized states [Fig. 3(g), orange region] corresponding to
the multistability shown in Fig. 2, appears. This implies
that our theoretical result derived with pure three-bodyinteraction is structurally stable and generic.
In this Letter, we have examined the behavior of phaseoscillators with three-body interactions. We have foundthat this system can take an inﬁnite number of synchronized states in a structurally stable manner [Fig. 3(g)].
We have derived the range of the order parameter R thatcan be attained by varying the initial condition. Our results are diﬀerent from the chimera state [9], because in
our model we can continuously control the order parameter of the steady state by choosing the initial condition.
In addition, our model system can be completely incoherent even in the K → ∞ limit (cf. [10]). There remain
several questions to be answered. Three-body interactions in real-world systems and their behavior should becompared to those of the present model. Neurophysiological experiments [11] have shown that some prefrontal
neurons keep their level of activity for several seconds.
It is believed that this persistent activity serves as working memory by encoding an analog quantity in the ﬁring
rate of multistable neuronal networks. Our results suggest the possibility that working memory uses the degreeof synchrony among neurons to encode an analog quantity. Finally, we should systematically investigate various
types of coupling function and the dynamical behavior oncomplex networks [12].
This work was supported by KAKENHI 21700250,
23115512, 19GS0208, 21120002, and 23115511 from
MEXT, and Global COE Program “Center for Frontier
Medicine”, MEXT, Japan.

∗
tanaka.takuma@gmail.com[1] B. A. Loiseau and Y. Nogami, Nucl. Phys. B2, 470(1967); H. P. Büchler et al., Nature Physics 3, 726 (2007).

0.2

0.6

q1(θ)
q2(θ)
uq2(θ)+(1-u)q1(θ)

0.4
0.2
r

0

0

(c)
2

0.2 0.4 0.6 0.8
R
1

0
-2
-4
-6
-8
-10
-π/4

(d)0.4

w2

0.5

w1

0

-0.5
-1
π/4

0
θ

0.2 0.4 0.6 0.8
R

S [r,K,q(θ)]
1

0.35
0.3
0.25
0.2
0.15
0.1
0.05
0
-π/4

(e)0.4

0.5
0

-0.5
-1
π/4

0
θ

(f) 1 Fig. 2(b)
0.8

0

1

q(θ)

0

1

S '[ r,K,q (θ)]
1

0.3
0.2
0.1
0
-0.1
-0.2
-π/4

0.5
0

-0.5
-1
π/4

0
θ

3

(g) Single sync.

Multi-sync.

2.5

Fig. 2(c)q(θ)

S[R,K,q(θ)]

S[R,K,q(θ)]

r1

0.4

0.8

W(θ,r,K)C(θ,r,K)

sin φji cos φki .

r

2
0.6
1.5

0.4

Bistable

1

Multi-sync. and desync.
1

0.2

R

N2

αq(θ)

0.6

K2

j

sin φji +

0.8

q(θ)
C(θ,r,K)

N

2K3 ∑

(b) 1

q(θ)

W(θ,r,K)φ̇i = ωi +

K2 ∑

r2

(a) 1

R

may be because the system size is too small or because
we assumed that the frequency of the average phase coincides with the mean of the distribution g(ω).
Finally, we should remark that interactions in realworld systems generally contain not only three-body but
also two-body interactions. We thus examine the behavior of the system described by

0.5

Fig. 2(d)0Desync.
0

0
0

1

2

K

3

4

5

0

0.5

1

1.5
K3

2

2.5

3

FIG. 3. (color online). (a) For some q(θ), the self-consistentequation R = S[R, K, q(θ)] have two solutions, r1 < r2 (brownline). Setting α appropriately makes R = S[R, K, αq(θ)]
have only one solution R = r < r2 (red line). Note thatS 0 [r, K, αq(θ)] = 1. (b) If we have S[r, K, q1 (θ)] = r andS 0 [r, K, q1 (θ)] = s1 ≤ 1 (blue line) and S[r, K, q2 (θ)] = r
and S 0 [r, K, q2 (θ)] = s2 ≥ 1 (green line), we can obtainq(θ) with which S[r, K, q(θ)] = r and S 0 [r, K, q(θ)] = 1
hold (red line). (c,d,e) q2 (θ) which maximizes S 0 [r, K, q(θ)](e) under the constraint S[r, K, q(θ)] = r (d) is given by
q2 (θ) = 2Θ[W (θ, r, K) − w2 ] − 1 (green line) where w2 is set
to satisfy S[r, K, q2 (θ)] = r (c). Under the same constraint,
S 0 [r, K, q(θ)] is minimized by q1 (θ) = 2Θ[w1 − W (θ, r, K)] − 1(blue line) where S[r, K, q1 (θ)] = r. (f) Attainable region of
the order parameter R (orange region). Note that the incoherent state R = 0 is stable for any K (red line). (g) Phase
diagram of the system of phase oscillators when the strengthof two-body and three-body interactions are changed. The
symbol in each region is a schematic representation of the attainable values of the order parameter R. Gray lines represent
the range of R from 0 to 1. The attainable values and rangesof R are indicated by black circles and boxes, respectively.

[2] Y. Shoham and K. Leyton-Brown, Multiagent systems: algorithmic, game-theoretic, and logical foundations (Cambridge University Press, 2009).
[3] K. S. Hsu et al., Brain Res. 690, 264 (1995); A. G. Carter
et al., J. Neurosci. 27, 8967 (2007); M. E. Larkum et al.,Science 325, 756 (2009).
[4] G. M. Shepherd, ed., The synaptic organization of thebrain (Oxford University Press, 2004); J. J. O’Connor
et al., Nature 367, 557 (1994); F. Saitow et al., J. Neurosci. 25, 2108 (2005).
[5] P. Dayan and L. F. Abbott, Theoretical Neuroscience(MIT Press, 2001).
[6] This model is described by V̇i = −gNa m3i hi (Vi − ENa ) −5
gK n4i (Vi − EK ) − gLeak (Vi − ELeak ) + Isyn + Ii , I˙syn =PP
2π[Ti0 (t)−Tk0 (t)]I− τsyn+ Ne, where Ii
j,n δ(t, Tj,n )
k cos
τ2
1
is the baseline input current of neuron i, Ti,n is then-th spike timing of neuron i, Ti0 (t) is the last spiketime of neuron i at time t, τ1 = 3 msec is the decay
time constant of synaptic current, e = 40 µA/mm2
is the maximum amplitude of synaptic current, andτ2 = 15 msec is the time scale of short-term plasticity.
The dynamics of gate variables follow those of the original Hodgkin-Huxley model [5]. Baseline input follows
Ii = 13 + 5(i/N − 1/2)3 + 5(i/N − 1/2).
[7] T. Danino et al., Nature 463, 326 (2010); J. Buck, Q.
Rev. Biol. 63, 265 (1988); S. H. Strogatz et al., Nature
438, 43 (2005).

[8] Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence (Springer-Verlag, 1984); F. C. Hoppensteadt
and E. M. Izhikevich, Weakly connected neural networks(Springer Verlag, 1997); A. T. Winfree, The geometry of
biological time (Springer Verlag, 2001); S. H. Strogatz,Physica D 143, 1 (2000); J. A. Acebrón et al., Rev. Mod.
Phys. 77, 137 (2005); B. Ermentrout, Neural Computation 8, 979 (1996).
[9] Y. Kuramoto and D. Battogtokh, Nonlinear Phenom.
Complex Syst. 5, 380 (2002); D. M. Abrams and S. H.
Strogatz, Phys. Rev. Lett. 93, 174102 (2004).
[10] H. Daido, Phys. Rev. Lett. 77, 1406 (1996).
[11] R. Romo et al., Nature 399, 470 (1999).
[12] S. Boccaletti et al., Phys. Rep. 424, 175 (2006).