Natural Frequency Adaptation as a Mechanism for Phase Synchronization in Neural Systems

Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011

Natural Frequency Adaptation as a Mechanism for Phase Synchronization in Neural Systems Ivan Y. Tyukin ∗ Viktor B. Kazantsev ∗∗ ∗

University of Leicester, Leicester, LE1 7RH, UK (e-mail: I.Tyukin@ le.ac.uk). ∗∗ Institute of Applied Physics, Ul’yanova 46 Str., Nizhny Novgorod, 603950 Russia (e-mail: [email protected]).

Abstract: We consider the problem of adaptive phase synchronization in a system of coupled phase oscillators. The system is limited to a pair of phase oscillators of which the natural frequencies are uncertain, yet they can be tuned. Coupling between the oscillators is allowed to depend on the natural frequencies nonlinearly. Three different mechanisms for the natural frequency self-tuning leading to the emergence of synchronous states are presented and analyzed. The problem is relevant for explaining and modelling phase synchronization and purposeful frequency adaptation in neuronal networks. Keywords: phase synchronization, neural systems, weak attractors, adaptive control 1. INTRODUCTION The problem of synchronization received substantial attention in the literature in the past. This is hardly surprising for the phenomenon of synchronization is widely observed in nature Strogatz (2003), and hence is a subject of intense research in physics, biology, mathematics, and engineering. The notion of synchronization itself is multifaceted since many types and definitions of synchronization are now available (see Blekhman et al. (1997) for review). Here we concentrate mostly on what is called phase synchronization Pikovsky et al. (2001): the case in which phases of interacting compartments carry information about dynamical state of the system but amplitudes of the signals are not relevant. The problem we study here is motivated by empirical observations of electrical activity in neuronal cultures, also known as animats. When left unattended the cultures grow and give rise to spontaneous electrical activity Pasquale et al. (2008), Kazantsev et al. (2010). After a number of days of development external stimulation can be used to induce in-phase synchrony in the network over a sustainable period of time. We hypothesize that the synch is a result of internal adaptation in the cultures to stimulation. The question, however, is that these adaptation mechanisms leading to synchronization are unknown. In this brief article we would like to make a step forward in answering to this question Instead of concentrating on detailed modelling of cells in neuronal cultures we consider general models of phase oscillators in which the coupling is specified by some functions of individual phases and, perhaps, natural frequencies. Since no specific models are considered we resort to rather general description of non-identical phase oscillators of which the natural frequencies are allowed to be different. Specifically, we wish 978-3-902661-93-7/11/$20.00 © 2011 IFAC

