PHYSICAL REVIEW E 81, 026203 共2010兲

Dynamics of coupled adaptive elements: Bursting and intermittent oscillations
generated by frustration in networks
Masayo Inoue1 and Kunihiko Kaneko1,2

1

Department of Basic Science, Graduate School of Arts and Sciences, University of Tokyo, 3-8-1 Komaba,
Meguro-ku, Tokyo 153-8902, Japan
2
ERATO Complex Systems Biology Project, JST, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan
共Received 7 July 2009; revised manuscript received 26 September 2009; published 4 February 2010兲

Adaptation to environmental change is a common property of biological systems. Cells initially respond toexternal changes in the environment, but after some time, they regain their original state. By considering anelement consisting of two variables that show such adaptation dynamics, we studied a coupled dynamicalsystem containing such elements to examine the diverse dynamics in the system and classified the behaviors on
the basis of the network structure that determined the interaction among elements. For a system with twoelements, two types of behaviors, perfect adaptation and simple oscillation, were observed. For a system withthree elements, in addition to these two types of dynamics, rapid-burst-type oscillation and a slow cycle werediscovered; depending on the initial conditions, these two types of dynamics coexisted. These behaviors are a
result of the characteristic dynamics of each element, i.e., fast response and slow adaptation processes. The
behaviors depend on the network structure 共in specific, a combination of positive or negative feedback among
elements兲. Cooperativity among elements due to a positive feedback loop leads to simple oscillation, whereas
frustration involving alternating positive and negative interactions among elements leads to the coexistence of
rapid bursting oscillation and a slow cycle. These behaviors are classified on the basis of the frustration indicesdefined by the network structure. The period of the slow cycle is much longer than the original adaptation timescale, while the burst-type oscillation is a continued response that does not involve any adaptation. We briefly
discuss the universal applicability of our results to a network of a larger number of elements and their possible
relevance to biological systems.
DOI: 10.1103/PhysRevE.81.026203

PACS number共s兲: 05.45.⫺a, 87.10.⫺e, 82.39.Rt

I. INTRODUCTION

Adaptation is a common property of biological organisms.
Individual organisms initially sense and respond to externalchanges in their environment; however, after some time, they
become habituated to the new environment, and hence, theresponse disappears. Adaptation can occur on a time scale ofgenerations, but it also progresses on much shorter timescales of seconds, minutes, or hours, i.e., within the lifetimeof an individual organism. Here, we focus on the adaptive
behavior on such short time scales; such behavior is generally observed not only in multicellular organisms but also inunicellular organisms like E.coli or paramecium and is regarded as a general property of a cell.
The study of cellular adaptation behavior has been pioneered by Koshland 关1,2兴 and Oosawa 关3,4兴. Koshland classified adaptation as absolute 共perfect兲 and partial adaptation.
Environmental change has an impact on intracellular statessuch as gene expression patterns through reaction networks,
for example, by signal transduction. In absolute adaptation,
variables representing the concentration of some chemicals
start to change in response to environmental change, andthen they regain their original values even though the environmental change is not reversed. Indeed, several examplesof simple biochemical circuits have been proposed to represent such adaptive behaviors 关5–7兴. Characteristic features as
dynamical systems are elucidated and general conditions in
the topology in gene regulation network to show such adaptive dynamics have been studied 关8,9兴. In addition to theoretical interest, experimental verification of such adaptation
1539-3755/2010/81共2兲/026203共14兲

processes has gained much attention, and recently there areextensive studies at a quantitative level 关10–13兴.
Adaptive response is not limited to a specific variable of abiological system. Not one but many variables regain their
original values within a certain time after an environmental
change has occurred. This type of response has been studied
in the microarray analysis of gene expression levels 关12,13兴.
Often, some partial units in a whole biological system have atendency to maintain their state, independent of the external
conditions—this is known as homeostasis. On the other
hand, they also have to respond to an external change inorder to survive. Hence, the adaptive behavior is expected to
be an inherent process that involves many partial units of abiological system. Then, the dynamic behavior of the wholesystem is determined by the behavior of several adaptation
units that interact with each other.
Now, in addition to revealing the processes and mechanisms in an adaptation unit, i.e., an adaptation module, it isalso important to study a system consisting of several adaptation units that are not isolated but interact with other units.
Here, such a unit of adaptive behavior can be a cell in acolony of interacting cells or in a multicellular organism. In
this case, the interaction occurs through chemical signals resulting in cell-cell communication. Alternatively, a certain
part of a global intracellular reaction network can function asa unit to exhibit adaptive behavior. In this case, each unit
interacts through reaction paths in the network. For example,
a signal transduction network or metabolic network is highly
complicated, involving many feedback and feedforward interactions or multistep reactions. Each adaptive or nonadap-

026203-1©2010 The American Physical Society

PHYSICAL REVIEW E 81, 026203 共2010兲

MASAYO INOUE AND KUNIHIKO KANEKO

tive unit is embedded within a cell and receives inputs not
only from outside the cell but also from other units within
the intracellular network. Such interaction between units
should play an important role in the global function of a cell.
In this study, we examine the behavior of a coupled system consisting of elements that exhibit adaptive behavior. If
each element is a module in a chemical reaction network or areaction motif, the coupled adaptive system represents a global intracellular reaction network. If each element represents
a cell, the system represents the cooperative behavior of cells
within a multicellular system.
As for dynamical systems with many degrees of freedom,
coupled oscillators 关14兴 and coupled chaotic systems 共suchas coupled maps兲 关15,16兴 have been extensively studied overdecades. Several concepts such as synchronization transition,clustering chaotic itinerancy, and collective chaos have beendeveloped in these studies. In contrast, systems of coupledadaptive elements have not been systematically studied thus
far. However, such a study may help develop concepts andreveal nontrivial interesting behaviors. In fact, dynamics thatfacilitate perfect adaptation have the following characteristics that distinguish them from other systems: when a specific parameter of the system is changed, the response to thischange appears first as a change in the values of the variables, but later, some variables regain their original values
independent of the parameter values. This independence results in the imposition of some constraints. Further, theabove behavior implies the existence of two time scales: onefor response and the other for relaxation to the original value.
Therefore, the manner of return to the original values and theexistence of the two time scales result in nontrivial dynamics
in a coupled system. Here, we classify the behaviors ofcoupled adaptive systems on the basis of the interaction network structure among elements. For the classification, we
will introduce the measure of frustration to characterize theinteraction among elements.
The paper is organized as follows. In Sec. II, we introduce
a model of coupled adaptive elements. In Sec. III, we present
results for the two-element case. The three-element case isstudied in Sec. IV. In this case, the coexistence of a slow
cycle and fast bursting oscillations is observed when frustration is taken into account. Extensions to cases with a larger
number of elements will also be discussed.
II. MODEL

Here, we introduce a model of coupled adaptive elements.
In this model, elements that show adaptive behavior to an
input interact with each other. Here, adaptation refers to dynamics by which some variables regain their original value,
independent of the value of the external inputs.
Now, we consider a minimal system of such adaptation.
When the value of a constant external input S is changed,
first, a state variable u changes in response to this input,
however, after some time, u regains its original value independent of the value of S. This system is a type of excitablesystem in which a state variable changes in response to a
stimulus and then returns to the prestimulated state. The simplest way to construct such dynamics is to introduce anothervariable 共v兲 that “absorbs” the changes in S. Hence, we consider the following two-variable dynamical system:du= f共u, v ;S兲,dt
dv
= g共u, v ;S兲,dt

共1兲

where the fixed point solution uⴱ , vⴱ given by f共uⴱ , vⴱ ; S兲
= 0 , g共uⴱ , vⴱ ; S兲 = 0 should be set to be stable. The above
condition is satisfied if uⴱ is invariant against the change inS, even though f increases with S. Then, the postulate for the
response to S and adaptation is satisfied because u first increases with S and then regains the value uⴱ.
As an example of such adaptive dynamics, we consider achemical reaction. In the reaction process, time evolution isgoverned by a set of rate equations in which variables u andv represent chemical concentrations. Indeed, there have been
several models that exhibit adaptive dynamics 关6,7,17,18兴. In
this paper, we adopt a system obtained by adiabaticallyeliminating a variable from a model originally introduced by
Levchenko and Iglesias 关19兴. The model was introduced torepresent adaptation process in biochemical signal transduction pathways underlying gradient sensing in amoebae andneutrophils. In this simplified model, a chemical U is either
active or inactive: u represents the concentration of thechemical U in the active state and ū represents that in theinactive state. We set the sum of u and ū to 1共u + ū = 1兲. An
external input 共S兲 activates U by catalyzing the process fromū to u, while it is also involved in the synthesis of anotherchemical V, which acts as a catalyst to deactivate U, asshown in Fig. 1. In fact, transformation between active andinactive forms of molecules is commonly observed in biological processes such as methylation 关6兴, phosphorylation,and so forth, while the activation or inactivation are generally controlled by external signals as well as internal states.
Hence the model in Fig. 1 is one of the prototype models forbiological adaptation.
When the degradation of V is also included, the kinetics isrepresented as
du S共1 − u兲 − uv
,
=
dt
␩

dv S − v
.
=
dt
␮

