Limit Cycles in Coupled Biological Systems

Copyright © IFAC Theory and Application of Digital Control New Delhi . India 1982

LIMIT CYCLES IN COUPLED BIOLOGICAL SYSTEMS R. Balasubramanian* and D. P. Atherton** .Department of Electrical Engineering, University of New Brunswick, Fredericton, New Brunswick, Canada •• School of Engineering and Applied Sciences, University of Sussex, Falmer, Brighton, Sussex, UK

Abstract. The paper investigates the effect of coupling on the synchronized oscillations of biological systems. The individual oscillators of the system are modelled as relay oscillators with the advantage that exact determination of the synchronized oscillations is possible using an extension of Tsypkin's method. An interesting result is the determination of asymmetrical oscillations in two symmetrically coupled oscillators. Keywords. systems.

Limit cycles ; relay control; biology; stability; multivariable

INTRODUCTION

estimate of the oscillation rather than the exact solution (Atherton, 1974) . Not surprisingly therefore where several limit cycle solutions exist in a coupled system, even when the filtering may be quite good, the OF method may give inaccurate indications of the stability of an oscillation.

The occurrence of oscillations in biological systems is well known and considerable research effort has been devoted to their investigation (Wever, 1979; Stucki, 1978; Rapp, 1976a) . Various biological oscillators, with different individual frequencies of oscillation, tend to interact with each other through biochemical and biological couplings. In general, these couplings manifest themselves in a 'synchronized rhythm ' wherein individual oscillators are forced to follow the frequency, normally, of the highest frequency oscillator . An investigation of the coupling requirements under which the oscillators will synchronize to a common mode, often called the 'entrainment conditions ' , is therefore of importance. Several factors make the investigation of these oscillations difficult. Good mathematical models are often not available so that a suitable description for the individual oscillators must be found. In this paper we choose to represent the oscillators by relay oscillators both for analytical reasons, which are discussed further, and also because there is strong evidence that on off type behaviour occurs in many biological systems (Kitney, 1979) . Previous analyses of coupled biological oscillators have used the harmonic balance or describing function (OF) method (Rapp, 1976b; Linkens , 1980), which being an approximate approach has disadvantages especially when a number of coupled oscillators is considered . The difficulty is not only that the OF method gives an approximate solution for an oscillation but to assess the stability of the oscillation an approximate method is used to determine the stability of the variational equation which contains the DF

643

The choice of relay oscillators for the models means that the multivariable extension of the exact analysis method of Tsypkin for relay systems (Atherton, 1975) can be used to evaluate synchronous limit cycle modes. In addition , an exact stability criterion has recently been derived (Balasubramanian, 1981) so that it is also possible to determine accurately which oscillations are stable. In this paper we briefly review the above theory and then consider two specific examples of its application . The first example investigates the entrained oscillations produced by the coupling of models of the cardiovascular and respiratory systems . The second example considers synchronization of two oscillators which model the myoelectric activity of the human large intestine .

ANALYTICAL APPROACH Figure 1. shows a block diagram for several relay coupled oscillators. An exact analytical procedure (Atherton, 1966), which is an extension of the original work by Tsypkin (1955), starts by assuming that when a synchronized oscillation exists the outputs, y (t), of the individual ideal relays will be s~mmetrical square waves with the same fundamental frequency but which switch at different time instants . If the square wave output th of the j relay switches from -h to h at time t -- a j /w then it can be written as the

R. Balasubramanian and D. P. Atherton

644

n

Fourier series y/t)

L

= 1:

When this signal is the input to the transfer th function Gij(s), the i-j element of the transfer matrix, G(s), the Fourier series of the output c . (t) is given by iJ

1:

cij(t) =

(4h/ nTI)gi ' sin(nwt - na. J J

n=l

(2)

+ 1J!ij) and that of its derivative by 00

L (4hw/TI)g cos(nwt - na . J n=l ij + 1J!ij)

(3)

where gi' and 1J! .. are now functions of nw and J

1J

we have assumed lim sG .(s) = O. When this s '>OO iJ condition is not satisfied a correction must be made Eq. (3) to account for the discontinuity in Cij(t) . If we define an A locus (Atherton, 1975), a function of the frequency, w, and an angle 8, for the transfer function G(s) by A (8, w) = Re A ( 8 , w) + j Im A ( 8 , w) G G G Re A ( 8,w) = G

1:

j=l

(1)

(4h/ nTI) sin(nwt - na . ) J n=l

(4)

ing say a