to be able to answer to the following set of questions: if there are any adaptation mechanisms exist making phase synchronization in the system possible, and how these mechanisms can be described and analyzed. The paper is organized as follows. Section 2 contains mathematical statement of the problem, in Section 3 we present main results of the paper, the results are discussed in Section 4, and Section 5 concludes the paper. 2. PROBLEM FORMULATION Consider the following systems of coupled phase oscillators ( φ˙ 1 = ω1 + f1 (φ1 , φ2 , t) Σ1 : (1) φ˙ 2 = ω2 + f2 (φ1 , φ2 , t) ( φ˙ 1 = ω1 + g1 (φ1 , φ2 , ω1 , t) Σ2 : (2) φ˙ 2 = ω2 + g2 (φ1 , φ2 , ω2 , t) where φ1 , φ2 stand for the phases of oscillators, ω1 , ω2 ∈ R are the frequencies, and f1 , f2 , g1 , g2 are some functions modelling coupling between phases of oscillations. Coupling in system Σ1 is assumed to be independent on the natural frequencies ω1 , ω2 , whereas coupling between phases of oscillations in system Σ2 is supposed to depend explicitly on ω1 , ω2 . Let us suppose, in addition, that both systems, Σ1 and Σ2 , are forward-complete, i.e. that following assumption holds: Assumption 1. Let ω1 (t), ω2 (t) ∈ C 0 [0, T ], then solutions φ1 (t, φ0,1 ), φ2 (t, φ0,2 ), φ0,1 , φ0,2 ∈ R are defined for every t ∈ [0, T ], T > 0. If the values of natural frequencies, ω1 , ω2 , are known then asymptotic properties of φ1 and φ2 can in principle be derived from equations (1), (2). This can be followed by applying a compensatory feedback ensuring that phase

11755

10.3182/20110828-6-IT-1002.02101

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

synchrony occurs. The problem, however, is that the values of ω1 , ω2 may not be known exactly. Thus in order to ensure successful synchronization, adaptation or finetuning of natural frequencies of oscillators may be needed. Here we consider and analyze the case in which natural frequency of one of the oscillators, say ω1 , is allowed to vary, and the natural frequency of the other, ω2 , is supposed to be fixed and unknown. In particular, we are looking for a function ω1 (φ1 , φ2 , t), and a set Ω such that if ω1 is replaced with ω1 (φ1 , φ2 , t) in (1), (2) then the following should hold lim φ1 (t, φ0,1 ) − φ2 (t, φ0,2 ) = φe , ∀φ0,1 , φ0,2 ∈ Ω, t→∞

where φe is a (possibly unknown) phase difference value. For the sake of simplicity, in what follows we suppose that the functions f1 , f2 , g1 , g2 already incorporate any compensatory feedbacks required for synchronization, and that the systems would synchronize provided that the value of ω1 would be chosen properly. We present three different mechanisms ω1 (φ1 , φ2 , t) leading to the emergence phase synchrony the systems. Two mechanisms are presented for system Σ1 , and one mechanism is provided for Σ2 . These mechanisms vary in their requirements for the functions f1 , f2 , g1 , g2 , and also they lead to systems the dynamics of which is qualitatively different. 3. RESULTS The following notational agreements are used throughout the paper. Symbol C k denotes the space of functions that are at least k times differentiable; K denotes the class of all strictly increasing functions κ : R≥0 → R≥0 such that κ(0) = 0. If, in addition, lims→∞ κ(s) = ∞ we say that κ ∈ K∞ . Further, Ke (or Ke,∞ ) denotes the class of functions of which the restriction to the interval [0, ∞) belongs to K (or K∞ ). Symbol KL denotes the class of functions β : R≥0 × R≥0 → R≥0 such that β(·, s) ∈ K and β(r, ·) is monotonically decreasing for each s, r ∈ R≥0 . Finally, ϕ denotes the phase difference between oscillators ϕ(t) = φ1 (t, φ0,1 ) − φ2 (t, φ0,2 ) We will start with considering the case in which the phase oscillators are given by (1).

are defined for all t ≥ 0, and moreover limt→∞ ϕ(t) = φe . Proof. Adaptation scheme (3) is essentially the velocity gradient algorithm Miroshnik et al. (1999), and we just rephrase the argument here for consistency. Since the right-hand side of the combined system is continuous, there exists an interval [0, T ], T > 0 such that solutions of the combined system are defined on [0, T ]. Consider the following function 1 V1 (ϕ, φe , ω1 , ω2 ) = V (ϕ, φe ) + (ω1 − ω1∗ (ω2 ))2 2γ The function V1 is positive definite wrt ϕ−φe , ω1 −ω1∗ (ω2 ), and its derivative wrt t is non-positive: V (ϕ, φe ) V (ϕ, φe ) V˙ 1 = ϕ˙ − (ω1 − ω1∗ (ω2 )) ∂ϕ ∂ϕ (4) ∂V (ϕ, φe ) = ϕ˙ ≤ −ρ(|ϕ − φe |) ∂ϕ ∗ ω1 =ω1 (ω2 )

Hence (ω1 (t) − ω1∗ (ω2 ))2 ≤ 2γV1 (ϕ(0), φe , ω1 (0), ω2 ), and ω1 (t) is bounded on [0, T ]. Moreover, according to Assumption 1 solutions φ1 (t), φ2 (t) are defined for all T > 0 (since if this would not be the case then there would exist T > 0 and t′ ∈ [0, T ] such that ω1 (t) → ∞ as t → t′ ; this, however, would contradict to that (ω1 (t) − ω1∗ (ω2 ))2 ≤ 2γV1 (ϕ(0), φe , ω1 (0), ω2 )). Given that ϕ(t) is bounded and ρ is continuous, the proposition now immediately follows from (4). Indeed, we can see that V1 (ϕ(t), φe , ω1 (t), ω2 ) − V1 (ϕ(t0 ), φe , ω1 (0), ω2 ) Z t ≤− ρ(|ϕ(τ ) − φe |)dτ. 0 Rt Hence limt→∞ 0 ρ(|ϕ(τ )− φe |)dτ = B ≤ ∞, and (without loss of the generality we can assume that ρ(·) is Lipschitz) the result follows from Barbalat’s lemma.  Frequency-tuning mechanism (3) ensures stable phaselocking of the oscillators. Yet, in order to implement this mechanism knowledge of φe and the corresponding Lyapunov function is needed. In the next section we provide an alternative that does not necessarily require availability of this information. 3.2 Asymptotic phase locking for system Σ1 Proposition 2. Consider (1), and suppose that Assumption 1 holds. Let

3.1 Globally stable phase locking Proposition 1. Consider system (1) and suppose that Assumption 1 holds. Let us also suppose that A1) for every ω2 ∈ R there exists ω1∗ (ω2 ) ∈ R such that φe is a globally asymptotically stable equilibrium of ϕ˙ = ω1∗ (ω2 ) − ω2 + f1 (φ2 + ϕ, φ2 ) − f2 (φ2 + ϕ, φ2 ), A2) a corresponding Lyapunov function V (ϕ, φe ) be known: ∂V (ϕ, φe ) ϕ˙ ≤ −ρ(|ϕ − φe |), ∂ϕ ω1 =ω ∗ (ω2 ) 1