共2兲

We rescale time as t̃ = t / ␩ in all the following discussions,˜ = ␮ / ␩, Eq. 共2兲 is rewritten as
and by setting ␮dudt̃= S共1 − u兲 − uv,

dv
dt̃

=

S−v
.
˜
␮

共3兲

This equation has a stable fixed point vⴱ = S, uⴱ = ūⴱ = 0.5,
which is independent of S. When the input S increases 共decreases兲, the concentration of U in the active 共inactive兲 stateincreases; therefore, u increases 共decreases兲 from its steadystate value 共uⴱ = ūⴱ = 0.5兲 according to the change in S, but
after some time, it regains the original value uⴱ as v relaxes
to the value S. The dynamics represented by Eq. 共3兲 have twotime scales. The time scales are defined by the reciprocals of
the two eigenvalues of the dynamics. The eigenvalues are

026203-2

DYNAMICS OF COUPLED ADAPTIVE ELEMENTS: …

(a)

PHYSICAL REVIEW E 81, 026203 共2010兲

S

u

u

v
(b)
1
0.9
10

0.8

v(t), S(t)

0.7
u(t)

0.6
0.5
5

0.4
0.3
0.2
0.1
0

0

20

40

60

80
time

100

120

140

0

FIG. 1. 共Color online兲 共a兲 Schematic of a reaction system represented by Eq. 共2兲. The area enclosed by the dashed line contains
one element. 共b兲 A typical example of a response that Eq. 共2兲 obeys;
u 共red ⫻ point, plotted on the vertical axis to the left兲 represents
adaptive dynamics and v 共dashed blue line, plotted on the verticalaxis to the right兲 represents relaxation according to the change in S
共chained cyan 共light gray兲 line, plotted on the vertical axis to theright兲.

calculated by linear stability analysis around the steady state.
The first time scale is given by ␶s = 1 / 共2S兲, which corresponds to the response time to the change in S, and the sec˜ , which is the relaxation time to the
ond is given by ␶a = ␮
original value 共i.e., adaptation兲. Here, we postulate faster re˜ . Seesponse and slower adaptation so that ␶s = 1 / 共2S兲 ⬍ ␶a = ␮
关3,17兴 for the relevance of this condition to biological adaptation.
Now, we consider a coupled system of elements each ofwhich follows the above reaction mechanism with the two
chemicals 共u , v兲 that temporally evolves according to Eq. 共3兲.
We assume that the interaction between two elements j and ican be modeled by using the input term S, i.e., an output u jfrom an element j triggers an input to the element i共Si兲. This
process u1 , u2 , ¯ → Si is represented by
Si = h共u1,u2, . . .兲.

共4兲

If each element i represents a cell, this interaction represents cell-cell communication through signal molecules. In
fact, cells mutually transmit and receive chemical signals towhich each cell shows adaptive dynamics. For example, cellular slime mold, Dictyostelium, communicate with cyclic
AMP 关5,20兴. Response to cAMP of each cell shows adaptive
dynamics of the above type, while with cell-cell interaction,
oscillatory and relaxation dynamics in intracellular cAMP
concentrations are recently measured 关21兴. Another example
in coupled adaptive systems is seen in neural networks where
each neuron adapts to external signal and ceases to respondafter some time, the neurons form a network with activationand inhibition of each other. On the other hand, if each element i represents a set of chemical components in an intracellular reaction, the interaction represents signal transduction in an intracellular reaction network. Indeed biologicalnetworks often consist of several molecules with active andinactive forms, as seen, for example, in a cascade reaction ofkinase-kinase-… enzymes 关22兴.
To choose a specific type of interaction, we take the following two facts into account. First, in a biological reaction,
input-output relations often take a sigmoidal form with somethreshold, as in an enzyme reaction in which relations aremodeled by a Michaelis-Menten type or Hill function. Further, cell-cell interactions can often be represented by such
threshold functions, e.g., in binding at receptors or in electricinteraction in nerve cells. Hence, we adopt a such sigmoidal
reaction to represent interactions between elements. The input S is transmitted according to the concentration of u in theactive state. When the balanced state 关u = ū = uⴱ共=0.5兲兴 is set
to the threshold, the signal is transmitted depending on thesign of u − uⴱ. Hence, each element receives a signal as afunction of 兺 jCij共tanh兵␬关u j共t兲 − uⴱ兴其兲, where Cij represents aninteraction factor from element j to i, and the element of theinteraction matrix Cij is 0, 1, or −1 depending on whether thesignal is nonexistent, excitatory, or inhibitory. We exclude
self-catalysis and define Cii = 0. Second, catalytic reactions
often occur in succession to relay signal transduction, andtherefore, the effects of the components are cascaded. Notealso that the change in the signal concentration often occurs
on the “logarithmic scale.” Without loss of generality, we canset S = 1 when there is no input from other elements, S= Smax ⬎ 1 if the received signal is excitatory, and S = Smin−1
= Smax
if it is inhibitory; thus, we obtain冋兺

Si = exp册␴
Cij共tanh兵␬关u j共t兲 − uⴱ兴其兲 ,
Ni j

共5兲

where Ni is the number of inputs that the element i receives,
i.e., Ni = 兺Nj=1兩Cij兩, and thus, comparison with cases involving
different networks with a different number of paths isstraightforward. With this form, Smax is given by the couplingstrength ␴ as Smax = e␴. Without loss of generality, we can set
Smax = 10 and ␴ = ln 10 in order to limit the range of S to10−1 – 10+1; then, the form of the interaction isⴱ

Si = 10共1/Ni兲兺 jCij共tanh兵␬共u j共t兲−u 兲其兲 .

共6兲

We use ␬ = 50 unless otherwise mentioned. The responsetime scale of each element is given by ␶s = 1 / 共2S兲 and that˜ . Setting ␮˜ = 10, we can use ␶s
for adaptation by ␶a = ␮= 0.05– 5 and ␶a = 10 in order to satisfy the condition for fast
response and slow adaptation 共␶s ⬍ ␶a兲. There is no qualitative change in our results as long as the threshold-type interaction is chosen and ␶s ⬍ ␶a irrespective of the parametervalue. Even if the exponential coupling form as in Eq. 共5兲 is
not adopted, most of the behaviors to be reported here are
obtained.
In the present model, each element communicates withother elements according to its state, i.e., depending onwhether it is active 共u ⬎ uⴱ兲 or inactive 共u ⬍ uⴱ兲. If Cij ⬎ 0, an

026203-3

PHYSICAL REVIEW E 81, 026203 共2010兲

MASAYO INOUE AND KUNIHIKO KANEKO

inactivate

S i activate(a)10

positive interaction

1

negative interaction
0.1
0

0.1

0.2 0.3

0.4

inactive state(b)

0.5 0.6

u j (t)0.7 0.8

0.9

1active state

1
5

2

6
4

3

FIG. 2. 共Color online兲 共a兲 Interaction function represented by
Eq. 共6兲. The solid red line represents a positive interaction and thedotted blue line represents a negative interaction. 共b兲 Conceptual
scheme of a coupled system. Each numbered circle represents anelement with 共u , v兲 adaptive dynamics. The solid red and dottedblue arrows represent positive and negative interactions,
respectively.

active 共inactive兲 element j activates 共inhibits兲 the following
element i, respectively, and the converse is true for Cij ⬍ 0
共Fig. 2兲. We use the above normalization 共summation Ni兲 sothat elements receiving a positive signal have identical input
values and identical response time scales, whereas those receiving a negative signal have smaller input values 共equal forall such signals兲 corresponding to a longer time scale.

If one interaction coefficient is positive and the other isnegative 共C12C21 ⬍ 0兲, the system exhibits perfect adaptation.
Without loss of generality, we choose C21 ⬎ 0 and C12 ⬍ 0.
When element 1 is activated 共u1 ⬎ uⴱ兲, element 2 is also activated by positive interaction 共C21 ⬎ 0兲. On the other hand,negative interaction with element 2共C12 ⬍ 0兲 inactivates element 1. Since the inputs received by the elements have opposite signs, excitatory output is not sustained. The systemthen undergoes damped oscillations and eventually exhibits
perfect adaptation.
In contrast, the system shows a sustained oscillation if
C12C21 ⬎ 0, for a range of parameter values to be analyzed
below. Depending on the sign of C12, elements 1 and 2 oscillate in phase 共if both C12 and C21 are positive兲 or in antiphase 共if both C12 and C21 are negative兲 between the active
and inactive states, as shown in Fig. 3. The period of thelimit cycle is ⬃␶a.
To understand the different behaviors for different networks, we first perform linear stability analysis around thesteady state 共ui = uⴱ = 0.5 and vi = vⴱ = 1兲. Setting ui = uⴱ + ␦ui
and vi = vⴱ + ␦vi, the linearized form of Eq. 共3兲 is given by
d␦ui

1
␬␴
= − 2␦ui − ␦vi +兺 Cij␦u j ,
2
2Ni j
dt̃
d␦vi
dt̃

1
␬␴
␦vi +
兺 Cij␦u j .
˜
˜ j
␮
2Ni␮

共7兲

