J. them. Biol. (1979) 77, 235-251

Modulation

Department

Analysis of Forced Non-linear Biological Modelling D. A. LINKENS University

of Control Engineering,

(Received 12 July 1938, and in rezlisedform

Oscillators for

of ShefJield, England 18 September 1978)

The synchronization of self-sustained oscillations by a forcing oscillation is of interest in a number of biological models. It has been considered for circadian rhythm modelling, heart-rate variability studies and forced breathing experiments. Outside the range of synchronization, conditions of almost-entrainment occur in which changes in amplitude and/or frequency are apparent. It is shown in this paper that such conditions can be analysed as modulation phenomena using the analytical method of harmonic balance. The degree of non-linearity in the self-sustained oscillation affects the nature of modulation, in that increasing distortion gives a trend towards frequency rather than amplitude modulation. The analytical results compare favourably with spectral analysis of simulated oscillators. 1. Introduction The ability of non-linear self-sustained oscillations to be synchronized by an external unidirectional forcing oscillatory function has been known since the classic circuit experiments by Appleton & van der Pol(1922). The analysis of these so-called “silent” regions was performed using the now famous van der Pol non-linear differential equation (van der Pol, 1927, 1934). Van der Pol’s equation has been used and modified for a number of biological modelling studies. In circadian rhythm research it has been used as the basis of the endogenous free-running oscillator (Wever, 1965). It has been studied as a model of the human thermal vasomotor control system with forcing functions applied thermally to the hand (Kitney, 1975). Because breathing is also of a non-linear oscillatory nature, van der Pal’s equation with external forcing functions has been considered as a model for forced experiments (Fincham & Liassides, 1978). Although respiration

synchronization or entrainment is feasible in each of these examples a common feature is that of almost-entrainment, and it is to the analysis of these conditions

that this paper is devoted. 235

0022-5

193/79/070235

+ 17 SO2.00/0

0

1979 Academic

Press Inc. (London)

Ltd.

236

D.

A.

LINKENS

Almost-entrainment conditions have frequently been encountered in circadian rhythms (e.g. Aschoff, 1965: Swade & Pittendrigh, 1967). The phenomenon has been variously described as “relative co-ordination” and “oscillatory free-run”. Wever (1972) used the term “relative co-ordination” for weak forcing conditions when “beating” effects predominate. He used the phrase “relative entrainment” for medium forcing conditions just outside the region of synchronization. It is in this region that changes in amplitude and phase occur periodically and produce modulation phenomena. In the next section some examples of the frequency spectra of heart-rate variability signals will be given to illustrate these points. The analysis of relative co-ordination in which beating effects are caused by the addition of two frequency components (referred to as “combination tones”) has been performed by van der Pol (1934). The fluctuations in amplitude caused by this condition have been demonstrated by Wever (1972) for small non-linearities, and the spectra shown to contain two major frequency components. Such an approach is not valid for the relativeentrainment region nor for larger degrees of non-linearity. A further objection to the classic combination tones analysis is that the individual frequencies are taken to be that of the forcing function and the self-sustained frequency with zero forcing. It will be shown experimentally and theoretically that the self frequency is pulled towards that of the forcing function during relativeentrainment. To account for both amplitude and frequency fluctuations the oscillation will be assumed to have three frequency components with unknown separation between them. 2. Modulation

in Heart-Rate

Variability

Signals

The human heart-rate is not constant with time and can show oscillatory behaviour reflecting different control mechanisms depending on the frequency of oscillation observed. The heart-rate variability is normally derived from the QRS complex of the electrocardiogram in terms of the R-R interval. The examples shown here use a method referred to as lowpassfiltered event series (Hyndman & Mohn, 1973) and instrumented by Coenen (see Rompelman et al.. 1977). In Fig. 1 the spectrum produced by an FFT of 5 12 data points at a sampling rate of 2 Hz clearly shows three major components. The central frequency corresponds to the self oscillation within the control loop involving blood pressure and which affects the heart-rate (Sayers, 1973) while the lower frequency is an injected component from the fifth harmonic of the square wave thermal stimulus. The other component is “reflected” about the central component as a “sideband”. As shown in amplitude form this spectrum could