where ρ : R≥0 → R≥0 is a strictly monotone continuous function with ρ(0) = 0.

Then solutions of system (1) in which ω1 evolves according to ∂V (ϕ, φe ) ω˙ 1 = −γ , γ ∈ R>0 (3) ∂ϕ

A3) system ϕ˙ = f1 (φ2 + ϕ, φ2 , t) − f2 (φ2 + ϕ, φ2 , t) + ξ(t), (5) where ξ(t) ∈ C 0 [0, t], be L2 [0, t)-to-L∞ [0, t) stable wrt ξ(t) and ϕ respectively. Finally, suppose that the function ω1 be defined as ω1 (ϕ, t) = −γ(ϕ + z(t)), γ ∈ R>0 (6) z˙ = f1 (φ2 + ϕ, φ2 , t) − f2 (φ2 + ϕ, φ2 , t). Then solutions of the adapting system are defined for all t ≥ 0, and ϕ(t) approaches the union of all limit sets of (5) where kξ(t)k22,[0,∞) ≤ kω1 (ϕ(0), 0) − ω2 k2 /(2γ). In particular, if φe is a globally stable equilibrium of (5) admitting a Lyapunov function ∂V V (ϕ, φe ) : [f1 (φ2 + ϕ, φ2 , t) − f2 (φ2 + ϕ, φ2 , t)] ∂ϕ ≤ −ρ(|ϕ − φe |)

11756

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

where ρ : R≥ → R≥0 is a strictly p monotone and non∂V decreasing function, with | ∂ϕ | ≤ ρ(·)D, D ∈ R>0 , and p ρ(·) is Lipschitz, then limt→∞ ϕ(t) = φe