Defining ␣ = 共 ␬2␴ 兲2 N121N212 = 共 ␬2␴ 兲2C12C21, the four eigenvalues ␭ in the two-element case are given as solutions of thecharacteristic equation
C C

冉 冊
␭+

yielding冉− 2+III. RESULTS FOR TWO-ELEMENT CASE

␭=First, we present results for the case involving two elements 共1 and 2兲. When one of the Cij values was 0, thesystem simply converged to the fixed point with uⴱ = 0.5 andvⴱ = 1. When one of the elements does not receive any signal,
it returns to the original fixed point, and then, the other element does so as well. In addition to such a trivial case, there
are three cases for different C12 and C21: positive-positive,
positive-negative, and negative-negative interactions. Depending on the type of interaction, two types of behaviors,
perfect adaptation and oscillation, are observed 共Fig. 3兲. In
perfect adaptation, each element eventually shows perfectadaptation, and therefore, the overall system also adapts perfectly and regains the original steady state. In the case of
oscillation, each element cannot adapt to the steady state andrepeats the cycle of response and adaptation. The output ofelement 1 causes a change in the input of element 2 according to the interaction function 共6兲, and element 2 shows response and adaptation. This output of element 2 in turncauses a change in S1. Thus, response and adaptation are
repeated continuously.

=−

1 2
共␭ + 2兲2 − ␣␭2 = 0,
˜
␮

冊 冑冉

1
⫾ 冑␣ ⫾
˜
␮

2

2+

共8兲

冊

2
1
8
⫾ 冑␣ −
˜
˜
␮
␮

.

共9兲

˜ ⬎ 0 because ␮˜ means adaptive
Here we limit the case ␮˜ values,
time scale in this model. According to the ␣ and ␮
the four eigenvalues are classified as follows: 共a兲 ␣ ⬍ 0, twopairs of complex conjugate eigenvalues with the negative
real part, 共b兲 0 ⬍ ␣ ⬍ 共冑2 − 冑1␮˜ 兲4, four negative real eigenvalues, 共c兲 共冑2 − 冑1␮˜ 兲4 ⬍ ␣ ⬍ 共2 + ␮˜1 兲2, two negative real eigenvalues and a pair of complex conjugate eigenvalues with thenegative real part, 共d兲 共2 + ␮˜1 兲2 ⬍ ␣ ⬍ 共冑2 + 冑1␮˜ 兲4, two negative
real eigenvalues and a pair of complex conjugate eigenvalues
with the positive real part, 共e兲 共冑2 + 冑1␮˜ 兲4 ⬍ ␣, two negative
real eigenvalues and two positive real eigenvalues 共Fig. 4兲.
Therefore, the transition from the perfect adaptation to the˜ ⱖ 0兲. If
simple oscillation occurs when ␣ = 共2 + ␮˜1 兲2共␮
C12C21 ⬍ 0, the system belongs to the case 共a兲, and thus, theoriginal steady state is stable and perfect adaptation is guaranteed. In contrast to the single-element case, these eigenvalues are complex, which leads to damped oscillation, as
shown in Fig. 3共a兲. When C12C21 ⬎ 0, the system shows Hopf

026203-4

DYNAMICS OF COUPLED ADAPTIVE ELEMENTS: …

PHYSICAL REVIEW E 81, 026203 共2010兲

C21
1

2

C12
(a) 1

10

0.6

v(t)

u(t)

0.8

1

0.4
0.2
0

0

1

2

time

3

4

5

0.1

(b) 1

0

1

0

50

2

time

3

4

5

200

250

300

200

250

300

10

0.6

v(t)

u(t)

0.8

1

0.4
0.2
0

0.1
0

50

100

150

200

250

300

100

time

150

time

(c) 1

10

0.6

v(t)

u(t)

0.8

1

0.4
0.2
0

0.1
0

50

100

150

200

250

300

0

time

50

100

150

time

FIG. 3. 共Color online兲 Time series of u 共left兲 and v 共right兲 for a two-element system. 共a兲 共C12 ⫻ C21 ⬍ 0兲 corresponds to perfect adaptation
behavior. 共b兲 共C12 ⬎ 0 , C21 ⬎ 0兲 and 共c兲 共C12 ⬍ 0 , C21 ⬍ 0兲 correspond to oscillations in phase and in antiphase, respectively. In each picture,
the dotted red line shows the dynamics of element 1; the solid blue line, the dynamics of element 2.

bifurcation to create a limit cycle, as ␣ is increased from theregion 共c兲 to 共d兲. For the parameters adopted here, the system
is at the region 共e兲 and the original fixed point is a saddle,
whereas the limit cycle still remains. This generation of alimit cycle is independent of the sign of C12 as long as
C12C21 ⬎ 0.
The sign of C12 determines whether the limit cycle isin-phase and antiphase oscillations between the two elements
共as shown in Fig. 3兲. This is explained as follows. Because ofa Z2 symmetry of the system, only the two types of periodic
solutions are possible according to the sign of C12 关23兴. First,

as observed in Fig. 3共b兲 共C12 ⬎ 0 , C21 ⬎ 0兲, the oscillationsof two cells are synchronized, i.e., u1共t兲 = u2共t兲 for all time t.
This solution is obtained by setting u1 = u2. Second, in Fig.
3共c兲 共C12 ⬍ 0 , C21 ⬍ 0兲, such in-phase solution does not exist. Following the argument of Z2 symmetry, only the timeperiodic solutions with a half-period out of phase exist, i.e.,
u2共t兲 = u1共t + T / 2兲, where T is the period of the limit cycle.
IV. RESULTS FOR THREE-ELEMENT CASE

For a three-element system, there are six directed interactions between elements that are given by Cij, which takes the026203-5

PHYSICAL REVIEW E 81, 026203 共2010兲

MASAYO INOUE AND KUNIHIKO KANEKO

1~
µ
(b)

0.6

v(t)

0.4

(c)

1

0.2

2

0

(a)

(d)

(b)

1

0

10

time

20

0.1

30

0

1

10

20

time

30

10

20

25

v(t)

0.4

1

0.2
0

FIG. 4. 共Color online兲 Phase diagram of the two-element systemwith regards to the eigenvalues of the Jacobian matrix around the
˜ 共ordinate
fixed point. Plotted against ␣ 共abscissa axis兲 and 1 / ␮˜ areaxis兲. The definitions of area 共a兲 ⬃ 共e兲 and the indices ␣ and ␮
given in the text. The shaded area with slashes shows existence ofcomplex eigenvalues and area with dashed horizontal lines shows
existence of eigenvalues with positive real part.

100

200

time

300

400

0.1

500

(c) 1

0

100

200

time

300

400

500

10

0.8

u(t)

0.6
0.4

1

0.2
0

0

500

1000
time

1500

2000

1

(d)

0.1

0

500

1000
time

1500

2000

0

500

1000
time

1500

2000

10

0.8

u(t)

0.6
0.4

1

0.2
0

0

500

1000
time

1500

2000

(e) 1

0.1

10

0.8

u(t)

value 1, −1, or 0. Hence, the total number of possible networks is 36 = 729. By disregarding the cases that are identicaldue to symmetry and those in which one or two elements are
disconnected, we obtain 78 possible networks. For example,
if the third element is driven by only one of the other twoelements and does not transmit to other elements, this obviously results in the same behavior as the two-element case,and hence, such networks are excluded. We have studied allof these 78 cases to find three distinct behaviors: 共i兲 perfectadaptation, 共ii兲 simple oscillation 共as in the two-elementcase兲, and 共iii兲 coexistence of rapid-burst-type oscillation anda slow limit cycle 共Fig. 5兲. Among the 78 systems, 7 networks show perfect adaptation and 46 show simple oscillation. These two cases are basically understood as straightforward extensions of the corresponding two-element cases.
The distinction between in-phase and antiphase oscillations
are also explained, based on the analysis on the symmetry
关24兴. The remaining 25 networks show a nontrivial behavior
that first appears in the three-element case 共Fig. 6兲.
There are two distinct types of dynamics depending on
the initial conditions. In the rapid-burst-type oscillation, eachvi tends toward a fixed value vifix, with tiny-amplitude oscillations around the value, while u components show a largeamplitude oscillation with a period on the order of ␶s. Notethat the values vifix are distinct from the original fixed pointvalue vⴱ = 1. In the slow cycle, u components approach thevalues of the original fixed point, remain close to these values for some time, and then depart from it. This process isrepeated periodically on a time scale much longer than ␶a
共the adaptation time scale兲.
Now, we analyze the rapid-burst-type oscillation. In thiscase, the vi’s stay almost constant and no adaptation occurs
in ui’s, which oscillate on the order of the fast response timescale. The existence of the rapid-burst-type oscillation is as˜ tends to infinsociated with the time-scale difference. As ␮
ity, or in other words, as the time-scale ratio ␶s / ␶a between uand v tends to zero, the influence of u on v through S is
averaged out because u oscillates too fast, and hence, theoscillation in v disappears. In fact, we have numerically confirmed that the amplitude of oscillation in the v’s vanishes as

0

v(t)

15

0.6

v(t)

10

α

0.6

v(t)

4 5

u(t)

0.8

(e)
0

10

0.8

u(t)

3

1

(a)