MODULATION

.ANALYSIS

FOR

BIOLOGICAL

MODELLING

237

H-R

I

1 0.0

I

0.4

0.2 HZ

llll~llll~l~l~~~ 0

12

24 CPm

(b)FFT

1. Heart-rate variability data and FFT spectrum showing relative entrainment onto the 5th harmonic of square wave thermal stimulus. The heart-rate data were obtained by measurement of R-R intervals, while the thermal stimulus was produced by the subject placing a hand alternately in water baths at different temperatures. FIG.

indicate either moderate amplitude and/or frequency modulation of a sinusoidal nature. From the data itself it can be seen that fluctuations in frequency were the predominant feature in this case. It will be seen later that this is a feature of oscillations with higher non-linearity factors and in fact, the appreciable distortion in the heart-rate waveform can be seen in the data. The apparent entrainment of the self oscillation onto a stimulus harmonic can be seen in Fig. 2, which also illustrates the point that biological rhythms are often labile in frequency and amplitude. This is indicated by the broad spread of the FFT spectrum. It is further illustrated in the unstimulated spectrum of Fig. 3, in which endogenous modulation sidebands are visible on either side of the major self-oscillation peak. Because of the labile nature of the rhythms the occurrence of relative-entrainment is a common feature when using fixed frequency stimuli in physiological experiments. For example, in

238

II.

4.

LINKENS

(a i Data

H-R

Shm l”.rnL

I 0.2

I 0.0

I 0.4 Hz

(b)FFl FIG.

2. Labile

entrainment

of heart-rate

variability

onto

5th harmonic

of stimulus

Fig. 2(b) of Kitney (1975) the spectrum of an entrained digit blood flow signal does in fact show the characteristic sidebands of a system having small modulation effects. Kitney & Rompelman (1978) have used the term “metastable entrainment” for the case of large phase and amplitude fluctuations which give rise to significant sidebands in the frequency spectrum. 3. Derivation

of Harmonic

Balance Equations

The system model comprises a van der Pol oscillator with sinusoidal forcing function giving the following differential equation .i-&(a2-X2)i+O~1X-Fcos(02t+~)

= 0,

(1)

where the amplitude of the self-oscillation is related to uz (i.e. approximately equal to 2a) and the frequency is related to wfl (i.e. approximately equal to w1 1). a2 is taken in the succeeding analysis to be unity. The phase angle of the

MODULATION

ANALYSIS

FOR

BIOLOCilCAL

239

MODELLING

(a) Data

1 0.0

I

I

0.2

0.4 Hz

LIIIIII11111111’ 0

12

24 cm

(b)FFT FIG. 3. Unstimulated modulation.

heart-rate

variability

spectrum

showing

small

amount

of frequency

forcing function is taken to be #Jrather than zero, to facilitate the presentation of modulation vector results. For small values of the non-linear factor E the solution to equation (1) may be assumed to have three spectral components to account for modulation phenomena as follows x = A cos(O~t+X)+BCOS(W1t+~)+CCOS(O~t+~).

(2)

In equation (2) B represents the major component with frequency o1 related to, but not necessarily equal to, o1 1. C represents the injected term at the stimulus frequency 02, and A is the “reflected” sideband component with frequency o3 given by W) = 204 -c$.

(3)

Differentiating (2) twice with respect to t and considering only the steady-state conditions (i.e. ignoring time derivatives of amplitudes and phases) gives 1 = -Am3 .iZ= -Am;

sin (oJt+a)-BoI cos (c+t+x)-Lko:

sin (c~~r+fi)-C~~

sin (c+t+~),

(4)

cos (o,t+/l-Cw;

cos (w2f+y).

(5)

240

I).

A.

L.INRENS