= 0, and the inequality must then l be checked to test that the solution is valid. This is not as difficult as it may appear at first, since closed form solutions can be obtained for Eqs. (5) and (6) for various transfer functions G(s) (Atherton 1975) . In general, however, the .resul t ing nonl inear algebraic equations given by Eq. (12) do not yield simple analytical solutions and must be solved numerically. One should, in principle, also check that the x (t) do not pass through zero at other than i the above switching times, but this does not appear to be necessary when the elements of G(s) have 'smooth' frequency responses (Atherton, 1981) . 2.1 Stability of a limit cycle solution

An exact method for determining the stability of an oscillation in a feedback loop with a single relay characteristic has recently been given by Balasubramanian (1981). The method extends directly to the multivariable system of Fig . 1. Writing a state space description for the system, we have

VG(nW)sin n 8

z

(l / n){V (nw)cos n 8 G

(6) where the summations are over odd terms only and G(jnw) = UG(nw) + jVG(nw) , then it is easily shown that Eqs . written

(2) and (3) can be (a. - wt, w) ij

(7)

J

= (4hw/TI )Re A th

Gij

( a . - wt, w) (8) J

relay , xi(t), is given

n

=-

L c .. (t) j=l 1J

(9)

and for the relay to give the assumed output square wave the switching conditions

o

(10)

xi ( a i /w ) > 0

(11)

xi ( a i / W)

(14)

n(x)

y

and it is easily shown that if the system has a periodic solution with x(t) = x*(t) then the A matrix, A(t), of the variational equation is A + Bn'(x*)C, where n'(x*) = diag {n'.(x .*) } with n'. (x . ) = dn . (x . ) / dx .. Since J J J J J J J n.(x . ) is an ideal relay characteristic, J

and

x (t) i

= Cz

J

n~(x . )

(4h/ TI)lm AG

Now the input to the i by

Az + By

x

-c

(5)

1:

( 13)

Eq. (12) and inequality (13) exist for all i = l,2 ... n , so that the former can be used to solve for the n unknowns w and a i' assum-

n=1(2)

n=1(2)

Re A ( a - a ,w ) < 0 Gij j i

is zero except over the infinitesimal J J time, 6 t, it takes the relay to switch. Thus the matrix , A(t) , is equal to A except at the switching instants . It is known (Willems , 1970) that the null solution of a time varying dif·ferential equation with a piecewise constant periodic matrix , A(t), is uniformly asymptotically stable if and only if the eigenvalues of the matrix m

(15)

TI e x p A. (t . - t _ ) i l i=l 1 1

have magnitudes less than unity , apart from one which will have unit magnitude . A. is the 1

must be satisfied. Substituting Eqs. (7), (8) and (9) in (10) and (11) respectively, gives the equation

switching at time to is exp BKiC where Ki

n

L j=l

constant value of A(t) in the time interval ti - t _ and T = tm - to is the period of i l A(t). The contribution to expression (15) frOm relay , i, with input xi(t) which causes

Im A

G

(a

ij

and the inequality

j

a. , w) 1

o

(12)

diag {k .} with k

= 0 for i ~ j and k

j 2h/ xt(t ) for i = j . When the relays have o symmetrical square wave outputs expression (15) may be taken over the half period , T/ 2. J

j

Limit Cycles in Coupled Biological Systems Thus for a system with two relays, having a limit cycle solution with the relay outputs switching positive at times tl and t , 2 respectively, expression (15) becomes exp BKIC exp A(t -t ) exp BK C 2 l 2 exp A(t -t +T/ 2) (16) l 2 assuming t2 > tl and t -t < T/ 2 . A similar 2 l expression follows directly for any number of relays. The stability of a limit cycle solution given by Eq. (12) can thus be evaluated, for a two relay system, by determining the eigenvalues of expression (16) with the solution parameters K , K , tl and t2 substituted. Kl and l 2 K2 are seen to depend upon the derivative of the relay input at the switching instant, which from Eqs . (8) and (9) and inequality (13), is seen to be (4hW/TI) times the left hand side of inequality (13), this is evaluated with the limit cycle solution, for relay i switching at time xi /W o

APPLICATIONS

645

parameters of the transfer matrix G(s) as input . For a given value of ~ l the mini mum value of ~2 required for stable entrained oscillations is shown in Table 1 together with the corresponding simulation results . Table 1.

Conditions for synchronous rhythm in Example 1 Minimum value of ~2 Simulation Theorv

~l

2.51 2 . 29 2.08 1.87 1.68 1.63

0.0 0.1 0.2 0. 3 0.4 0.42

2.50 2 . 30 2 . 03 1.85 1.67 1.65

Fig. 3 shows the region of synchronization graphically . For ~ > 0.42 no theoretical reFults for a synch~onized limit cycle could be found for any value of ~2 whilst in the analogue simulation complex oscillations were found for this condition. Fig. 4 shows the synchronized limit cycles obtained in the simulation for ~l = 0.1 and ~2 = 2.3 .

Two applications to biological systems of the analysis described in the previous section are given below

DF analysis produced poor results , no doubt due to the nonsinusoidal nature of the waveforms shown in Fig . 4.

3 . 1 Example 1

3 . 2 Example 2

Here we consider the conditions of coupling required for entrained oscillations in models of the cardiovascular and respiratory systems. In terms of Fig. 1 the relays have h = 1 and

The effect of symmetrical coupling on two identical relay oscillators is investigated. Previous studies (Linkens, 1980) on this interesting biological problem have used Van der Pol oscillator models and harmonic balance analysis. Here the relays of Fig . 1 are again assumed to have h = 1 and the transfer function matrix, G(s) , is given by

G(s)

When

~l

= ~2 = 0

G(s) there is no coupling between

the systems and the transfer function of the -S T l cardiovascular loop is Gl(s)e with Gl(s)

=

22(1+12s) (1+92s)(1+2s)

and Tl

=

3 . 3 secs.

The transfer function of the respiratory loop -s T 2 G (s)e bas G (s) = Gl(s) and T2 = 1.2 sec . 2 2 The limit cycles in the uncoupled loops obtained by analogue simUlation are shown in Fig. 2 . The time delays , since they act on square waves, were simula-ted exactly Using parallel logic. The measured frequencies of the limit c y cles were within experimental accuracy of the exact calculated values of 0.6727 and 1 . 5747 rads / sec . for the respiratory and cardiovascular systems . The corresponding DF values are 0 . 643 and 1 . 534 rads / sec . To investigate the region of entrained oscillations, solutions for Eq. (12) were found for various combinations of ~ l and ~ 2 using a general computer program which requests the

= 1-0.3~ (1+s)

(l / S

(18)

~

The non-minimum phase nature of the elements of G(s) is typical of many biological OSCillators. With no coupling, i . e . ~ = 0, the exact frequency 6f oscillation is 0.7641 rads / sec and the DF solution 0.791 rads / sec. The limit cycle waveform obtained from an analogue simulation is given in Fig . 5 and is seen to be almost sinusoidal . DF analysis of the coupled system yields a synchronized limit cycle solution for five values of e , where e/w is the time difference between the positive switchings of the two relays. These are for e = 0 and ±TI when the corresponding limit cycle frequency is given by (19)

e

with the + sign for = 0 and the - sign for e = ±TI , and for e = ±TI/2 when the limit cycle frequency is 0 . 791 rads / sec. independent of ~. The solutions for e = ±TI / 2 surprisingly, since the system is symmetrical, yield unequal oscillation amplitudes at the inputs to the two relays . Analysis of the system using

R. Balasubramanian and D. P. Atherton

646

Eq . (12) yields seven exact solutions for e = 0, ±1T, ± By and ± 8 2 , where of course 8 1 and 8 depend upon~ . For comparison Table 2 shows 2 toe = 0 and ±1T solution frequencies obtained by the exact and OF methods .

e

Table 2

Comparison of OF and exact solutions

e =0 ~

0.0 0.2 0.4 0.6 0.8 1.0

e = ±1T

CONCLUSIONS

OF

Exact

OF

Exact

0 . 791 0.917 1.101 1. 351 1.619 1 . 805

0.764 0.878 1.034 1. 231 1 . 435 1.607

0 . 791 0.701 0.637 0.585 0 . 545 0.513

0 . 764 0.679 0 . 6l4 0.563 0 . 521 0.486

The stability of the exact solutions evaluated using Eq . (16) showed that the solutions for e = ±1T and e = ±8 2 were unstable for ~ in the range 0 ~ ~ ~ 1 considered . The limit cycle frequencies for e = ±1T decreased from 0.7641 as ~ increased and had values extremely close to those e = ±8 , which are given in Table 3. l The frequencies for e = ± 8 increased slowly 2 from 0 . 7641 as ~ increased. All the solutions for e = ±8 and those for e = 0 with ~ < 0.66 l were found to be stable and corresponded to similar results found in the analogue simulation. These results are compared in Table 3. Table 3 . Comparison of exact solutions and simulation results

e= e

~

Exact

= 0

l

Exact

Sim .

0.1 0.816 0.816 0.2 0.878 0.873 0.3 0 . 950 0.943 0 . 4 1 . 034 l.013 0 . 5 l.128 l.112 0.6 1 . 231 l.2l0 0.7 No stable 0.8 oscs. for 0.9 \.1 >0.66 f~0.63 l.0

±8

w 0 . 729 0.692 0.657 0.625 0.596 0.570 0.546 0.525 0.506 0.488

8 1 3 . 856 3 . 748 3.664 3 . 601 3 . 552 3.513 3.482 3.458 3.437 3.420

Sim. w 0.731 0.694 0.668 0.633 0 . 593 0 . 571 0.543 0.526 0.503 0.486

8 1 3 . 86 3.76 3.65 3.61 3.55 3 . 52 3 . 49 3.45 3.43 3.37

The values of w and e in the simulations were accurately measured using a counter to obtain the times between relay switchings . The e = ±8 l solutions are unsymmetric in the same sense mentioned earlier since they yield unequal oscillation amplitudes at the relay inputs. Figures 6, 7 and 8 show various limit cycles recorded from the simulations . For ~ = 0.3 stable oscillations exist for e = 0 , as shown in Fig. 6, and for = ±3.65 where the waveforms of Fig. 7 clearly show the unsymmetric behaviour . For ~ = 0.8 the only solution is the unsymmetric case and the waveforms are shown in Fig . 8. It was also found from the simulation that from zero initial conditions the low frequency

e

unsymmetric oscillation behaviour evolved and initial conditions were required to excite the higher frequency, e = 0, oscillation . In addition, the sign of 6 1 could be changed by variation of the initial conditions .

In this paper it has been shown that where individual oscillators can be represented as relay oscillators it is possible to determine accurately conditions for synchronized oscillations when they are interconnected . The approach is particularly useful for the study of biological rhythms which can often be modelled in this manner and two specific applications have been given . The fact that the stability of a synchronized oscillation can also be evaluated by the method presented is extremely important, since several oscillation solutions may exist when oscillators are interconnected but only a limited number of the solutions will be stable .

REFERENCES Atherton, O.P . (1966) . Conditions for periodicity in control systems containing several relays . Paper 28E 3rd IFAC Congress, London . Atherton, O. P. (1974). The stability of oscillations. Automatic Control Theory and Applications, ~, 47-49. Atherton, O. P . (1975). Nonlinear control engineering , Van Nostrand Reinhold, London . Atherton, D.P. (1981) . Stability of nonlinear oscillations, Research Studies Press, John Wiley, Chichester. Balsubramanian, R. (1981). Stability of limit cycles in feedback systems containing a relay, Proc . lEE , l28D, 24-29 . Kitney, R.I. (1979). A nonlinear model for studying oscillations in the blood pressure control system. J.Biomed.Eng., l, 89-99 . Linkens, O. A. (1980). Pulse synchronization of intestinal myoelectric models. IEEE Trans . Biomed. Eng., 27, 177-186 Rapp, P. (1976a) . Mathematical techniques for the study of oscillations in biochemical control loops . Bull . Inst . Math . Applic . , ~, 11 ..,. 21 Rapp , P. (1976b) . Analysis of biochemical phase shift oscillators by a harmonic balancing technique . J . Math. BioI. , 3, 203-224 . Stucki, J . W. (1978). Stability analysis of biochemical systems - a practical guide. Prog . Biophys. Molec. BioI ., 33, 99..,.187 . Tsypkin, Ya. Z. (1955). Theory of relay control systems . State Press, MOECOW . Wever, R. A. (1979) . The Circadian system of man. Springer-Verlag, Berlin . Wi11ems, J . L . (1970) . Stability theory of dynamical systems . Nelson, London.

647

Limit Cycles in Coupled Biological Systems

Entrainmeat Region Transfer Matrix

G(s}

----0------0. - ----0-

o c

o

-----0.0

Theoretical Simulation

n

o 0.0

0.1

0.]

0.2

0.4

III Fig. 1. Coupled oscillator block diagram Fig. 3. Region of synchronized oscillations _

Fig. 2 . Uncoupled limit cycles Fig. 4. Synchronized limit cycles for ~ 1 = 0.1 and ~2 = 2.3

R. Balasubramanian and D. P. Atherton

648

I

i

I

Fig . 5. Limit cycle with

~

0

Fig. 7. Synchronized oscillations for ~ = 0.3 with e = 3.65 rad.

Fig. 6. Synchronized oscillations for ~ = 0.3 with = 0°

e

Fig. 8. Synchronized oscillations for ~ = 0.8 with = 3.45 rad.

e