0.4
0.2
0
1999 1999.2

1

~

~

∆t = 1
1999.4∆t = 1

1999.6

1999.8

2000

0.1
1999

time

1999.2 1999.4 1999.6 1999.8
time

2000

C21
C12

1

2

C31
C13

C23
3

C32

FIG. 5. 共Color online兲 Time series of u 共left兲 and v 共right兲 forthe three-element system. Panel 共a兲 共C12 = −1 , C13 = −1 , C21= 1 , C23 = −1 , C31 = 1 , C32 = 1兲 corresponds to perfect adaptation
behavior; 共b兲 共C12 = −1 , C13 = 0 , C21 = 0 , C23 = 1 , C31 = −1 , C32= 1兲 corresponds to an example of simple oscillation. 共c兲 共C12
= −1 , C13 = −1 , C21 = −1 , C23 = 1 , C31 = 1 , C32 = 0兲 and 共d兲 共C12
= −1 , C13 = −1 , C21 = 0 , C23 = 1 , C31 = 1 , C32 = −1兲 correspond toa slow cycle and 共e兲 关same network as in 共d兲兴 corresponds to a rapid
oscillation. Note the difference in the scales of the abscissa fordifferent figures. In each picture, the dotted red line shows the dynamics of element 1; the solid cyan line, the dynamics of element 2;
and the dashed blue line, the dynamics of element 3.

˜ increases.
␮˜Now, we consider this limiting behavior at ␶s / ␶a ⬃ 1 / ␮→ 0. From the condition dv / dt̃ = 0 in Eq. 共3兲, vifix is proved to
be a long-time average of Si, i.e.,

026203-6

DYNAMICS OF COUPLED ADAPTIVE ELEMENTS: …

PHYSICAL REVIEW E 81, 026203 共2010兲

1

0.7

0.8

u3 0.60.6
u(t) , du/dt = 0

0.4
0.2
0

1

0.8

u2

0.6

0.4

0.2

0

0.2

0.4

0.6

0.8

1

u1

0.5
0.4
0.3

10

0.2
1000

v3
10

v1
1
0.1
10

v2

1

0.1

FIG. 6. 共Color online兲 Trajectories in the phase space ofu1 , u2 , u3 共up兲 and of v1 , v2 , v3 共down兲 for three-element systems oftype 共i兲–共iii兲. Each orbit is drawn after neglecting a sufficiently longtransient time. The fixed point 共black⫻ point兲 represents perfectadaptation and the limit cycle shown by the dotted blue line is an
example of simple oscillation. An example of rapid oscillation isplotted as a bold dashed red line for u1 , u2 , u3 and as a + point for
v1 , v2 , v3, as the oscillation amplitude in v is negligibly small. Notethat this fixed value of the v’s differs greatly from the original fixed
point that corresponds to perfect adaptation. The solid cyan 共lightgray兲 line for plots of both u and v denotes a slow cycle. The results
for 共a兲, 共b兲, 共d兲, and 共e兲 in Fig. 5 are plotted.

1
T→⬁ T

vifix = lim冕T

Si共t̃兲dt̃ = 具Si典⬁

共10兲

0

共Fig. 7兲. Because dv / dt̃ = 0, no adaptation occurs and thesystem is driven only by u variables with the fast responsetime scale ␶s⬘ = 1 / 共S + v fix兲, which is obtained fromdt̃

= S共1 − u兲 − uv fix = − 共S + v fix兲u + S.

共11兲

10

v(t) , <S(t)>

5

1

0.5
0

1400

1600

1800

2000

time
1

du

1200

500

1000

1500
2000
time

2500

3000

FIG. 7. 共Color online兲 Comparison of the solution v共t兲 共points兲
with an approximation solution of 具S典⬁ obtained from dv / dt̃ = 0
共line兲. The parameters are the same as those in Fig. 5共e兲. The resultsfor elements 1 共dotted line, red ⫻兲, 2 共solid line, cyan +兲, and 3
共dashed line, blue ⴱ兲 are shown. Notice that each approximation
共line兲 gives a good solution for the behavior 共points兲.

FIG. 8. 共Color online兲 Comparison of the solution u共t兲 共points兲
with the following solution obtained from an adiabatic approximation: u⬘ = S / 共v fix + S兲 from du / dt̃ = 0 共line兲. The parameters are thesame as those in Fig. 5共d兲. The results for elements 1 共dotted line,
red ⫻兲, 2 共solid line, cyan +兲, and 3 共dashed line, blue ⴱ兲 are shown.
Notice that each approximation 共line兲 gives a good solution for the
behavior 共points兲.

Of course, the above solution should be valid only in thelimit ␶s / ␶a → 0, but as shown in Fig. 7, the solution agrees
well with the present case in which ␶s / ␶a ⬃ 共0.5− 0.005兲. In
this rapid-burst-type oscillation, the system enters a state inwhich the average input from all the other elements is cancelled; consequently, v constantly remains in a state that isaway from the original perfect adaptation state.
The rapid-burst-type oscillation occurs when successive
inputs throughout the loop of three elements change sign, or,in other words, there is “frustration” in the loop interactions
共frustration is a term in spin glass theory 关25兴兲. As a simple
example, consider the case in which C21 = 1, C32 = 1, andC13 = −1. When element 1 is active, element 2 is activatedthrough positive interaction with C21 = 1 and element 3 isalso activated through positive interaction with C32 = 1; however, negative interaction with C13 = −1 inactivates element 1.
Hence, the input inconsistent with its state exists. This argument is true even if we assume that element 1 is inactive. For
any value of ui, such opposite inputs remain, and hence,“frustration” among elements exists, as studied in spin statistical mechanics, e.g., spin-glass theory. Unlike the case of
perfect adaptation, this inconsistency in two-to-one antiphaserelationship cannot be eliminated, and hence, some outputs
remain and destabilize the original adapted state. The v variables that cannot change as fast as this opposite input stayalmost constant at v fix = 具S典⬁, whereas only u variables show
fast oscillations.
Next, we study the slow oscillation in which both u and v
components show a periodic change with a time scale much
longer than ␶a 共the original adaptation time scale兲. As vchanges slowly, the adiabatic approximation is valid. In particular, this approximation holds under the limit ␶s / ␶a
˜ → 0. Slow v variables gradually relax to S values with⬃1/␮␶a, while fast u variables relax to the equilibrium values,
defined by du / dt̃ = 0, according to the values of the v variables at each instant within ␶s 共Fig. 8兲. Thus, ui’s approach
the original adapted value and remain close to this value over
the time scale ␶a. For a majority of the time, ui’s are adapted
to the input S from other elements and show only intermit-

026203-7

PHYSICAL REVIEW E 81, 026203 共2010兲

MASAYO INOUE AND KUNIHIKO KANEKO

10
v2 initial

Tsys

1000

100

10

1

10
κ

1

100

0.1
0.1

FIG. 9. 共Color online兲 Dependence of Tsys in the slow cycle on
the strength of the interaction function 关␬ in Eq. 共6兲兴. The symbols⫻ and + represent the results corresponding to networks in Figs.
5共c兲 and 5共d兲, respectively. As the interaction becomes stronger
共␬ becomes larger兲, the characteristic time scale gets longer. The
time scale of the response 共␶s兲 is 0.05–5, and the time scale ofadaptation 共␶a兲 is equal to 10.

101
v1 initial

10

1v2 initial10

v3 initial

10

1

0.1
0.1
10
v3 initialtent, periodic responses, while vi’s change gradually. By
adiabatically eliminating the variable ui, the motion of slow
v variables is obtained. In Fig. 8, we have plotted the present
orbit and the solution obtained by this adiabatic elimination;
we find that the agreement is rather good. Note that in thepresent case, the input S varies between 0.1 and 10, andindeed, when S ⬃ 0.1, the effective response time ␶s= 1 / 共2S兲 ⬃ ␶a / 2. Despite this change in the time scale, theadiabatic approximation agrees well with our simulation
results.
In the slow cycle, the variables u and v undergo a repetitive process in which they approach, remain nearby, and depart from the original fixed point 共u = uⴱ = 0.5, v = vⴱ = 1兲; theduration of the slow cycle is much longer than ␶a. This longtime scale is a result of interaction. To verify this claim, westudied the dependence of the period on several parametersin the interaction function 共5兲. We found that the period depends on the strength of the interaction 关␬␴ / 2 in Eq. 共8兲 or共12兲兴. Since ␶s changes with the value of ␴, we study thedependency of the period on the value of ␬. Up to somevalue of ␬, the interaction function h does not have a sufficiently steep threshold, and hence, the original fixed point isstable and the slow-cycle state does not exist. Beyond somecritical value of ␬, the slow cycle appears as a result of Hopf
bifurcation. With a further increase in ␬ or the interaction
strength ␬␴ / 2, the period increases even though the adaptation time ␶a remains fixed 共Fig. 9兲. The period of burst-type
oscillation, in contrast, remains at ␶s despite the change in ␬.
For most networks that we have numerically studied,
rapid-burst-type oscillation and the slow cycle coexist. Thereexist only two attractors, depending on the initial conditions
共Fig. 10兲. How the two solutions are separated in phase space
˜ → 0. Underis clear when we consider the limit ␶s / ␶a ⬃ 1 / ␮such a limit, the slow cycle is obtained by taking an adiabatic
approximation given by du / dt̃ = 0, according to the value ofv at each instant. On the other hand, rapid-burst-type oscillation is obtained in the condition dv / dt̃ = 0. The two solutions coexist without canceling each other, since the rapidburst-type oscillation traces around the manifold spanned byui共i = 1 , 2 , 3兲 whereas the slow cycle traces around the manifold spanned by vi; the two oscillations occur in spaces or-1
v1 initial