Using the method of harmonic balance (e.g. Lawden, 1959, p. 349) the assumed expressions (2), (4) and (5) are substituted into the governing equation (1). Simplification of the non-linear term expressions into linear Fourier series leads to a number of harmonic and sum and difference terms. All terms different to (jj3. OJ~ and (fJ1 are ignored. and then the method of harmonic balance equates the various sinusoidal functions dependent on time to leave a set of algebraic functions instead of the original differential equation. The following relationships are thus obtained (1) equating coefficients

of cos Ojj t 2

A(o~,-oJ~)COsa+eAo~J

i -ET

B2C

co3 sin (2p--y)

(2) equating coefficients

sinr

I-!j(A’+B”+Cz)f%

= 0.

of sin o3

(6)

t

-A(o:r-o:)sinrX+aAw3

A2

I-!2(Az+B2+C.‘)+1m-

cos.~

i -&c0S(2~+0, (3) equating coefficients B(w:,

- ~~)cos/?+EBco~

-~EABC~,

sin(y+r-/I)

(4) equating coefficients -B(d,

-u$)

(7) of cos or f I-i(A”tB’+C’)+$ (

sin/I >

(8)

= 0, of sin or t sin /?-t~Bq

l-$A2+B2+C2)+;

>

cosfi

-~EABCo,cos(y+a-~)=O. (5) equating coefficients

(9)

of cos w2 t

( -i:$Aw2sin(ZB-1)-Fcos9=0,

sini

I-k(A2+B2+C2)+$

C(Wf,-U:)COSy+ECCO2

) (10)

MODULATION

ANALYSIS

FOR

BIOLOGICAL

241

MODELLING

(6) equating coefficients of sin o2 t - cbd,

-o:)siny+cCWZ

(

l-$Az+B”+C2)+_1

C2

~0s; >

2

-sqAo12

cos

(2/I-2)fFsin

d, = 0.

(II)

Equations (6)-( 11) form a set of six non-linear algebraic equations whose solution gives directly the amplitudes and phases of spectral components under modulation conditions. In these equations there are three amplitudes, four phases and one frequency which are unknown if forcing at a fixed amplitude F is considered. To give a unique solution two of the phases are eliminated by setting ,!?= 7 = 0”. This is feasible since it merely sets a particular time origin when the self and injected components in the oscillator are in-phase. The resulting equations have been solved using a hill-climbing technique due to Rosenbrock, which has been used successfully in other harmonic balance studies on coupled oscillators (Linkens, 1974. 1976). 4. Results The harmonic balance equations were solved for two types of parameter scan. Firstly, for fixed detuning [i.e. difference between the intrinsic frequency parameter o1 1 and forcing frequency o2 in (1 1)] the magnitude of the forcing function was incremented between solutions. Secondly, for fixed forcing amplitude the frequency detuning was scanned by varying the forcing frequency o2 for fixed w1 1. By using small increments in the parameter the starting values for the variables in the hill-climb were well approximated by the last values in the previous solution. Using this technique convergence to the solution was rapid, and no problems were encountered with convergence to a false minimum. A number of experimental results were also obtained using simulation on an Applied Dynamics AD4 analogue computer. A van der Pol oscillator with small non-linearity (E = 0.1) was used as the forcing function to a second van der Pol oscillator. The output from the forced oscillator was measured on a real-time spectrum analyser Systems Dynamics type SD335. This produces a 1024 point FFT with readout of spectral component amplitudes graphically. The simulations were run 1000 times faster than real-time in equation (1) for (tiI1 = 1 rad s-l, giving oscillations of about 160 Hz. This was a convenient frequency for the spectral analyser and gave fast experimentation. The frequency of simulated oscillations was important since to study modulation phenomena hundreds of cycles of oscillation were required to perform accurate spectral analysis.

242

D.

A.

LINKENS

For a detuning of 5 “(, between the intrinsic and forcing frequencies and E = 0.1, Fig. 4 shows the harmonic balance results for varying forcing function amplitude. When F = 0.16 the self component had dropped to zero indicating that complete entrainment had occurred. As forcing increased the selfcomponent frequency decreased from its intrinisic value of 1 .O showing that for modulation studies it cannot be assumed constant as in combination tone analysis. The injected and reflected sidebands increased with forcing strength almost linearly up to F = 0.1. Beyond that, strong modulation existed and the forcing frequency component increased rapidly and the reflected sideband fell. The reflected sideband was lower than the injected component showing that the behaviour is a mixture between modulation and beating phenomena. The phase angle of the reflected sideband was about + 130” which indicates a tendency towards frequency rather than amplitude modulation (for equal sideband amplitudes, 0” sideband angle indicates pure amplitude modulation _._._.

-

-.-.

_

-.-.-

-.a

0.--.______ 0 L’ ; -----____ --..._ ..*.8‘. ..\* \\. o 1,”

O-

V

OOI

0.00 0.12 0.16 F FIG. 4. Harmonic balance results for E = 0.1, 5Y0 detuning versus forcing amplitude. The continuous lines represent the analytical harmonic balance results, while the blocked symbols are corresponding experimental values from simulation studies.

MODULATION

ANALYSIS

FOR

BIOLOGICAL

MODELLING

243

and 180” indicates almost pure frequency modulation with minimal amplitude modulation). The simulated results compare well with analytical results except for the self component amplitude B. Discrepancies are likely in amplitude estimates because of sidelobe leakage (Durrani & Nightingale, 1972) with the poor frequency resolution of finite length spectral analysis. For F = 0.095 the simulated data and the FFT spectrum are shown in Fig. 5 and considerable amplitude modulation can be seen. and further sidebands are evident since the system was not far from complete entrainment. For F = 0.04 the degree of modulation was simply calculated from the harmonic balance results giving an amplitude modulation index of IO’?,,, an instantaneous frequency deviation of 2.10,, and an FM index (frequency deviation/modulating frequency) of 42 oo. Equivalent results for E = 0.1 and detuning of 15 ‘(, are shown in Fig. 6. Entrainment occurred at F = 0.4 for both analytical and experimental results. The fluctuations in amplitude can be seen in Fig. 7(a) which represents forcing of F = 0.1. The shorter period between fluctuation peaks indicates wider sideband separation as shown in the FFT of Fig. 7(b). In this case no additional sidebands were evident which illustrates the validity of the assumed solution of equation (2). For this detuning large fluctuations occurred as evidenced by an AM index of 36.5 9, when F = 0.2 together with a frequency deviation of 10 “*a and an FM index of 747;. It should be noted that for FM indices about this value more than two sidebands have significant amplitude for sinusoidal FM (van der Pol, 1930). With about 100/o detuning and E = 1.0 the harmonic balance results are shown in Fig. 8. The experimental and analytical results show the same trends although agreement is not so good since harmonics cannot really be ignored for these values of the non-linearity. Entrainment occurred at F = 0.35 experimentally and was predicted from the harmonic balance results at the same value from the poor hill-value of 10e5 (cf. 10-l 5 at F = 0.3) and the obviously invalid results beyond F = O-35. The sideband amplitudes were almost equal and the sideband angle approached 180”. This indicates very small amplitude fluctuations and this is verified by Fig. 9(a). Instantaneous frequency fluctuations are harder to detect visually, but the basic three component spectrum for frequency modulation indices less than 50 9, (van der Pol, 1930) is clearly shown in Fig. 9(b). For F = 0.25, which is near entrainment. the harmonic balance components gave an AM index of only 4 “,,. whereas the frequency deviation was 76 ‘,, which is much larger than for the nearly sinusoidal waveforms of c = 0.1. The corresponding FM index was 909, which suggests that further sidebands should be considered for this strength of forcing and is a further

244

I).

A.

l.INKENS

E 8 -1 T

2 -8

FIG. 5. Simulated results for E = 0.1. 5”, detuning fluctuations. (a) output waveform, (b) spectrum.

0

and F = 09!95

showing

large

amplitude

MODULATION

AKALYSIS

FOR

BIOLOGICAL

245

MODELLJNG

!.O-

n e $ z a

.o-

(

o-

0

04

0.2 F

0.3

0.

FIG. 6. Harmonic balance results for E = 0.1. 15 O,,detuning, versus forcing amplitude. The continuous lines represent the analytical harmonic balance results. while the blocked symbols are corresponding experimental values from simulation studies.

reason for discrepancies between analytical and experimental results for F. = 1.0. The analytical results for a parameter scan of forcing frequency from 0.8 to I .O and F.= 0.1 is shown in Fig. 10. The large effect of forcing amplitude on entrainment frequency is seen in Fig. 10(a), giving values of 6 O0detuning for F = 0.2, and almost zero for F = 0.04. The amount of frequency “pulling” from the intrinsic value depended on forcing strength as shown in Fig. 10(b). The sideband angle approached 180” for small detuning, giving tendency towards frequency modulation. The forcing function phase angle C$approached 90” for small detuning, but had a low value of about 10” for large detuning. The equivalent results for c = l-0 are shown in Fig. 11. For large detuning the method gave good accuracy but towards entrainment the very strong frequency modulation gave additional sidebands which invalidate the assumed solution of equation (2). The amplitude curves of Fig. 1 l(a) should show a monotonic shape which suggests that the method is valid up to frequencies of O-98, 0.96 and 0.92 for F of 0.04, 0.1 and 0.2 respectively. The

-0

FIG. 7. Simulated spectrum.

results

force

6: = 0.1. 15 “(, detuning.

and F = 0.1. (a) output

waveform,

(b)

MODULATION

ANALYSIS _

_

_,----.-

FOR BIOLOGICAL -.--a

-._

MODELLING -

-.-.-

247

_._,_

I--/

FIG. 8. Harmonic balance results for E = 1.0. IO”, detuning, versus forcing amplitude. The continuous lines represent the analytical harmonic balance results, while the blocked symbols are corresponding experimental values from simulation studies. The notation is: o- ~~ ---L&o-ww,, rJ--C.V---A.

phase angles showed the same trends as for E = 0.1. while the amount of frequency pulling was greater. 5. Conclusions

It has been shown that the relative-entrainment phenomenon encountered in forced biological systems is amenable to analysis via the harmonic balance method. The assumption of a three spectral component solution enables output fluctuations due to mixtures of beating, amplitude and frequency modulation effects to be analysed. Analysis is possible because of the way in which sidebands are reflected symmetrically, giving tractable algebraic equations whose solution gives information directly about spectral components.

0

FIG. 9. Simulated soectrum.

results

for /. _ 14.

IO”,,

detuning.

and F=

0.1. (a) output

waveform,

(b)

MODULATION

ANALYSIS

FOR

BIOLOGICAL

MODELLING

249

(0)

0.9

C

I-O w2

FE. 10. Harmonic balance frequency and phases.

results for E = 0.1 and detuning

parameter

scan, (a) amplitudes.

(b)

For non-linear parameter E = 0.1 in van der Pal’s equation with normalized frequency of o = 1 .Othe method has shown that the complex fluctuations due to beating and modulation effects can be simply described by the two sideband amplitudes and phases. For E = 1.O the oscillation waveform shows considerable distortion [e.g. see Fig. 9(a)] but has strong amplitude stability. This is shown clearly under relative entrainment conditions in that frequency modulation predominates over beating and amplitude modulation, as evidenced simply from the harmonic balance results by equal sideband amplitudes and sideband angle approaching 180”. In the heart-rate variability recordings of Figs l-3 the oscillation clearly is non-sinusoidal and the frequency modulation can be detected visually. It should be noted that small instantaneous frequency fluctuations are not readily seen visually, but give significant sideband amplitudes in FFT analysis and in the harmonic balance approach.

I).

A,

LINKENS 3.C

Cl)

_-- ___-______ -.

__\ ‘.\

F=O-04 ----- F=O., - - fzO.2

G B $ 2 a

---__ B 1. --

2.c

----- -

WI

F=0.04

--.----------

,r=o., F=O.2

-

A._

___________

;--

#

-.-.-_-

o-

c / ..--

I-(

C

.__-_______

0 0.8

FIG. 11. Harmonic frequency

0.9

balance and phases.

I.0

,

results for E = 1 .O and detuning

1.8 parameter

--

0.9

scan. (a) amplitudes,

IX (b)

Although harmonics should not be ignored when E = 1 a0it is encouraging that the analysis here has shown correct trends and entrainment limits when compared with simulation experiments. This is important since many biological systems are non-linear and show stable amplitude non-sinusoidal oscillations for which the classic methods (van der Pol, 1934; Kitney, 1975) are not valid under relative-entrainment conditions. The analysis is thus relevant in the study of forced circadian rhythm models such as those of Wever (1972) where E in the range of 0.1 to 10.0 have been considered. Apart from the unidirectional forcing of biological systems such as circadian rhythms, artificial respiration and heart-rate variability mentioned here, the techniques of this paper can be applied to bidirectionahy coupled oscillator systems. Such structures have been suggested as populations of oscillators for circadian modelling (Winfree, 1967 ; Pavlidis, 1969) and as one

MODULATION

ANALYSIS

FOR

BIOLOGICAL

251

MODELLING

or two-dimensional arrays for digestive tract electrical activity (RobertsonDunn & Linkens, 1974; Sarna et al., 1972). It appears that the occurrence of relative entrainment is a common biological phenomenon, and it is of interest to note the large amount of information contained in the spectral amplitudes and phases under these conditions. Thus, beating behaviour can be detected via non-equal sideband amplitudes. Frequency rather than amplitude modulation can only be detected via the difference in sideband phases for which the harmonic balance method is particularly well suited. It is suggested therefore that phase information is of considerable importance in biological systems and should be emphasised more than is done in the majority of engineering applications of spectral analysis, as evidenced by the fact that the majority of real-time analysers of the SD335 type give amplitude measurements only. The data on heart-rate

variability

from which Figs l-3 were obtained

was made

available by 0. Rompelman of Delft University of Technology, Netherlands. and R. I. Kitney of Imperial College and Chelsea College. London. REFERENCES APPLETON. E. V. & VAN DER POL. B. (1922). Phil. Ma,q. 43. 177. AX-HOFF. J.. (1965). Circadian Clocks (J. Aschoff. ed.), p 95. Amsterdam: North-Holland co.

I-‘&l.

DURRANI, FINCHAM.

T. S. & NIGHTINGALE, J. (1972). Proc. IEE, Mar. 343. W. & LIASSIDES. C. (1978). The role of the van der Pol oscillator in the control of breathing in the human, UKSC Conf. on Compu/er Simulurion. Chester. April. HYNDMAN, B. W. & MOHN. R. K. (1973). 10th Int. Con/. @‘Medical and Biological Engineering, Dresden, 223. KITNEY, R. I. (1975). J. theor. Biol. 52. 231. KITNEY. R. 1. & ROMPELMAN. 0. (1978). 8th Int. (‘orzf: on Recent Advances in Biomedical Engineering, Sheffield, April. LAWDEN, D. F. (1959). Mathematics of Engineering Swtems (linear and non-linear), 2nd edn.

London : Methuen. LINKENS, D. A. (1974). IEEE Trans. Cct.Th. 21, 294. LINKENS, D. A. (1976). IEEE Trans. Ccr and SW. 23. 113. PAVLIDIS, T. (1969). J. theor. Biol. 22. 418. ROBERTSON-DUNN, B. & LINKENS, D. A. (1974). Med. & Biol. Eng. 12, 750. ROMPELMAN. 0.. COENEN, A. J. R. M. & KITNE~. R. I. (1977). Med. Biol. Eng. Comput.

15,233.-

239. SARNA, S. K.. DANIEL, E. E. & KINGMA, Y. J. (1972). .4m. J. Dig. Dis. 17. 299. SAYERS, B. McA., (1973). Ergonomics 16, 17. SWADE. R. H. & PITTENDRIGH, C. S. (1967). Am. .Va/. 101. 431. VAN DER POL, B. (1927). Phil. Mag. 3. 65. VAN DER POL, B. (1930). Proc. Inst. Radio Engrs.. N.Y. 18, 1194. VAN DER POL, B. (1934). Proc. Inst. Radio. Engrs., N.Y. 22, 1051. WEVER. R. (1965). In Circadian Clocks(J. Aschoff. ed.). p. 47. Amsterdam: North-Holland

co. WEVER, R. (1972). J. theor. Biol. 36. 119. WINFREE, A. T. (1967). J. theor. Biol. 16.

15.

Publ