Proof. The argument is similar to Tyukin et al. (2007). The right-hand side of the combined system is continuous, hence Peano’s existence theorem assures existence of a non-empty interval [0, T ] for which solutions of the system are defined. Let us consider the function 1 W (ω1 , ω2 ) = (ω1 − ω2 )2 . 2γ Its derivative is ˙ = −(ω1 − ω2 )2 ⇒ kω1 (ϕ(t), t) − ω2 k2 W 2,[0,T ] ≤ W (ω1 (ϕ(0), 0), ω2 ). We can therefore conclude that solutions of the combined system are defined for all t ≥ 0 (if this would not be the case then there would exist t′ > 0 such that either of φ1 , φ2 , ϕ, ω1 be discontinuous or would approach infinity; the former cannot occur by construction, and the latter is not possible too since ω1 is bounded for all t ∈ [0, t′ ], and the system Σ1 satisfies Assumption 1). Variable ϕ(t) satisfies (5), and hence it will approach the union of the omegalimit sets of (5) with kξ(t)k22,[0,∞) ≤ W (ω1 (ϕ(0), 0), ω2 ). The last conclusion of the proposition Rfollows from the ∞ analysis of V2 (ϕ, t) = V (ϕ, φ2 ) + D2 /4 t (ω1 (ϕ(τ ), τ ) − 2 ω2 ) dτ .  Propositions 1,2 formulate conditions for stable / asymptotic adaptive phase-locking in a system of coupled oscillators. These conditions include an assumption that the coupling functions f1 , f2 do not depend explicitly on the natural frequencies ω1 , ω2 . While this assumption may be plausible for systems of which the natural frequencies are sufficiently close, this assumption does not hold in general. In the next section we present a result that allows to develop the required adaptation mechanisms for systems in which the coupling between phases of oscillations depends on the natural frequencies ω1 , ω2 explicitly. 3.3 Asymptotic phase locking for system Σ2 In order to formulate the main result of this subsection would like to remind the notion of (ε, δ)-input-to-state stability Teel and Praly (1995). The system x˙ = f (x, u), f (xe , 0) = 0, x ∈ Rn , u ∈ Rm , is called (ε, δ)-input-tostate stabile at xe iff there exist β ∈ KL, γ ∈ K such that kx(t) − xe k ≤ β(t, kx(0) − xe k) + γ(ku(τ )k∞,[0,t] ) for all kx(0) − xe k < δ, ku(τ )k∞,[0,t] < ε. In what follows we will require a slightly stricter version of this property, namely we will say that the system x˙ = f (x, u, t), f (xe , 0, t) = 0, x ∈ Rn , u ∈ Rm , is called linearly (ε, δ)-input-to-state stabile at xe iff there exist a non-increasing function β, limt→∞ β(t) = 0 and a constant c ∈ R≥0 such that kx(t) − xe k ≤ β(t)kx(0) − xe k + cku(τ )k∞,[0,t] for all kx(0) − xe k < δ, ku(τ )k∞,[0,t] < ε. Now we are ready to formulate the following Proposition 3. Consider (2), and suppose that Assumption 1 holds. Let A4) system ϕ˙ = g1 (φ2 + ϕ, φ2 , ω2 , t) − g2 (φ2 + ϕ, φ2 , ω2 , t) (7) + u(t),

be linearly (ε, δ)-input-to-state stable at φe , A5) the function g1 (φ2 + ϕ, φ2 , ω1 , t) be Lipschitz in ω1 : |g1 (φ2 + ϕ, φ2 , ω1 , t) − g1 (φ2 + ϕ, φ2 , ω1′ , t)| ≤ c|ω1 − ω1′ |, c ∈ R≥0 A6) ω2 ∈ [a, b], and the values of a, b are known Consider

x1 + 1 ω1 (x1 ) = a + (b − a) , 2  x˙ 1 = γ|ϕ − φe |x2 , x (0)2 + x2 (0)2 = 1. x˙ 2 = −γ|ϕ − φe |x1 1

(8)

Then for every ω2 ∈ [a, b] there is a set Ωγ of initial conditions corresponding to the trajectories of the combined system converging to φe provided that the following condition is satisfied (b − a)(c + 1) γ· · G < 1, (9) 2 where       d k κ −1 G = βt βt (0) 1 + +1 κ k−1 1−d for some d ∈ (0, 1), κ ∈ (1, ∞). Proof. The proof is based on the results from Tyukin et al. ˙ (2008). Consider ϕ(t): ˙ = g1 (φ2 + ϕ, φ2 , ω1 , t) − g2 (φ2 + ϕ, φ2 , ω2 , t) ϕ(t) + ω1 − ω2 = = [g1 (φ2 + ϕ, φ2 , ω2 , t) − g2 (φ2 + ϕ, φ2 , ω2 , t)] [g (φ + ϕ, φ2 , ω1 , t) − g1 (φ2 + ϕ, φ2 , ω2 , t)] + 1 2 +(ω1 − ω2 ) | {z } u(t)