1

0.1
0.1

FIG. 10. 共Color online兲 Basin structure of rapid oscillation 关red
共dark gray兲 ⫻兴 and slow cycle 关blue 共light gray兲 +兴 for different
initial perturbations. The parameters are the same as those in Figs.
5共d兲 and 5共e兲. The perturbation is introduced only in each of theindicated variables, and the remaining variables are set to the steadystate value 共uⴱ = 0.5, vⴱ = 1兲. In order to classify the attractors, Tsysis computed from the return time after a sufficiently long period that
includes transients has been disregarded.

thogonal to each other. It is remarkable that two attractors
whose periods differ by more than two digits coexist in thesame network.
Now, we classify the networks into classes of networks inwhich one of the following occur: 共i兲 perfect adaptation, 共ii兲
simple oscillation, and 共iii兲 coexistence of rapid-burst-type
oscillation and a slow cycle; the classification is based on the
network structure. For cases 共ii兲 and 共iii兲, we have computed
the period of the system 共Tsys兲. The time periods are plottedin Fig. 11 as a function of a network index 共abscissa兲. Since
perfect adaptation corresponds to a fixed point solution, net-

026203-8

DYNAMICS OF COUPLED ADAPTIVE ELEMENTS: …

PHYSICAL REVIEW E 81, 026203 共2010兲

β

1000

1

100
Tsys

3

( κσ2 )

10

1

-1

2

( κσ2 )

2γ γ

0.1

(ii)’

(iii)
0

10

20

β=− 3 3

(ii)

30
40
50
network number

60

70

-1

FIG. 11. 共Color online兲 Periods Tsys of three-element systems
共ordinate兲 plotted against the network number 共abscissa兲. The network numbers are determined arbitrarily but are sorted so that
No.1–25 correspond to case 共iii兲, where rapid oscillation and a slow
cycle coexist, and No. 26–71 correspond to the simple oscillation.
The networks that show perfect adaptation are not included here.
Among the networks that show simple oscillation, No. 26–31 are
networks corresponding to the shaded area in Fig. 12. The results
for five different initial conditions for each network are overlaid.
The values ␶s = 0.05– 5 and ␶a = 10 are chosen here. The periods arecomputed from the return time after a sufficiently long period that
includes transients has been disregarded.

works of type 共i兲 are not included. For the parameter values
chosen here, the time scales for a single element are ␶s
˜ = 10 for slow
= 1 / 共2S兲 = 0.05– 5 for fast response and ␶a = ␮adaptation. The networks with indices less than 25 belong tocategory 共iii兲, where two types of limit cycles coexist. The
rapid oscillation has a period 0.05ⱕ Tsys ⬍ 5, which implies
that the period is within the range of ␶s. The slow cycle hasa period much longer than ␶a: Tsys ⲏ 10␶a for the present
parameter values. The coexistence of the two solutions on
distinctly separated time scales is clearly observable in the
figure. The simple oscillation has a period that is slightly
larger than, but on the same order as, the adaptation timescale ␶a, i.e., ␶a ⱗ Tsys ⱗ 10␶a. For these networks, the simple
oscillation is the unique attractor; as seen in the figure, theperiod of oscillation always lies between the time periods ofrapid and slow oscillations. The classification of the type ofoscillations on the basis of the order of Tsys, ␶s, and ␶a is not
influenced by the time scale for each element.
To classify the cases 共i兲, 共ii兲, and 共iii兲, we again performed
a linear stability analysis of Eq. 共3兲 around the original, perfectly adapted fixed point solution. The six eigenvalues ␭ aresolutions of the following characteristic equation:

冉 冊
␭+

γ

0

1

冉 冊

1 3
1
共␭ + 2兲3 − ␥␭2 ␭ +共␭ + 2兲 − ␤␭3 = 0,
˜
˜
␮
␮
共12兲

where

␤=

冉 冊冉

冊

␬␴ 3 C12C23C31 C13C32C21+,
2
N 1N 2N 3
N 1N 3N 2

共13兲

FIG. 12. 共Color online兲 Phase diagram of the three behaviors ina three-element system characterized by two indices ␤ 共ordinateaxis兲 and ␥ 共abscissa axis兲. ␤ indicates the consistency around thethree elements and ␥ shows the average cooperativity between theneighboring two elements. The definitions of ␤ and ␥ are given in
the text. The green box represents a network that exhibits perfectadaptation behavior 共i兲; the blue 䊊s represent simple oscillation 共ii兲;and the red ⫻s represent the coexistence of a slow cycle and a rapid
oscillation 共iii兲.

␥=

冉 冊冉

冊

␬␴ 2 C12C21 C23C32 C31C13++.
2
N 1N 2
N 2N 3
N 3N 1

共14兲

Note that the two indices ␤ and ␥ characterize a dominant
network structure; ␤ indicates whether the system has a positive or negative loop over the three elements. When ␤ ⬎ 0,each element receives an input traversed the entire loop without changing its sign. Inputs to support and suppress activity
are received through the loop when an element is active andinactive, respectively. Hence, three elements can behave inunison or in two-to-one antiphase synchronization. If ␤ = 0,
the input through the loop is cancelled, while for ␤ ⬍ 0, as
previously mentioned, the frustration already exists, andthrough the loop, each element receives an input inconsistent
with its state. The index ␥, on the other hand, represents theaverage cooperativity between the two neighboring elements.
When ␥ ⬎ 0, two elements can behave in a synchronizedmanner, while when ␥ ⬍ 0, they cannot do so because offrustration between the two elements.
We study the sign of the real part of the eigenvalues from
Eq. 共12兲 that determine the stability of the steady state about
each of 78 possible networks and examine the dependence ofthe sign on ␤ and ␥. First, if ␤ = 0 and ␥ ⱕ 0, there are six
negative eigenvalues, where the perfect adaptation state isstable. If ␤ ⬎ 0 or ␤ = 0 and ␥ ⬎ 0, there are two positive real
eigenvalues and four eigenvalues with negative real parts. If␤ ⬍ 0, there are two negative real eigenvalues and four eigenvalues with positive real parts. At the lattice region in
Fig. 12 共where either ␥ ⬍ 0 , ␤ ⬍ 0 or 0 ⱕ ␥ , ␤ ⬍ − 23␥ 冑 ␥3 兲,
there are two pairs of complex conjugate eigenvalues withpositive real parts. On the other hand, at the shaded region in
Fig. 12 共where 0 ⱕ ␥ , − 23␥ 冑 ␥3 ⱕ ␤ ⬍ 0兲, there are two realpositive eigenvalues and a pair of complex conjugate eigenvalues with positive real parts 关26兴. The limiting case of˜ → 0 leads to a double zero eigenvalues. In this sense, the1/␮

026203-9

PHYSICAL REVIEW E 81, 026203 共2010兲

MASAYO INOUE AND KUNIHIKO KANEKO

bifurcation is similar with Takens-Bogdanov 共TB兲 bifurcation, even though there is no saddle-node bifurcation in ourmodel 共as there exists only one fixed point兲. Recall that unfolding of TB bifurcation typically leads to a homoclinic
orbit, which often introduces a burst-type oscillation and aslow cycle 关27兴. Even though direct application of TB bifurcation analysis is not possible here, appearance of the bursttype oscillation and a slow cycle may be understood througha similar analysis.
Interestingly, this classification on the basis of the eigenvalues of the stability matrix corresponds to the classificationon the basis of the three behaviors 共i兲–共iii兲. We classify the
networks that show a particular behavior on the basis of theindices ␤ and ␥ and plot each type of behavior in a phase
diagram in Fig. 12. As shown by the linear stability analysis,
perfect adaptation occurs when ␤ = 0 and ␥ ⱕ 0. Simple oscillation appears under the condition ␤ ⬎ 0 or ␤ = 0 and ␥⬎ 0 and also in the shaded region corresponding to− 23␥ 冑 ␥3 ⱕ ␤ ⬍ 0 and ␥ ⱖ 0 关28兴. The coexistence of rapidburst-type oscillation and the slow cycle occurs at ␤ ⬍ 0,
except under the conditions corresponding to the shadedarea.
In other words, phase 共ii兲 appears when the system has nofrustration 共␤ ⬎ 0兲. Because the elements can cooperate evenwhen the original fixed point is unstable, the system belongs
to the same category as that for ␣ ⬎ 0 in the two-elementcase, and thus, a limit cycle is generated. Phase 共iii兲 appears
when the system has frustration. There are two pairs of complex conjugate eigenvalues. One pair has a large positive realpart and the corresponding eigenvector is orthogonal to theadiabatic space. The real part of the other pair is positive but
has a much smaller value; the eigenvector of this pair is notorthogonal to the adiabatic space, and hence, in the adiabaticcase, the instability of the fixed point solution is weak, suggesting the emergence of a slow oscillatory mode. Therefore,
the coexistence of the two types of limit cycle attractors, onewith a short period and the other with a much longer period,
is a natural outcome of the appearance of phase 共iii兲.
In addition to a pair of complex conjugate eigenvalues
with a positive real part, the shaded area in Fig. 12 corresponds to two more positive eigenvalues that are real. Unlike
the region for phase 共iii兲, there is no additional complex conjugate pair. In fact, only a single limit cycle arising from thesingle pair of complex eigenvalues with a real part exists,and consequently, a simple oscillation is developed, as inphase 共ii兲. The modes with two positive real eigenvalues do
not generate a different attractor. Even though three-body
frustration exists, as indicated by ␤ ⬍ 0, the cooperativity of
two-body interactions cancels the frustration. However, there
still seems to be a difference between the behavior when ␤⬍ 0 and the simple oscillation when ␤ ⬎ 0. The period in thisregion is shown in Fig. 11 as the period of network 共ii兲⬘.
V. CASES WITH A LARGER NUMBER
OF ELEMENTS

Now, we study a case involving a larger number of elements 共N ⱖ 4兲. In such cases, we once again observed three
types of behaviors; 共i兲 perfect adaptation, 共ii兲 simple oscilla-tion, and 共iii兲 coexistence of a rapid-burst-type oscillation
and a slow cycle. In addition to these basic behaviors, more
complicated behaviors, e.g., quasiperiodic or aperiodic oscillation, are also observed, as N is further increased 共Fig. 13兲.
To classify the network behaviors, we again performed
linear stability analysis of Eq. 共3兲 around the original, perfectly adapted fixed point solution. In this case, the dominant
structures in the interaction network are again characterized
by the coefficients of the characteristic equation, which indicate the n-loop frustration in the network. The coefficientsare given by loop structures over n elements
共n = 2 , 3 , ¯ , N兲 containing a product of smaller loops as follows:

冉 冊兺
冉 冊兺

␬␴ 2 CijC ji,
2 具ij典 NiN j

共15兲

␬␴ 3
CijC jkCki,
2 具ijk典 NiN jNk

共16兲

L共2兲 =

L共3兲 =

L共4兲 =

冉 冊 冉兺冊

␬␴ 4
CijC jkCklCli
CijC ji CklClk
− 兺
⫻
,
2
N
N
N
N
N kN li j k l具ijkl典具ij典具kl典 NiN j
共17兲

where the summation by 具ij典 runs over all possible pairs withi ⫽ j and that by 具ijk典 runs over all triplets with i ⫽ j, j ⫽ k
and k ⫽ i; the pairs and triplets are chosen from among the N
elements.
For the case in which more than four elements are involved, the n-loop frustration is characterized by the L共n兲that is determined by a product of L共m j兲共m j ⱖ 2 , 兺 M
j=1m j= n兲, as follows:
L共n兲 =

C i1i2C i2i3 ¯ C ini1

兺 Ni Ni ¯ Ni
具i ¯i 典
1

n

1

+兺 ¯兺
m1

mM

2

共− 1兲

n

¯ 共− 1兲mM +1
L共m1兲 ¯ L共m M 兲.
共− 1兲n+1

m1+1

共18兲
For N = 4, we again draw a phase diagram for the fourtypes of behaviors, which are plotted with respect to thethree L共n兲 indices, by classifying the networks on the basisof the period of the attractor共s兲. The classification of the
networks on the basis of the sign of L共n兲 seems to be valid
as in the case with N = 3 共Fig. 14兲. Here networks for perfectadaptation always satisfy L共2兲 ⱕ 0 , L共3兲 = 0 , L共4兲 ⱕ 0.
Conditions for L共2兲 and L共4兲 are further specified as either
L共4兲 ⬇ 0 and L共2兲 ⱕ 0, or L共4兲 = L共2兲 + 1 with L共4兲 ⱕ 0.
These conditions for perfect adaptation are explained in the
Appendix.
We have also studied a case with N = 10 elements. The
three basic behaviors are often observed, whereas complicated behaviors are also observed, including aperiodic oscillations 共Fig. 15兲. As N becomes larger, there is a larger variety in introducing frustrations among elements due todifferent signs in the interaction coefficients Cij. Even withregards to the linear stability analysis, the combinatorial va-

026203-10

DYNAMICS OF COUPLED ADAPTIVE ELEMENTS: …
(a)1

PHYSICAL REVIEW E 81, 026203 共2010兲

10

1

1

0.5

0.6

v(t)

u(t)

0.8

0.4
0.2
0
0

100

200

300
time

400

500

0.1
0

(b) 1

100

200

300
time

400

L(3) 0

500

10

-0.5

0.6

v(t)

u(t)

0.8

0.4

1

-1
-2

0.2
0
0

500

time

1000

1500

0.1
0

1

(c)

500

time

1000

-1

0

1

2

1

2

L(2)

1500

10

1

1

0.5

0.6

v(t)

u(t)

0.8

0.4
0.2
0
0

(d)

500

time

1000

1500

0.1
0

1

500

time

1000

L(4) 0

1500

10

-0.5

0.6

v(t)

u(t)

0.8

0.4

1

-1
-2

0.2
0
500

(e)

500.2

500.4 500.6
time

500.8

501

0.1
500

1

500.2

500.4 500.6
time

500.8

501

v(t)

u(t)

0.8

0.4

1

0.2
0
0

1000

time

2000

3000

0.1

0

0

L(2)

10

0.6

-1

1000

time

2000

3000

FIG. 13. 共Color兲 Time series of u 共left兲 and v 共right兲 for the
four-element system. 共a兲 共C12 = 1 , C13 = −1, C14 = 1 , C21 = 0, C23= −1 , C24 = 1 , C31 = −1 , C32 = 0 , C34 = −1 , C41 = 0 , C42 = 0, C43= −1兲 corresponds to oscillatory behavior. 共b兲 共C12 = −1 , C13= 0 , C14 = −1 , C21 = −1 , C23 = 1 , C24 = 0 , C31 = 0 , C32 = −1, C34= 0, C41 = 1 , C42 = −1 , C43 = −1兲 and 共c兲 共C12 = 0 , C13 = 1 , C14= −1 , C21 = 0 , C23 = 1, C24 = 0 , C31 = 0 , C32 = −1 , C34 = −1, C41= 1 , C42 = 1 , C43 = 0兲 correspond to slow-cycle behavior and 共d兲关same network as in 共c兲兴 corresponds to rapid oscillation. 共e兲 共C12
= 1 , C13 = −1 , C14 = −1 , C21 = −1 , C23 = 0 , C24 = 1, C31 = 1,
C32 = −1, C34 = −1, C41 = −1 , C42 = 1 , C43 = 1兲 corresponds to aperiodic oscillatory behavior. Note the difference in the scales of theabscissa for different figures. In each picture, the red line shows the
dynamics of element 1; the cyan line, the dynamics of element 2;
the blue line, the dynamics of element 3; and the green line, thedynamics of element 4.

riety in defining L共n兲 is huge. Although the existence ofperfectly adapted fixed point solution is always assured, its
instability leads to complex dynamics with multiple coexisting attractors. Here, frustration-induced chaos in a coupled
dynamical system with a large degree of freedoms seems to

FIG. 14. 共Color兲 Classification of four-element systems on thebasis of L共2兲 共abscissa; upper figure兲 and L共3兲 共ordinate; upper
figure兲 as well as L共2兲 共abscissa; lower figure兲 and L共4兲 共ordinate;
lower figure兲. The indices are defined in the text. The green box
represents a network that exhibits perfect adaptation behavior; theblue 䊊s represent simple oscillation; the red +s represent the coexistence of a slow cycle and a rapid oscillation; the black ⴱs represent aperiodic oscillation. 1,000 networks are generated randomly
and each network is labeled as a particular type on the basis of theperiod Tsys.

be rather general, as is also discussed in 关29,30兴.
We have made a phase diagram for N = 10 as in Fig. 14 toclassify the behaviors with regards to the frustration indices,
L共2兲 and L共3兲. The diagram has similarity with the caseN = 4 共data not shown兲. Interestingly, networks with perfectadaptation are always located at the region L共2兲 ⱕ 0 and
L共3兲 = 0. The region for aperiodic oscillations is embedded
within the region for rapid and slow oscillations, and thecondition to distinguish the former region from the othertypes is not obtained just by the lowest two frustration indices.
VI. DISCUSSION AND CONCLUSION

In this study, we have introduced a system of coupledadaptive elements, studied the behavior of the dynamics ofthese elements, and classified the behavior on the basis of its
relationship with the network structure. For a system withtwo elements, two types of behaviors, perfect adaptation and026203-11

PHYSICAL REVIEW E 81, 026203 共2010兲

MASAYO INOUE AND KUNIHIKO KANEKO

(a)1

10

0.8

v(t)

u(t)

0.6
1

0.4
0.2
0

(b)

0

100

200

time

300

400

500

1

0.1

0

100

200

time

300

400