According to assumptions of the proposition we have that |ϕ(t) − φe | ≤ β(t)|ϕ(0) − φe | + (c + 1)kω1(x1 (t)) − ω2 k∞,[0,t] (10) whenever |ϕ(0) − φe | < δ and kω1 (x1 (t)) − ω2 k∞,[0,t] < ε. Notice that Rt sin(γ 0 |ϕ(τ ) − φe |dτ + h1 ) + 1 ω1 (x1 (t)) = a + (b − a) 2 h1 ∈ R and that ω2 can be represented as sin(h0 ) + 1 ω2 = a + (b − a) , 2 Clearly, |ω1 (x1 (t)) − ω2 | ≤ ((b − a)/2)|h(t) − (h0 − h1 )|, Rt where h(t) = γ 0 |ϕ(τ )− φe |dτ . This implies that (10) can be rewritten as |ϕ(t) − φe | ≤ β(t)|ϕ(0) − φe | (b − a)(c + 1) + kh(t) − (h0 − h1 )k∞,[0,t] 2 Letting β(0)|ϕ(0) − φe | + ((b − a)c/2)|h0 − h1 | < δ, |h0 − h1 | < ε we ensure that |ϕ(t) − φe | ≤ β(t − t1 )|ϕ(t1 ) − φe | (b − a)(c + 1) + kh(t) − (h0 − h1 )k∞,[t1 ,t] 2 holds whenever t1 ≥ 0 and h(t) ≤ h0 −h1 , and choosing the value of γ according to (9) guarantees that h(t) ≤ h0 − h1 for all t ≥ 0 (see Corollary 4.1, Tyukin et al. (2008)). The result now follows immediately from Barbalat’s lemma and Rt that limt→∞ h(t) = 0 |ϕ(τ ) − φe |dτ ≤ ∞. .

11757

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

4. DISCUSSION

Pre-synaptic

The main advantage of the results of the previous section, i.e. generality of the scope of models to which they apply, have a clear shortcoming. The shortcoming is that physical realization of these schemes may not be immediately obvious. While there is no doubt that at the level of circuits comprised of several elements adaptation schemes involving integration of phase differences may occur, the question of how the schemes can be realized by a separate physical organ needs a particular care to answer. In what follows we provide an example of a system having certain degree of biological realism and to which intuition developed earlier can apply.

Post-synaptic

Feedback