500

10

0.8

v(t)

u(t)

0.6
1

0.4
0.2
0

(c)

0

500

time

1000

1500

1

0.1

0

500

1000

time

1500

10

0.8

v(t)

u(t)

0.6
1

0.4

∆t=10

0.2
0
300

(d)

302

304

time

306

308

310

1

0.1
300

302

304

time

306

308

310

10

0.8

v(t)

u(t)

0.6
1

0.4
0.2
0

(e)

0

500

time

1000

1500

1

0.1

0

500

time

1000

1500

10

0.8

v(t)

u(t)

0.6
1

0.4
0.2
0

0

200

400

time

600

800

1000

0.1

0

200

400

time

600

800

1000

FIG. 15. 共Color online兲 Time series of u 共left兲 and v 共right兲 forthe N = 10 system. Panel 共a兲 corresponds to a simple oscillation. 共b兲
corresponds to a slow cycle and 共c兲 corresponds to a rapid oscillation, where 共b兲 and 共c兲 are examples from a same network. 共d兲 and共e兲 are examples of aperiodic behaviors. Note the difference in thescales of the abscissa for different figures. In each picture, the dynamics of each element is overwritten with a thin line and theaverage of them is shown with a black bold line.

simple oscillation, were observed. For a system with threeelements, rapid-burst-type oscillation and a slow cycle werediscovered in addition to the above two behaviors. The rapid
oscillation and the slow cycle may coexist depending on theinitial conditions. These dynamics are a result of the existence of two distinct time scales for the behavior of the adaptive element, i.e., fast response and slow adaptation processes; hence, the dynamics are characteristic of coupledadaptive systems. The dynamics of each network are classified on the basis of the network structure 共in specific, a combination of positive or negative feedback among elements兲.
Cooperativity among elements leads to simple oscillation,
whereas frustration among elements leads to the coexistenceof a rapid bursting oscillation and a slow cycle. The index forfrustration is defined by the product of signs of the interaction matrix elements over the entire loop. This index alsocharacterizes the eigenvalues in the linear stability matrix
around the steady state.
The period of the rapid-burst-type oscillation is on theorder of the response time in a single element, implying thetotal absence of adaptation; in contrast, the period of theslow cycle is much larger than the original time scale ofadaptation, where frustration among elements caused by interaction was important. It is remarkable that the attractors
with such large differences in their periods generally coexist.
Coupled oscillators have attracted much attention andmany general concepts as synchronization have been developed 关14兴. In oscillator networks, frustration causes synchronization 关31兴 or leads to an ordered state characterized byquasientrainment 关32兴. Synchronization in coupled chaotic
oscillators have also been studied 关33,34兴. However, acoupled system of adaptive elements has not been studied
thus far. Our study here demonstrates that cooperative dynamics generally arise as a result of frustration in coupledadaptive systems.
Differences in the response time scales are commonly observed in biological systems. The frequency of the calcium
oscillations in a cell often depends on the cellular state, and
it is proposed that this frequency is related to cell type differentiation 关35兴. Neuroscience is another example, where
two-state fluctuations in neural activity have recently attracted much attention 关36–38兴; where a neuron spontaneously switches between the up state with rapid-burst firingand down state with rare firing. Some models have been
proposed to explain the mechanics of two-state fluctuations
by taking into account the recursive excitable interaction
关39,40兴; however, the mechanism is not yet fully understood.
Although our model has not been developed specifically forthe neural system, it shows excitatory response that is common in the neural system. Our assertion in this context is
simple: frustration in interactions among adaptive elements
generates bistability between burst firing and rare firing.
The extension of the three-element system to a larger
number of elements is straightforward. Once again, we observed the coexistence of burst-type oscillations and a slow
cycle. The behaviors are classified on the basis of the frustration in a loop structure over n elements 共n = 2 , 3 , . . . , N兲 in
an interaction network; the frustration in the interaction network is characterized by the coefficients in a characteristic
equation in the linear stability analysis around the steadystate. These frustration indices are again relevant to the classification.
The generation of long-term dynamics and differences intime scales is essential in determining cellular and neural
behaviors. Studies focusing on the time-scale difference, a
coupled chaotic system with a variety of time scales 关41兴 orcoupled phase oscillators with diversification of the timescales 关42,43兴 have been carried out. Depending on the network structure, such behavior is inherent to a coupled adaptive system in a network, despite the simplicity of the modelthat is utilized in this case. Here, we have only studied thecase in which the coupling between identical elements is thesame. In biology, elements are often heterogeneous and thenumber of elements can be very large. Such complex cases

026203-12

DYNAMICS OF COUPLED ADAPTIVE ELEMENTS: …

PHYSICAL REVIEW E 81, 026203 共2010兲

are currently being investigated, and these investigations will
be reported elsewhere.

冕

ACKNOWLEDGMENTS=i

The authors would like to thank A. Awazu, S. Ishihara, K.
Fujimoto, and D. Shimaoka for their valuable comments, and
M. Tachikawa for his advice on the analysis of the stability
condition in the Appendix. M.I. was partially supported by
the JSPS.
APPENDIX: Conditions for Perfect Adaptation for N = 4

For N = 4 case, we perform a linear stability analysis of
Eq. 共3兲 around the original, perfectly adapted fixed point
solution under the limit of ␶s / ␶a → 0. The four eigenvalues ␭are solutions of the following characteristic equation:
共␭ + 2兲4 − L共2兲共␭ + 2兲2 − L共3兲共␭ + 2兲 − L共4兲 = 0. 共A1兲
Here, we write L共n兲 = 共 ␬2␴ 兲nL⬘共n兲 and set 2共␭ + 2兲 / 共␬␴兲 = x,
then, the Eq. 共A1兲 becomesx4 − L⬘共2兲x2 − L⬘共3兲x − L⬘共4兲 = 0.

共A2兲Let H共x兲 = x − a2x − a3x − a4 ⬅ h共x兲 − a3x. For a system of theperfect adaptation, the real part of all eigenvalues 共␭兲 cannot
be positive. Since ␭ = −2 + ␬␴2 x and 共␬␴兲 / 2 Ⰷ 2, the condition
for the stability means that the real part of the above solution
cannot be positive. To obtain the condition, we first showthat at least one root of the equation H共x兲 = 0 has a positive
real part as long as a3 ⫽ 0.
First we study the case with a4 ⫽ 0. Consider the integral
on complex plane
4

1dx,
C1+C2 H共x兲

共A3兲

where C1 = ri共−R ⬍ r ⬍ R兲 is a path on the imaginary axis andC2 is a semicircle with radius R with a positive real part ofx共C2 = 兵Rei␪ 兩 −␲ / 2 ⱕ ␪ ⱕ ␲ / 2其兲. First of all, we note H共x兲⫽ 0 on C1 from the assumption that a3 ⫽ 0. Indeed, when x= ri is a pure imaginary number, h共ri兲 = r4 + a2r2 − a4 is real
while −a3ri is pure imaginary. Thus H共x兲 cannot vanish
when x is a pure imaginary number.
Following the residue theorem, Eq. 共A3兲 should be zero if
there were no x satisfying H共x兲 = 0 inside the C1 + C2. Bytaking the limit R → ⬁,

冕

1
dx =
H共x兲
C2

On the other hand,

冕

−␲/2

␲/2

d␪

i Rei␪
1⬀→ 0.
H共Rei␪兲 R3

R

dr

0
R

= 2i

0

dr

1
+
h共ri兲 − a3ri

冕

R

0

dr

1
h共− ri兲 + a3ri冎

r 4 + a 2r 2 − a 4.
共r + a2r2 − a4兲2 + 共a3r兲2
4

At least for a certain set of parameter values of a2 and a4,
this integration does not vanish. For this set of parameters,
Eq. 共A3兲 cannot be zero, which means that there exists atleast one root of H共x兲 = 0 with positive real part.
Now we consider the parameter change in a2 and a4. If
this root of H共x兲 can be removed from the positive real region, the solution should get across the imaginary axis andmove to the negative real region 共because the zero-pointshould continue to exist as it is a root of a polynomial equation兲. This means at certain parameter values of a2 and a4,
the root of H共x兲 should be pure imaginary. However, as already mentioned, H共x兲 cannot be zero if x is pure imaginary
as long as a3 ⫽ 0. Therefore, as long as a3 ⫽ 0 and a4 ⫽ 0,
there must be a root of H共x兲 with a positive real part.
For the case with a4 = 0,
H共x兲 = x4 − a2x2 − a3x = x共x3 − a2x − a3兲 = xG共x兲.

2

冕

冕
再冕
冕

R
1
1
dx = idr
h共ri兲
− a3ri
H共x兲
−R
C1

In this case again, G共x兲 ⫽ 0 when x is a pure imaginary number as long as a3 ⫽ 0, and we can repeat the above argument
to show again that for a3 ⫽ 0, there is a root of H共x兲 = 0 witha positive real part.
Now, we have shown that there exists at least one root of
H共x兲 = 0 with a positive real part as long as a3 ⫽ 0. Hence, forall eigenvalues 共␭兲 to have zero or negative real part, a3 = 0 isnecessary when 共␬␴兲 / 2 Ⰷ 2. Next, we study roots of H共x兲under a3 = 0.
In this case, when x = −xr ⫾ ixi共xr , xi ⬎ 0兲 are solutions of
H共x兲 = h共x兲 = 0, x = xr ⫾ ixi are also its solutions. This means
H共x兲 = 0 can have only pure imaginary or 0 solutions, not tohave a root with a positive real part. From H共x兲 = 0, we getx2 =a2 ⫾ 冑a22 + 4a4
2

共A4兲

and x2 ⱕ 0 is necessary according to the above discussion. It
is straightforward to show that there are only two possibilities. 共i兲 a2 ⱕ 0 , a4 = 0 where x = 0 , ⫾ i冑兩a2兩 and 共ii兲 a2 ⱕ
−1 , a4 = a2 + 1 where x = ⫾ i , ⫾ i冑兩a2 + 1兩. Hence, for perfect adaptation, L共4兲 = 0 and L共2兲 ⱕ 0, or L共4兲 = L共2兲 + 1 with
L共4兲 ⱕ 0, in addition to L共3兲 = 0. These results are consistent
with Fig. 14. A few examples that might violate the above
condition with a slightly negative L共4兲 are due to finitenessof ␬␴, as the real part of ␭ = −2 + ␬␴2 x can be negative even if
the real part of x is slightly positive.

026203-13

PHYSICAL REVIEW E 81, 026203 共2010兲

MASAYO INOUE AND KUNIHIKO KANEKO
关1兴 D. E. Koshland, Jr., Science 196, 1055 共1977兲.
关2兴 D. E. Koshland, Jr., A. Goldbeter, and J. B. Stock, Science
217, 220 共1982兲.
关3兴 F. Oosawa and Y. Nakaoka, J. Theor. Biol. 66, 747 共1977兲.
关4兴 Y. Nakaoka, H. Tokui, Y. Gion, S. Inoue, and F. Oosawa, Proc.
Jpn. Acad., Ser. B: Phys. Biol. Sci. 58, 213 共1982兲.
关5兴 B. E. Knox, P. N. Devreotes, A. Goldbeter, and L. A. Segel,Proc. Natl. Acad. Sci. U.S.A. 83, 2345 共1986兲.
关6兴 S. Asakura and H. Honda, J. Mol. Biol. 176, 349 共1984兲.
关7兴 N. Barkai, U. Alon, and S. Leibler, C. R. Acad. Sci., Ser. IV
Phys. Astrophys. 2, 871 共2001兲.
关8兴 P. Francois and E. D. Siggia, Phys. Biol. 5, 026009 共2008兲.
关9兴 W. Ma, A. Trusina, H. El-Samad, W. A. Lim, and C. Tang, Cell
138, 760 共2009兲.
关10兴 S. M. Block, J. E. Segall, and H. C. Berg, J. Bacteriol. 154,
312 共1983兲.
关11兴 E. Karatan, M. M. Saulmon, M. W. Bunn, and G. W. Ordal, J.
Biol. Chem. 276, 43618 共2001兲.
关12兴 A. P. Gasch, P. T. Spellman, C. M. Kao, O. Carmel-Harel, M.
B. Eisen, G. Storz, D. Botstein, and P. O. Brown, Mol. Biol.
Cell 11, 4241 共2000兲.
关13兴 S. Stern, T. Dror, E. Stolovicki, N. Brenner, and E. Braun,
Mol. Syst. Biol. 3, 106 共2007兲.
关14兴 Y. Kuramoto, Chemical Oscillation, Waves, and Turbulence
共Dover, New York, 2003兲.
关15兴 K. Kaneko and I. Tsuda, Complex Systems: Chaos and Beyond
共Springer-Verlag, Berlin, Heidelberg, 2001兲.
关16兴 K. Kaneko, Physica D 41, 137 共1990兲.
关17兴 M. Inoue and K. Kaneko, Phys. Rev. E 74, 011903 共2006兲.
关18兴 R. Erban and H. G. Othmer, SIAM J. Appl. Math. 65, 361
共2004兲.
关19兴 A. Levchenko and P. A. Iglesias, Biophys. J. 82, 50 共2002兲.
关20兴 S. Sawai, P. A. Thomason, and E. C. Cox, Nature 共London兲
433, 323 共2005兲.
关21兴 T. Gregor, K. Fujimoto, and S. Sawai 共unpublished兲.
关22兴 B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P.
Walter, Molecular Biology of the Cell, 5th ed. 共Garland Science, New York, 2008兲.
关23兴 I. Stewart, M. Golubitsky, and M. Pivato, SIAM J. Appl. Dyn.
Syst. 2, 609 共2003兲.
关24兴 M. Golubitsky and I. Stewart, The Symmetry Perspective:
From Equilibrium to Chaos in Phase Space and PhysicalSpace 共Birkhäuser, Basel, 2003兲.
关25兴 H. Nishimori, Statistical Physics of Spin Glasses and Information Processing: An Introduction 共Oxford University Press,
New York, 2001兲.
关26兴 Taking the adiabatic approximation, Eq. 共12兲 goes to ␭3兵共␭
+ 2兲3 − ␥共␭ + 2兲 − ␤其 = 0, i.e., the three eigenvalues are 0 and theother three are roots of 共␭ + 2兲3 − ␥共␭ + 2兲 − ␤ = 0, which are2␥ ␥classified as follows. 共a兲 ␥ ⬎ 0 , ␤ ⱕ 8 − 2␥ and ␤ ⱕ 3 冑 3 and2␥ 冑 ␥
− 3 3 ⱕ ␤, three negative real eigenvalues. 共b兲 ␤ ⱕ 8 − 2␥ and4␥ − 64ⱕ ␤ and excluding 共a兲 area, one negative real eigenvalue and a pair of complex conjugate eigenvalues with nega2␥ ␥
tive real part. 共c兲 8 − 2␥ ⬍ ␤ and 3 冑 3 ⬍ ␤, one positive real
eigenvalue and a pair of complex conjugate eigenvalues with2␥ ␥
negative real part. 共d兲 8 − 2␥ ⬍ ␤ and ␤ ⱕ 3 冑 3 , one positive2␥ ␥and two negative real eigenvalues. 共e兲 ␤ ⱕ 8 − 2␥ and − 3 冑 3ⱕ ␤, one negative and two positive real eigenvalues. 共f兲 ␤2␥ ␥⬍ 4␥ − 64 and ␤ ⬍ − 3 冑 3 , one negative real eigenvalue and apair of complex conjugate eigenvalues with positive real part.
Therefore, the boundary between shaded and dotted area in2␥ ␥Fig. 12 is given by ␤ = − 3 冑 3 .
关27兴 G. Dangelmayr and E. Knobloch, Philos. Trans. R. Soc. London, Ser. A 322, 243 共1987兲.
关28兴 With the present normalization by Ni = 兺Nj=1兩Cij兩, all networks2␥ ␥with ␥ ⱖ 0 , − 3 冑 3 ⱕ ␤ ⬍ 0 locate just on the boundary of theshaded area. We used another type of normalization by Ni= 1 / 2 and studied the behaviors, where networks inside theshadedarea
共e.g.,
C12 = −1 , C13 = −1 , C21 = −1 , C23= −1 , C31 = −1 , C32 = 0 with ␤ = −1 / 8 , ␥ = 1 / 2兲 showed a
simple oscillation.
关29兴 H. Bersini and V. Calenbuhr, J. Theor. Biol. 188, 187 共1997兲.
关30兴 M. H. Jensen, S. Krishna, and S. Pigolotti, Phys. Rev. Lett.
103, 118101 共2009兲.
关31兴 D. H. Zanette, Europhys. Lett. 72, 190 共2005兲.
关32兴 H. Daido, Phys. Rev. Lett. 68, 1073 共1992兲.
关33兴 M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, Phys. Rev.
Lett. 76, 1804 共1996兲.
关34兴 M. G. Rosenblum, A. S. Pikovsky, and J. Kurths, Phys. Rev.
Lett. 78, 4193 共1997兲.
关35兴 R. E. Dolmetsch, K. Xu, and R. S. Lewis, Nature 共London兲
392, 933 共1998兲.
关36兴 C. J. Wilson and Y. Kawaguchi, J. Neurosci. 16, 2397 共1996兲.
关37兴 E. A. Stern, A. E. Kincaid, and C. J. Wilson, J. Neurophysiol.
77, 1697 共1997兲.
关38兴 J. Anderson, I. Lampl, I. Reichova, M. Carandini, and D. Ferster, Nat. Neurosci. 3, 617 共2000兲.
关39兴 A. Compte, M. V. Sanchez-Vives, D. A. MaCormick, and X. J.
Wang, J. Neurophysiol. 89, 2707 共2003兲.
关40兴 S. Kang, K. Kitano, and T. Fukai, Neural Networks 17, 307
共2004兲.
关41兴 K. Fujimoto and K. Kaneko, Chaos 13, 1041 共2003兲.
关42兴 D. Armbruster and P. Chossat, Phys. Lett. A 254, 269 共1999兲.
关43兴 M. Tachikawa and K. Fujimoto, Europhys. Lett. 78, 20004
共2007兲.

026203-14