Consider a pair of neural oscillators connected through a dynamic synapse affecting the value of a threshold variable in one (or both) of the oscillators. Activity of the synapse is made dependent on the timing of spikes, and thus the connection mimics one of the documented plasticity mechanisms. For simplicity suppose that the dynamics of the system is given by ( V˙ 1 = I(V1 ) − W1 − I1 + ǫc(V1 , V2 ), (11) ˙ 1 = W∞ (V1 ) − W1 , W (

V˙ 2 = I(V2 ) − W2 − I2 − εz, ˙ 2 = W∞ (V2 ) − W2 , W

(12)

z˙ = −τ (z − z0 ) − k sin(2φ(Φ − Φc )) + λ (13) where Vi stand for the membrane potential, I, W∞ , c are continuous functions, and Wi are currents through ionic channels in the membrane. Functions W∞ , I were chosen as in the standard Rowat-Selverstone model Rowat and Selversone (1993). Variables I1 and I2 are constants modelling excitation of the membrane, ϕ is the phase difference, ǫ, ε, τ, k > 0, z0 ∈ R are parameters. The function c is instantaneous, fast component coupling between cells: 1 c(V1 , V2 ) = (V1 − V2 ). (Θ−V 2 )/ksyn 1+e Variables Θ, ksyn characterize midpoint and slope of the synaptic activation respectively. Variable z is a dynamic component of the coupling. System (11) stands for the postsynaptic neuron, and (12) models the presynaptic one (see Fig. 1). Variable Φ is the spiking relative phase defined as Φ=

tpost (i) − tpre (i) , 0 < Φ < 1, Tpre

(14)

where tpre (i) is the time instant of the i-th presynaptic spike, and tpost (i) is the time of the closest response spike generated in the postsynaptic neuron. Hence the variable Φ may be viewed as a sampled version of the relative phase, ϕ. Variable Φc is the reference post-to-pre- synaptic phase shift, and λ is a baseline of neural excitation. We assume that solutions of (11), (12) exist and bounded. The assumption will obviously hold for semi-passive systems, provided of course that the growth rate of the function c(V1 , V2 ) is constrained by an appropriate condition; and Steur et al. (2009) show that most canonic models of neural oscillators are indeed semi-passive.

Fig. 1. Connection diagram for (11)–(13). Postsynaptic cell, (11), is connected with a presynaptic one, (12), via dynamic synapse, (13). With regards to expressing the phase dynamics of (11), (12) phases of oscillators can be formally defined as functions φ1 = φ(V1 , W1 , I1 ) and φ2 = φ(V2 , W2 , I2 ) such that ∂φ1 ˙ ∂φ1 ˙ V1 + W1 = ω(I1 ) + ǫg(φ1 , φ2 , t) ∂V1 ∂W1 ∂φ2 ˙ ∂φ2 ˙ V2 + W2 = ω(I2 + εz) ∂V2 ∂W2 where g(·) reflects the influence of instantaneous coupling c(V1 , V2 ), and ω(·) is the function modelling natural frequency of oscillations. In the first approximation, if the value of ǫ is negligibly small, dynamics of the relative phase, ϕ = φ1 − φ2 , can be described  as ϕ˙ = ω(I1 ) − ω(I2 + εz) z˙ = −τ (z − z0 ) + k sin(2φ(Φ − Φc )) + λ.

The function ω(·) for the chosen parameter values of (11) is shown in Fig. 2. Simplifying the model further by replacing discrete variable Φ with its continuous substitute, ϕ, and denoting ω(I1 ) − ω(I2 ) = ω0 , ω(I2 + εz) − ω(I2 ) = f (z), we arrive at  ϕ˙ = ω0 − f (z) Σ3 : (15) z˙ = −τ (z − z0 ) + k sin(2π(ϕ − Φc )) + λ where f is a monotone, continuous, and strictly increasing function with f (0) = 0, and ω0 is the variable reflecting the difference in the natural frequencies of oscillators (11), (12). Let ϕ∗ , z ∗ be an equilibrium of (15), and let the function f be differentiable. Then the equilibrium is (locally) asymptotically stable as long as df /dz > 0 at z = z ∗ and cos(2π(ϕ∗ − Φc )) > 0. The conclusion follows

11758

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

V1 , V2

Fig. 2. Frequency of oscillations in (11) for ǫ = 0 as a function of I1 . 1.0 0.5 0.0 -0.5 -1.0 -1.5 0

200

400

600

800

1000

1200

1400

F

1.0 0.5 0.0 0

200

400

600

800

1000

1200

time

Fig. 3. Upper panel: trajectories of variables V1 (t), V2 (t) as functions in a system of two neural oscillators in a plane. The value of λ was set to zero. Spikes from a post-synaptic oscillator are shown by solid lines, and spikes from the pre-synaptic one are depicted by dashed line. The relative phase ϕ is approximated by Φ, (14). For t ≤ 600 the synaptic connection is turned off, and relative phase of the oscillators (Lower panel) is drifting. At t = 600 the synaptic connection, (13), is activated. We can see that, after a period of transient dynamics, phases of two oscillators are converging to the phase-locked state. from the analysis of the time-derivative of the following Lyapunov candidate: Z z V (ϕ, z) = (f (s) − f (z ∗ ))ds+ ∗ Z ϕz k sin(2π(s − Φc )) − sin(2π(ϕ∗ − Φc ))ds. ϕ∗

And this is indeed the behavior we observed in simulations (see Fig. 3). The phenomenon is not limited to systems (11) – (13); it can be straightforwardly generalized to substantially more complex and realistic conductancebased models (see e.g. Izhikevich (2007)) in which the frequency of spikes can vary in a relatively large interval as a function of the threshold parameter. An interesting and unexpected observation is that regulatory scheme (15) closely resembles structurally abstract adaptive phase locking schemes considered in Sections 3.2 and 3.1. The differences being the presence of extra stabilizing terms k sin(ϕ − Φc ) and −τ (z − z0 ) in (15) as compared to (6) and (3) respectively.

Fig. 4. Φ∗ bifurcation diagram (for varying k). The value of τ is 0.01, I2 = −0.05, ǫ = 0.008 the reference relative phase, Φc , is 0.6. We would like, however, to notice that despite these similarities the original scheme is not entirely adaptive in the usual control-theoretic sense: phase locking in system (11), (12) generally occurs with an error. Indeed, the value of ϕ∗ satisfies: τ (f −1 (ω0 ) − z0 ) − λ sin(2π(ϕ∗ − Φc )) = . k Thus if the values of I1 and I2 are not identical, and if z0 , λ 6= 0 then ϕ∗ 6= Φc + n (n is the period). Hence a compensatory mechanism is needed if higher accuracy of locking is desired. Since the values of I1 , I2 , τ , and z0 are assumed to be unknown, a rather restricted range of possibilities for compensating mismatches between ϕ∗ and Φc is available. The most obvious one is to increase the value of k in the dynamic component of the synaptic connection. The alternative is to modify the value of λ adaptively so that the influence of τ (f −1 (ω0 ) − z0 ) is canceled. Both possibilities have been studied. We found that operational range of the first possibility, i.e. increasing the value of k in (12), is restricted by relatively small values of k. This is because (largely due to that Φ is discrete) equilibrium ϕ∗ , z ∗ loses stability if k exceeds some critical value (see Fig. 4). Thus if further fine-tuning of ϕ is required then adaptation of the baseline parameter λ is needed. One of the most simplest ways to implement such an adaptation is to follow the idea very similar to that discussed in Section 3.3 and analyzed in Tyukin et al. (2008). The rationale behind this choice being that the right-hand side of the equations governing relative phase dynamics is not known. Thus making use of the gradient of derivative of the Lyapunov function (such as e.g. in Section 3.1) or exploiting knowledge of the right-hand side explicitly (Section 3.2) in an adaptation scheme is not always plausible. Instead we suppose that the relative phase variable, ϕ, locally satisfies the following inequality:



ϕ(t + T )



≤ β(T ) ϕ(t)

z(t + T )

z(t)

+ c · max λ(t) − τ (z0 − f −1 (ω0 )) , t∈[t,t+T ]

where β(·) is a decreasing strictly monotone function such that limT →∞ β(T ) = 0, and λ(t) − τ (z0 − f −1 (ω0 )) is a

11759

18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011

and/or stability of the phase-locked (synchronous) states is only local. These results, forming prototypical scenarios for adaptive phase synchrony, may serve as simple models of phase adaptation and tuning in networks of neural oscillators and cultures provided that the corresponding adaptation mechanisms can be realized in local neural circuits in the network. We have shown how one of these mechanisms can be implemented as a component of dynamic coupling in a model of coupled neural oscillators. Finding additional mechanisms and deriving their plausible realizations in neural circuits is the subject for future study. Fig. 5. Adaptive compensation of phase error via (extracellular) adaptation feedback. Parameter values: τ = 0.01, I2 = −0.05, ǫ = 0, Φc = 0.1, k = 0.0005, γ = 0.001. term characterizing the amplitude of the relative phase’s fluctuations around desired values at t → ∞. The value of λ(t) is set to vary according to the simple adaptation mechanism below: λmax − λmin λ = λmin + (1 − sin(ζ)) 2 (16) ζ˙ = γ|Φ − Φc |, γ ∈ R>0 , λmax > λmin , λmax , λmin ∈ R. This adaptation mechanism is in essence a slow fluctuation of the excitation thresholds. Frequency of these fluctuations increases if absolute values of the relative phase are far away from the desired ones; the fluctuations’ frequency slows down when relative phase approaches Φc . Asymptotic properties of the system can be analyzed similarly to our earlier assessment of the dynamics of (2) (see Proposition 3), and due to the space limitation we omitted such analysis here. Numerical simulations of (11), (12) with dynamic synapse and excitation baseline adaptation (16) is illustrated in Fig. 5. As the figure illustrates, activation of adaptation feedback ensures that the value of relative phase, ϕ, is slowly approaching the desired reference value Φc asymptotically. The picture remains the same for a broad range of Φc . The adaptation process above can be thought of as slow fluctuations of “state of extracellular matter”. At the present level of biophysical details used in our simulations we could not associate such process explicitly with particular extracellular molecular cascades. Nevertheless, we can speculate that certain characterizations of the process (e.g. low strength, relatively slow time scale, and integration effect) are quite similar to the influence of glia and extracellular matrics on synaptic transmission described in the literature. 5. CONCLUSION In this paper we considered the problem of adaptive phase synchronization in a system of two phase oscillators with nonlinear coupling. In particular, we studied existence and asymptotic properties of adaptation mechanisms aimed at self-tuning of natural frequency in one of the oscillator so that phase synchronization or phase-locking could occur in the system. Three different mechanisms were presented and analyzed including settings in which the phase couplings could depend on the natural frequencies nonlinearly,

ACKNOWLEDGEMENTS This work is supported by the Royal Society UK-Russia International Research Grant, and Russian Foundation for Basic Research (Grants No 09-04-12254, 09-02-92611, 0904-97090, 10-01-00690). REFERENCES I.I. Blekhman, A.L. Fradkov, H. Nijmeijer, and A.Yu. Pogromsky. On self-synchronization and controlled synchronization. Systems and Control Letters, 31:299– 305, 1997. E. M. Izhikevich. Dynamical Sytems in Neuroscience. MIT Press, 2007. V. Kazantsev, I. Mukhina, Y. Zaitsev, and A. Pimashkin. Generation of repeated patterns of excitation transmission and cirulation in the dynamics of neuornal circuits of brain: data analysis and modelling. In Proc. of MEPHI-2010, January, Moscow. 2010. I. Miroshnik, V. Nikiforov, and A. Fradkov. Nonlinear and Adaptive Control of Complex Systems. Kluwer, 1999. V. Pasquale, P. Massobrio, L. Bologna, M. Chiappalone, and S. Martinoia. Self-organization and neuoronal avalanches in networks of dissociated cortical neurons. Neuroscience, 153:1354–1368, 2008. A. Pikovsky, M. Rosenblum, and J. Kurths. Synchronization: a universal concept in nonliear sciences. Cambridge University Press, 2001. P.F. Rowat and A.I. Selversone. Modelling the gastric mill central pattern generation of the lobster with a relaxation oscillator network. Journal of Neruophysiology, 70(3):1030–1053, 1993. E. Steur, I. Tyukin, and H. Nijmeijer. Semipassivity and synchronization of diffusively coupled neural oscillators. Physica D: Nonlinear Phenomena, 21:2119–2128, 2009. S. Strogatz. SYNC: The Emerging Science of Spontaneous Order. Hyperion, 2003. A. Teel and L. Praly. Tools for semiglobal stabilization by partial state and output feedback. SIAM J. Control Optim., 33:1443–1488, 1995. I.Yu. Tyukin, Prokhorov, and C. van Leeuwen. Adaptation and parameter estimation in systems with unstable target dynamics and nonlinear parametrization. IEEE Trans. on Automatic Control, 52(9):1543–1559, 2007. I.Yu. Tyukin, E. Steur, H. Nijmeijer, and C. van Leeuwen. Non-uniform small-gain theorems for systems with unstable invariant sets. SIAM Journal on Control and Optimization, 47(2):849–882, 2008.

11760