A feedback model for biological rhythms

J. theor. Biol. (1972) 36, 153-174

A Feedback Model for Biological Rhythms I. Mathematical Description and Basic Properties of the Model A. JOHNSON AND H. G. KARLSON Department of Electrical Measurements, Lund Institute of Technology, PO Box 725, S-220 07 Lund, Sweden (Received 29 November 1971) To describe sustained biological rhythms-with special interest focused on 24-hour rhythms-a feedback model in block form is presented. A simple linear model containing a time delay serves as a starting point for a more complex and refined model. The ‘oscillating variable in the model, called c(c), is assumed to reflect the measureable biological rhythm. c(r) oscillates around a reference value c&r) and the difference cr.&) -c(t) constitutes an “error signal”. The refined model has three essential features: (i) the error signal is transformed in a non-linear block; (ii) the transformed error signal is delayed and smoothed; (iii) the transformed and delayed error signal affects the oscillating variable itself through an integrator. The introduction of external perturbations into the model is discussed. It is shown that both phase and amplitude of the oscillating variable is affected by perturbations. Methods to test the model are outlined.

1. rntroduction During

recent years much interest within the life sciences has been focussed phenomena. Theoretical and experimental work in the biophysical-biochemical field has led to a relatively good understanding and description of some biologically interesting oscillatory systems (see examples in Pye & Chance, 1972, and Hess & Boiteaux, 1971). At the present stage other such systems seem to be so complicated that a theoretical assessment of these must be rather vague. For example, the number of reactions involved in the observed oscillator is often so large that a point by point description of the whole oscillatory system is out of the question and remains for the future. For the reasons mentioned, biological oscillations will be dealt with here from a fairly general point of view. This also imphes that the results can be on oscillatory

153

154

A.

JOHNSSON

AND

H.

G.

KARLSSON

applied to both slow and rapid oscillations after a simple and appropriate transformation of the time scale. The class of biological oscillations which is called circadian rhythms (daily rhythms) has attracted widespread interest. Nowadays, these circadian rhythms are almost exclusively looked upon as reflecting the self-sustained oscillations of one or more basic oscillators. It is possible to perturb such oscillators by changing the illumination, the temperature conditions, etc. Excellent reviews in the field of circadian rhythms have been given by Sweeney (1969), Btimting (1967), Aschoff (1965) and in the Proceedings of Cold Spring Harbor Symposium on Quantitative Biology (1960). In the case of circadian rhythms their frequency is remarkably constant at different temperatures. This is a necessary condition for the use of these rhythms as time-measuring instruments by the organisms. The nature of the temperature-compensating mechanism in the circadian rhythms is fairly unknown, even if interesting suggestions have been published (e.g. Pavlidis & Kauzmann, 1969). Therefore, temperature effects are not yet incorporated into the feedback model to be presented here. Several models for the circadian rhythms have been published. Some important contributions are those of Btinning (1960), Pavlidis (1967a,b, 1968), Pittendrigh (1960, 1966), Pittendrigh & Bruce (1957), Wever (1964, 1965a,b) and Winfree (1970a, 1972). Each model emphasizes some essential features of the circadian rhythms. The diversity of models might to some extent be explained by the diversity of organisms studied. Generally speaking, Btinning favours a model where the basic oscillations can be regarded as a type of relaxation oscillators (Btinning, 1960). Btinning’s model accounts for a feature, observed in many biological rhythms, viz. one light “sensitive” and one light “insensitive” phase of its cycle. Pittendrigh’s model (Pittendrigh, 1960, 1966) also provides a context to describe this feature. His basic oscillator is measurable via a second coupled oscillator. The basic oscillator is phaseshifted “immediately” by a suitable light pulse but isinother respectsunaffected. Other models mentioned use different mathematical approaches of studying non-linear oscillators. Wever (1964, 1965a,b) starts with an extended van der Pol equation, which explains several features of circadian rhythms though it means a restriction to a rather special class of oscillators (cf. Pavlidis, 1967a). Kalmus & Wigglesworth (1960) and Pavlidis (1967a,b, 1968) give a topological description of the behaviour of the oscillator. Winfree (1970a,b, 1972) has the same general approach and focusses his interest on the possible existence of singularities in the phase plane representation of the eclosion rhythm in Drosophila pseudoobscura. The topological way of treating a non-linear oscillator seems particularly well suited for discussions about such singularities (see, e.g., Winfree, 1972).

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

I

155

Some of these models show evident drawbacks. For example, Pittendrigh’s model does not explicitly account for amplitude variations and both Pittendrigh and Pavlidis have used a second ascillator to explain the time required to reach steady state after phase shifts observed in some organisms. The present feedback model includes these features as properties of one single oscillator. The presentation of the feedback model considered here is made in two parts. The first part deals with the general outline of the model. The second one, called part II (Karlsson & Johnsson, 1972), contains applications of the tenets of the model to some specific cases and organisms and discusses phase and amplitude problems. The content of the present papers fully covers the rough outline of the model given a couple of years ago (Johnsson & Karlsson, 1969; Karlsson & Johnsson, 1969). 2. Abbreviations and NoOatiom The following symbols will be used throughout

this paper.

oscillating variable of the basic oscillator variable representing the value around which c(t) is oscillating variable representing a perturbation to the basic oscillator function of c(t) and c,,,(t) whose characteristics will be discussed later function of time (will also be referred to as a weighting function) notation for an integrator a constant parameter denoting a time delay period time and angular frequency of the oscillator constant parameter phase shift of the oscillations (“ + ” means advance, “- ” means delay) 3. Basic Consideratiuns

This paper deals mainly with the basic oscillator-the “central clock” in the sense of Btinning (Biinning, 1967). This oscinator is revealed and studied through its timing of different “output” reactions, for example, the eclosion of Drosophila pupae. It can also be a&ted through different receptor systems with, for example, a light pulse as the primary event. It should be emphasized that very little is known about the receptor systems. It is reasonable to assume that a light “pulse”, characterized by its intensity, its length, and its infinite slopes, will not keep this shape until entering the basic oscillator. Instead, for example, the pulse front is likely

156

A.

JOHNSSON

AND

H.

G.

KARLSSON

to be smoothed out by different kinetic steps. These steps are very complicated, however, and the details are unknown. They should, therefore, be avoided as long as possible in a model which attempts to describe the qualities of a basic oscillator. The following points might serve as guide lines when modelling a basic oscillator, a clock: (i) the oscillator in the model should be as simple as possible and yet (ii) explain as many characteristic features in the rhythms as possible, and (iii) the model should visualize as clearly as possible the components of the oscillator. The points mentioned imply that a minimum number of parameters and variables should be used and that the receptor systems and output mechanisms should be avoided, if possible, to explain basic characteristics of the oscillating system. The present model starts with some simple concepts which are used in control systems where time delays play a role. The model has been developed from a previously published idea, that in circumnutations of some plants (oscillations of the stem around the plumb line) the effective time delay between the stimulation and the reaction causes the oscillations (Israelsson & Johnsson, 1967; Johnsson & Israelsson, 1968, 1969; Johnsson, 1971). 4. Feedback Model of an Oscillator (A)

A SIMPLIFIED

MODEL

In Fig. 1 the principle of a simple linear model for the basic oscillator is given. The variable is denoted c(t), indicating that we probably can find a chemical substance varying with time which is involved in the biological oscillator. The concentration of this substance is called c(t). The oscillations

c,,,(tl-c(t) h kx[c,&)-c(t)]

C,dIfi

\

-+

kx[c,,,(t-t,)-c(t-t,il I

‘-

c(t)

:

i

I

1. Block diagram of a linear feedback model generating oscillations in c(t) with suitable values of the multiplier k and the time delay to. The value of the reference c&t) is determined outside the model and might vary with time. FIG.

of c(t) occur around a reference value which will be denoted c,,~. Sometimes c,,r is found to be time-variable, for example, because of a variable average light intensity. Then it will be called c,,,(t).

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

I

157

Suppose c(t) is observed directly or indirectly, dependent upon the biological material and the measuring method. The actual value of c(t) is in the basic oscillator compared with c,&, and the difference, constituting an “error” signal, is fed into the first one of the thee blocks depicted in Fig. 1. The arrows indicate the direction of signal flow. Each block in Fig. 1 represents a transformation of the signal entering that block with the transformed signal appearing at the output. Let us now, in this simple model, assume that the signal is transformed as follows: the input signal, i.e. c&)-c(t), is multiplied with a constant k; second block (to): the input signal, k x [c&)-c(t)] is delayed t,, t.u. (time units). The output signal from the block is then k x CcrcrO- fo) - CO- &,>I ; third block (I): the input signal, k x [c&t - to) - c(t - t,)] is integrated. Thus, the time derivative of the output signal, c(t), equals the input signal to this block. first block (k):

In section 5 these blocks will be further developed and their relation to biological conditions will be discussed. With suitable assumptions concerning k and t,, the system will oscillate with the frequency being dependent upon these parameters. This is studied in the following section. (B)

THE

SIMPLE

MODEL

IN

MATHEMATICAL

TERMS

To get a convenient description in mathematical form of the control system in section 4(~) it is suitable to look at the signal just before block I in Fig. 1. At this point the signal can be described as dc(t)/dt, since the integration in the block I should give c(t). On the other hand, the discussion in section 4(~) showed that the signal at this point equals k x [c,,& - to) c(t - t,)]. Thus, the following equation must hold :

W) -

= k x [c,,& - to) - c(t - t,)].

dt For simplicity, assume c,,&) to be a constant equal to zero (cf. the mathematical treatment in Appendix A). Then equation (1) can be written as -WO = -kx

dt

c(t-t,).

This equation has the solution (see Appendix A) c(t) = A(t)x cos (wt+rp), (3) where A(t) is the amplitude, o the angular frequency, and cp a phase angle.

158

A.

JOHNSSON

AND

H.

G.

KARLSSON

A study discloses that the amplitude A(t) is constant if k = n/2t,. If will increase, if k < n/2t0 it will decrease. The angular frequency o of this oscillation can be calculated and the following relation will hold when the amplitude is constant k > x/2t, the amplitude

w = 27+&,

or

t = 4t,.

Equation (l), which corresponds to the model in Fig. 1, can be used for numerical iterative calculations of the behaviour of this oscillator. Figure 2 shows such results. 50 40 T 30 0 20 I@ 0

0

100

200

300 400 500 t (timeunlls)

600

700

FIG. 2. Oscillations produced by the model in Fig. 1. In the curves A, B, and C the multiplier k is less than, greater than, and equal to a critical value, respectively. It is only in the critical case (curve C) that the amplitude is constant in time. - - - -, Curve A; -- , curve B; -, curve C.

Figure 2, curve A: k < n/2t,,. to has been chosen as 34 t.u. (time units) and k = O-0415 (t.u.)- *. The oscillator is seen to be free-running but with decreasing amplitude. The period length z is constant and equals 140 t.u. Figure 2, curve B: k > n/2&,. to as in curve A, k = 0.0488 (t.u.)-‘. The amplitude of the oscillations is now increasing. z is constant and equals 134 t.u. Figure 2, curve C: k = n/2t, = 0.0462 (t.u.)-‘. Under these circumstances the amplitude is constant. This is the “critical” case and a slight disturbance of the parameters will change the behaviour of this oscillation to either curve A or curve B. z = 136 t.u. Because of this sensitivity to slight changes in the parameters (e.g. k) the simple, linear model is rejected as a suitable description of a physiological oscillator. A non-linear model will therefore be presented below.

FEEDBACK

MODBL

FOR BIOLOGICAL

RHYTHMS,

I

159

5. Extemded Feedback Model for the’ Basic Oscillator Figure 3 shows an extension of the control scheme in Fig. 1. The new blocks in Fig. 3 introduce a non-linear function,f, into the model as well as a more realistic delaying block, H. The features off and H are discussed

FIG. 3. Block diagram showing a modified version of the feedback model in Fig. 1. The multiplier k has been replaced by a non-linear function f (see section 5(A) and Fig. 4) and the pure time delay to by a more realistic delaying unit in block H (see section 5(~) and Fig. 5). The numerals I-IV denote different positions where perturbations can be introduced,

and motivated below as well as the presence of the integrator in block I. Readers not specially interested in these partly mathematical details, are recommended to continue with section S(D), dealing with perturbations given to the oscillator. (A) BLOCK

f

This block contains a more realistic function f instead of the simple constant k in the example above. The relation between the input signal and the output signal from this block is described by J The input signal is q&t)- c(t) as seen in Fig. 3. The block f might, for example, transform this signal into some other chemical substance via one or more chemical reactions. In a physiological case, it is reasonable that the output signal of the block is a non-linear function of the concentration difference c&t)-c(t). In this paper the function f is often described by the coordinates of some characteristic points as shown in Fig. 4, curve A. f can be either symmetrical or asymmetrical. Symmetrical means that f(c) = -f( - c). Both cases will be discussed in this paper. From the preceding section 4(~) it is concluded that the slope off near the origin [where c(t) = c,,(t)] must not be less than a certain critical value if sustained oscillations should start spontaneously. Figure 4, curve B

160

A.

JOHNSSON

AND

H.

G.

KARLSSON

shows one example of f where the slope is too small oscillations.

-04

I M’ ,-

- 0 2,/’ 1111111111. -I_o- -.:-._-------,

2’,4

,‘/ / -. //’ IIIIIIII~ 5

4

/--- ’ ,’

to produce such

A

crttm slope

,----Et IO Cre, -c

-

j-02

FIG. 4. Some examples of the function f. The coordinates of the points 1, l’, 2, and 2’ shown to the left of curve A, are used to characterize the form of $ For a certain block H (Fig. 3) curve A will give self-sustained oscillations while the slope of curve B near the origin is so small that only damped oscillations occur. Curve C to the right contains a “dead zone” around the origin. Oscillations with small amplitude will then fade out but self-sustained oscillations will occur if the amplitude, through environmental perturbations, becomes large enough.

Figure 4, curve C gives an example off with a “dead zone”, i.e. the slope equals zero around the origin. Outside the “dead zone” the slope is rather high. This example can be said to represent a combination of curve A and curve B. In this case oscillations will not start spontaneously. A perturbation strong enough, however, will immediately give rise to oscillations with relatively large amplitude. The mean slope off, influencing these oscillations, will be sufficient to prevent a decaying amplitude and the result will be sustained oscillations. This special f will also be referred to in part II, section 4(~),t where a possibility of stopping a rhythm is discussed. (B)

DELAY

BLOCK H

In the simple model in Fig. 1 above, the block just delayed the signal t, time units. However, in a physiological system it is likely that a signal is both delayed and transformed when it appears at the output of this block. We can, for example, imagine this to be the case when fluid transport is involved in a system. An abrupt concentration step will be detected at another position in such a system as a delayed and a slowly increasing t Throughout

this paper part II refers to Karlsson & Johnsson, 1972,

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

161

I

concentration profile. Thus, the step has been both delayed and transformed. The function which describes the transformation of the input signal is called h(t’). This function is defined by the notations in Fig. 5. An exponential decay is very frequent in physical and chemical systems and is therefore a reasonable assumption. The rising part, tihe front, of h(t’) is assumed to be linear mainly for numerical reasons. Further, h(t’) is normalized, i.e. its time integral is unity. The input signal to block H, i.e. f[c,,,(t)-c(t)],

0 /(time

units

1

FIG. 5. An example of the function /I in block H. The points tI and la are the limits of the linear part of h(t’) and A determines the rate of decay of the exponential part. The slope of the linear part is always chosen in such a way that the integral of h(t’) equals unity. Thus, the parameters tl, rz, and /I fully characterize the function h(t’).

gives an output signal mathematically

described at the time t by

~~(~‘)./[c,Xt-t’)-e(t-t’)]

dt’.

(4)

Another (equivalent) way of looking upon the function h(t ‘) is obvious from equation (4). The argument for the concentration is (t- t ‘) which means that the integral (4) depends upon concentrations not only at time t but also at times before (t’ > 0). The concentration at a certain time t--t’ is weighted according to the value of h(t’) by integrating the product of h(t’) and f[c,,,(t-t’)--c(t-t’)]. Th eref ore h(t’) is usually called a weighting function. In a physiological system this corresponds to a situation where the response at time t is built up as a “sum” of stimulations prior to time t. Since old stimulations contribute to the response at time t through the action of the weighting function, this function can be interpreted biologically as a type of “memory”. Very old stimulations are likely to give small contributions to the total “sum” [integral (4)]. The weighting function is therefore chosen to be exponentially decaying to very small values for old stimulations as mentioned above. Theoretically, it is possible to determine the shape of h(t ‘) by giving a 11 T.B.

162

A. JOHNSSON

AND

H.

G.

KARLSSON

very short and intense perturbation to the feedback system itself. It can be shown that the system response as a function of time will then be proportional to h(t’). The technique is, however, often impossible to use in practice. A short and intense light pulse given to the organism studied might be smoothed by the light receptor system. The function h(t’) thus contains two important features of the model: (i) it determines how soon after a perturbation all blocks of the oscillator in Fig. 3 have been affected. The time necessary is due to the time delay involved; (ii) it describes how soon after a perturbation all the resulting transients in the oscillator have faded out. This time is due to the exponential decay of h(C). Whenfis symmetrical [cf. section 5(A)], the period time for the oscillations is roughly given by the “centre of mass” of the function h. Some results showing this are presented in Table 1. The period time 7 of the basic oscillator is tabulated for certain values off and h. The influence of h on z is pronounced and bigger than that off (cf. Appendix A). TABLE

1

The inJEuenceof the functions f(c,,f -c) and h(t’) on the free-running period t of the oscillator in Fig. 3. f(+ -c) was chosen to be symmetrical, i.e. f(x) = - f( - x). As is seenz is determinedmainly by the centre of massof h(t’) Non-linear f(Gd Notations (point 1) = - (point I’)

function - cl in Fig. 4 (point 2) = - (point 2’)

(6, 0.39) (6, 0.39) (6. 0.39)

(16,0.48) (16, 0.48) (16. 0.48) (20; o.soj (16,057)

cio;0.50, (2, 0.20)

Weighting function Centre of mass of h0’) Notations in Fig. 5 W’) tz t1 t.u. t.u. CU. (t.u!)- 1 ----30 34 0.3 35.58 10 40 0.3 32.42 10 40 0.03 59.89 30 34 0.3 35.58 30 34 0.3 35.58 (C) THE INTEGRATOR

Freerunning period r

t/(centre of mass)

t.u. .- ~--~ 142 128 213 141 143

3.99 3.95 3.56 3,96 4.02

1

The error signal, being transformed and weighted in the blocks f and H respectively, eventually affects the output signal c(t) through a third block. This block has in the present model been chosen as an integrator. Thereby, c(t) will be adjusted until the integral (the time average) of the transformed error signal is zero. Thus, in the long run, the mean level of c(t) cannot deviate from c,,,(t) through the influence of the environment but will be corrected after some time. In Appendix B the alternatives for the third

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

I

163

block are discussed in greater detail. The integrator is then found to have properties which for the present application are preferred compared with other alternatives. (D)

INTRODUCING

PERTURBATIONS

INTi3

THE

OSCILLATOR

Many experiments on physiological oscillators deal with the effects of perturbations, often light pulses. Giving a light pulse to the receptor system will eventually lead to a chemical signal which enters the oscillator. This signal might have a shape different from that of the light pulse, as mentioned before. A chemical perturbation which enters the signal path in the scheme will be denoted cP, where p stands for perturbation. In the model it is assumed that the blocks themselves are unaffected by cP. At the present stage it is possible to introduce c,, in different positions of the scheme in Fig. 3. These are denoted by the numerals I-IV. I. The c,-signal is transformed by the action of the block f before being delayed in the block H. The change in the output c(t) can be detected after a certain delay. II. cP is not transformed by f but delayed and integrated before it appears at the output. III. In this case c, will cause an immediate change of c(t). (Note, all the same, that h(C) in block H means a delay in the system, causing it to oscillate.) IV. ci, causes an immediate change in the output. (E) C,,,(t) The reference value, crcf, might of course be changing with time. It is reasonable that oscillations in constant dim light will occur around a new reference value if the light conditions are altered. Such an interpretation is discussed by Karvt, Engelmann 8z Schoser (1961) in the case of KuZuncho& blossfeldiana. Thus, a light treatment can have a dual effect: it can change the reference level, crcf, around which the oscillations occur, and it can enter and perturb the oscillating system in position I-IV of the model. The receptor systems for these two effects need not be the same and the absorption spectra can even be different (see KarvC, et al., 1961). (F)

OUTPUT

LIMITS

If c(t) is to describe a concentration, it has to be positive and restricted in amplitude. It will therefore be assumed that c(t) > 0. If, for example, perturbations are introduced in the system tending to send c(t) to negative values, this lower limit of c will prevent negative output concentrations

164

A.

JOHNSSON

AND

H.

G.

KARLSSON

occurring. In much the same way, c(t) is likely to be restricted by an upper limit. It might be worth stressing some implications of the limitations of the oscillations. For example, interesting features of the oscillator may not turn up if perturbations given to the oscillator are too big. Such perturbations may simply force the oscillator strongly towards some limit, thereby preventing features of the normal, unperturbed oscillator to be observed in such a case. This gives a hint that small perturbations normally should be preferred, if possible for practical reasons, as they do not force the system out of its normal working range. 6. Outline of Methods to Test the Model Two variables of an oscillator are normally studied, viz. the phase and the amplitude. Observation of phase changes due to single perturbations is the method most frequently used to study biological clocks. Amplitude studies are less frequent, since the oscillator is often inaccessible to direct measurements. Phase and amplitude changes of the oscillator are experimentally caused by, for example, light pulses which are introduced into the oscillator via a receptor system. The experiments should preferably be designed in such a way that the lack of information of the precise function of the receptor system is of minor importance. One such approach is to compare the effects of one or more identical light pulses given at different times in the cycle. Another type of such experiments is qualitative comparisons of the effects of pulses with different signs, for example, light pulses and dark pulses (measured relative to the normal behaviour under dim light). On the other hand, experiments with light pulses of different wavelengths will only help to reveal the nature of the photoreceptor system and will normally say nothing about the basic oscillator. Phase and amplitude studies will be mentioned here from a general point of view. Details as well as comparisons between the model and actual biological experiments will follow in part II. (A)

PHASE

STUDIES

Phase changes of biological oscillators, caused by perturbations, are often represented in phase response curves. As is well known, these curves show the phase shift of the rhythm as a function of the phase (time) of the rhythm when the perturbation was applied. For the present model phase shifts were determined as follows: a correspondence was established between characteristic points (maxima) of the perturbed oscillation and of an unper-

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

I

165

turbed oscillation. Except for some rather specific perturbations (see part II, section d(B)) this could be done unambiguously. By comparing corresponding points it was decided whether the perturbed oscillation was delayed or advanced relative to the control, defining a phwe delay or a phase advance, respectively. Thus, with the exceptions suggested above, a phase delay cannot be taken for a phase advance or vice versa. The magnitude of a phase shift is usually well defined experimentally. On the other hand, the time or the phase of the rhythm when the perturbation is applied is measured relative to some reference time or reference phase, chosen for practical purposes by the experimenter. To illustrate the phase shift concept, a result from the model is shown in Fig. 6. The theoretical time course of c(t) is presented for the parameters given in the legend and a perturbation applied in advance position. Maxima (A) and minima (v) of the unperturbed control are indicated at the top

t( time units )

Fro. 6. A result of the feedback model after a perturbation, shown as a broken bar, in position I in Fig. 3. The coordinates off were: 1 = (2; O*lO), 1’ = (-4;. -0.32), 2 = (5; 0*15), and 2’ = (-7; -0.37). The characteristics of h(r’) were: t1 = 30 t.u., 4 = 34 t.u., and h = 0.3 (t.u.)-‘. The triangles at the top of the figure denote maxima (A) and minima (V) of an unperturbed control oscillation. The effect of this perturbation on the steady-state oscillation was evidently a phase advance.

of the figure. The phase shift caused by the perturbation was in steady state, interpreted as a phase advance. A phase response curve can be constructed from a set of curves with the perturbation given in different phase positions. In Fig. 6 it is also observed that the amplitude of c(t) is affected by the perturbation. A drawback with the usual phase response curves as published in the literature is that no amplitude information is included. Phase response curves usually show the effects of single perturbations applied at different phases. In part II phase shift experiments are discussed where two pulses are given to the circadian rhythm. Experimental results and theoretical predictions for perturbations of different signs are also presented there.

166

A.

JOHNSSON (B)

AND AMPLITUDE

H.

G.

KARLSSON

STUDIES

Amplitude studies of biological oscillators are in principle as important as the phase studies. Though mostly difficult it is in some cases possible to observe amplitude variations directly in the system studied (e.g. Zimmer, 1962). Examples exist where the oscillating system also shows an interesting relation between phase shift and amplitude after a perturbation. Such a relation is predicted by the present model as distinguished from that of Pittendrigh. This difference might be used to discriminate between the two models (see part IT, section 3(c)). The connection in the model between amplitude and phase shift is illustrated in Fig. 6. It is also observed that when the amplitude is small the period length of c(t) is less than that of the control. Therefore, the phase advance in Fig. 6 increases with time until the steady-state amplitude is reached. As the rate of amplitude growth is mainly dependent upon the function f (cf. section 4(~)), this function will affect the theoretical steadystate response curve. As already mentioned, certain perturbations will reduce the amplitude of c(t) in the present model. In fact the amplitude can become even zero or very close to zero after a suitable perturbation given in a suitable phase. This situation will be discussed in more detail in part II, section 4(~), with reference to experiments reported in the literature. 7. Discussion The aim of part I of the present papers is to present a model for biological oscillations. Further tests and comparisons with experiments will be treated in part II. Several theoretical approaches to the subject of biological oscillations are possible. So far at least three different types can be distinguished, all being intimately related. The first one uses a more or less complicated differential equation to describe the basic oscillator. Wever’s application of van der Pal’s equation to circadian rhythms is the most elaborate representative of this approach (Wever, 196%~). The second type of approach is to use a number of ordinary first-order differential equations, which can be represented in a topological way (a state-space representation). For example, the van der Pol equation can be expressed as two coupled first-order differential equations, thus governing two variables. For each time c one variable can be plotted against the other in a two-dimensional diagram. An oscillating process in steady state will form a closed loop in this diagram. If three variables are involved a threedimensional diagram is needed, etc. It will be more and more difficult to

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

I

167

visualize the “plot” as the number of dimensions increases but the principle will be the same. In the case of circadian rhythms this approach was used by Kalmus & Wigglesworth (1960) and in a more elaborate version by Pavlidis (1967qb). These authors stress the generality of the topological approach. Winfree (19704b), in his penetrating Drosophila studies, focusses his interest on an important and, in the topological approach, rather conspicuous property of oscillating systems, viz. a singular state. This point will be dealt with in part II, section 4(B). The third type of approach is the present one using a feedback model and methods from control theory. This approach is very generally applicable and can often be used instead of the first and second type mentioned above. One advantage of a feedback model is that it can give a clear picture of the different biophysical subsystems (blocks) and of the relations between their input signal and output signal (“cause” and “effect”). The present model is by necessity simplified and for each organism further blocks and “blocks within blocks” can be added. One such extension is a c,,rvalue varying with the time, mentioned in section 4(~), which is certainly the case in many biological systems. Other modifications of the present block scheme are also conceivable: as mentioned above, perturbations could be introduced at different positions in the block diagram (Fig. 3), the order of the blocks with regard to each other could be changed, block properties could be affected by light, etc. The significance of one specific feature of the present model, viz. the pure time delay involved, also deserves a short discussion. This feature cannot be represented topologically unless an infinite number of variables are used, while it is easily included in the feedback model. Pure time delays could arise due to, for example, transport processes and are known to play an important role in certain rhythms (Grodins & James, 1963; Israelsson & Johnsson, 1967). In this case the feedback approach with a pure delay has an advantage. The pure time delay should not, however, be emphasized as a general feature of oscillating systems. As the length of the pure time delay is decreased in the feedback model [by a change in the weighting function h(t’)], the similarities with the other two approaches increase. Probably most biological rhythms can be more or less satisfactorily described with aid of any of the three approaches mentioned. Finally a few words should be said about relaxation vs. pendulum oscillations, a subject much discussed in the field of circadian rhythms. The feedback model is adjustable in this respect, i.e. oscillations can be obtained which should be judged as either pendulum or relaxation oscillations, the type being dependent on the properties of $ Two examples are given in Fig. 7, where the f-functions and the corresponding oscillations are shown to-

168

A.

JOHNSSON

AND

H.

G.

KARLSSON

fiC,,&C)

30-

! 34

1 02 -10

-5

5

ILLLLLU

4

/

IO Cre, -c

-02

20

1 -04

T

JAfmA

IO0 L

.LL--L-

I

1

I

I

1

1

I

I

1

;; f&,-c)

T-10

-5

-04 01 -02

5

'0 &t-c

3020 IO

OOY,-??

200

400 tftime

600

umts)

FIG. 7. Two curves from the feedback model demonstrating that both pendulum (upper curve) and relaxation (lower curve) oscillations can be obtained dependent upon the choice of J The f-functions used are in both cases shown to the left of its corresponding oscillation, c(f). The same weighting function h(r’) as in Fig. 6 was used.

gether. It is seen that a linearfwith a suitable slope can give very sinusoidal oscillations, like a swinging pendulum, while a more complicated f can give saw-tooth like oscillations. The properties of f are, of course, likely to be different in different rhythms. Wever (see, e.g., Wever, 196%) has pointed out that pendulum and relaxation oscillations are two extremes on a continuous scale with gradual transitions. This is found to hold also in the present approach. It is a pleasure to thank ProfessorC. H. Hertz of this institute for his encouragement during the work as well as for the reading of the manuscript. Dr W. Engelmann, Tiibingen, Dr A. T. Winfree, Cambridge, and Dr M. K. Chandrashekaran,Tiibingen, have given kind and valuable adviceat the completing of the manuscript. The work has been financially supported by the Swedish National Science ResearchCouncil, which is gratefully acknowledged.

REFERENCES (ed.) (1965). Circadian Clocks; Amsterdam : North Holland Publishing Co.

ASCHOFF,

J.

Aoc.

Frldajtqy

Smvnrr-

School

1964.

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

I

169

WNN~, E. (1960). In Coid Spring Harb. Symp. qaant. Biol. 25,l. BINNING, E. (1967). The PhysioIogicaI Clock. New York: Springer-Verlag. Cold Spring Harbor Symposia on Qaantitative BiorogV (1960). 25. GRODINS, F. S. & JAMES, G. (1963). Ann. N. Y. Acad. Sci. 109,852. HESS, B. & BORDEAUX,A. (1971). A. Rev. Biochem. 40,297. ISRAELSON, D. C JOHNSSON,A. (1967). Physiologia PI. w), 957. JOHNSSON,A. (1971). Physiologia PI. 24,419. JOHNSSON,A. & ISRAELSON, D. (1968). Physiologia Pl. 211,282. JOHNSSON,A. & ISRAELSSON,D. (1969). Physiologia PI. 22,1226. JOHNSSON,A. & KARLSSON, H. G. (1969). In Proc. Conf. of the Swedish National Committee for Physics, Manual B, 90. Stockholm. KALMUS, H. & WI~GLESWORTH, L. A. (1960). In Cold Spring Harb. Symp. qaant. Biol. 25,211. KARLSON, H. G., B~RXSSON, P. 0. & JOHNSSON,A. (1971). Med. % Biol. Engng. 9,721. KARLSSON, H. G. 8r JOHNSSON,A. (1969). In Proc. Conf. of the Swedish National Committee for Physics, Manual B, 105. Stockholm. KARLSSGN, H. G. & JOHNSSON,A. (1972). J. theor. Biol. 36,175. Ku&!, A., ENGELMANN, W. & SCHOSER,G. (1961). P&a 56,700. PAVLIDIS, T. (1967a). Bull. math. Biophys. 29,291. P~v~n>rs, T. (1967b). Bull. math. Biophys. 29,781. PAVLIDB, T. (1968). In Lectures on Mathematics in the Life Sciences, Vol. 1. (M. Gerstenhaber, ed.), p. 88. Providence, Rhode Island: American Mathematical Society. PAVLIDIS. T. & KAUZMANN, W. (1969). Archs Biochem. Biophys. 132,338. P~TTENDRIGH, C. S. (1960). In Cold Spring Harb. Symp. qaant. Bill. 25, 159. P~~~ENDRI~~H,C. S. (1966). Z. PpPhysioI. 54,275. ~NDRIGH, C. S. & BRUCE, V. G. (1957). In Rhythmic and Synthetic Processes in Growth. (R. Rudnick, ed.), p. 75. Princeton, N.J.: Princeton University Press. PYE, E. K. & CHANCE, B. (1972). Biochemical Oscillators; Proc. Johnson Res. Found. Symp., Prague, Czechoslovakia, 1968. New York: Academic Press. SWRENEY, B. M. (1969). Rhythmic Phenomena in Plants. London: Academic Press. Wavaa, R. (1964). Z. Tierphysiol. Tiererniihr. Futtermittt?Ik. 21, 359. Weve~, R. (19650). In Circadian Clocks. (J. Aschoff, ed.), p. 47. Amsterdam: North Holland Publishing Co. WEVER, R. (1965b). In Circadian Clocks. (J. Aschoff, ed.), p. 74. Amsterdam: North Holland Publishing Co. WINFREE, A. T. (197Ou). In Lectures on Mathematics in the Life Sciences, Vol. 2. (M. Gerstenhaber, ed.), p. 110. Providence, Rhode Island: American Mathematical Society. WINFREE, A. T. (19706). J. theor. Biof. 28,327. WINPREE, A. T. (1972). In Biochemical Oscillators; ptoc. Johnron Res. Found. Symp.. Prague, Czechoslovakia, 1968. (E. K. Pye & B. Chance, ed.). New York: Academic Press. ZIMMER, R. (1962). Pfanta 58,283.

Appendix A The properties of the following equation -W)

= k.[c,,&-to)-c(t-to)]

dt will be discussed here for different values of the constant parameters k and t, (k > 0, t, > 0). c&t - to) is assumed to be constant and equal to crCf.

170

A.

JOHNSSON

AND

H.

G.

KARLSSON

As (Al) is a linear equation each frequency component of its total solution can be studied separately and then superimposed whereby the final result is obtained. Therefore, suppose that c(t) = cref+ co cur cos (cot+ cp) 642) satisfies (Al). co, rr, o and cpare constants to be determined. (A2) and (Al) then give

% [cd+

co e”’ cos (cot + rp)] = - kc, eb(‘-*o) cos [o(t - to) + q]

or u cos (wt+rp)-w

sin (wt+cp) = - k e-O*” [cos ot, cos (wt + cp)+ sin ot, sin (ot + cp)], which must hold for all f. Identification of coefficients for cos (ot+ rp) and sin (wt + q) then gives u = -k eearo cos do, (A3) and u = k emat0sin ot,. (A4) The equations (A3) and (A4) give implicitly cr and o for all combinations of k and to. co and cpin (A2) can be determined through the initial conditions. Thus (A2) is an oscillating solution to (Al) with the amplitude A(t) = co cur and the angular frequency CD. The interest will now be focussed upon solutions with constant or almost constant amplitude A(t). A(t) is constant if r~ = 0. Then (A3) and (A4) give wt, = n/2 + nn, and w = k( - l)“, where 12is a positive and even integer since k, to and o must all be positive. Below, n will be given as a subscript of the actual variable, e.g. w,, a, and A,,(t). Thus kt, = 1112+ mr, n = 0, 2,4, etc. (A5) and 0, to = 7112+ nn, n = 0,2, 4, etc. 646) must hold for A,(t) to be constant with time. For n = 0 in (A5) k = n/2t,. (A7) This relation is used in the present paper. Since the period time r = 27r/o equation (A6) in this case can be expressed as z = 4t,. (‘48) To find out what will happen if the condition (A5) is almost but not quite

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

I

satisfied, the equations (A3) and (A4) oan be expanded in Taylor after kt, around n/2 + nn. After some calculations this gives

171 series

and

%to

(42 + nn)

for n = 0, 2, 4, etc. and kt, in a small neighbourhood of n/2+nz From (AS) it follows that the amplitude A,(t) = con e”” increases with time if kt, > z/2 +nx and decreases with time if kt, < z/2 +nlt. This can be proved to hold generally and not only with kt, near n/2+nn. (AIO) gives the dependence of frequency upon kt,. The equation (AS) is illustrated in Fig. 2 of the present paper. Let us now suppose that sustained oscillations are observed in a system governed by (Al). As all possible solutions must then be limited in amplitude, kt, < 7c/2+nn for n = 0, 2, 4, etc., i.e. kt, < ~12. Otherwise the solution A,(r) cos (o. t+ rpo) would force the system into saturation. Further, sustained oscillations require that A,(t) is roughly constant, implying that kt, R n/2. Then the solutions A,(t) cos (CD,?+ rp,) with n = 2, 4, etc. will all be heavily damped. Accordingly they will tiy affect the transients and not the steady-state solution, A,(f) cos (o. t + cpo).

Appendix B In control systems the deviation of the output signal from its reference value (the error) is used to influence the output signal in the direction desired. Three different regulatory mechanisms are frequently used in control theory for that purpose. These will be briefly studied here for the simple model in Fig. 1. The properties of the third block, denoted 1, will then be varied. The output signal of block I will here be called s(t). This signal is a function of k[c,,,(t- to) - c(t- to)]. The system is also assumed to be affected by environmental perturbations, c,(t), at the output of block I not shown in Fig. 1. The corresponding point in Fig. 3 is called position IV. Thus the system output, s(t), is the sum of s(t) and c,(t): c(t) = c&)+s(f). (Bl) In the present analysis c,,(t) and c,,,(t) are supposed to be constant and below will be denoted cP and cref, respectively. It can be shown that the mean level of c(t) under the present assumptions

172

A.

JOHNSSON

AND

H.

G.

KARLSSON

will reach a final value, E, in steady state. Mathematically

Eis defined through

C = lim f ~~‘c(r’) df’, WY t+m where r is the period time for c(t). Z is interesting since the mean level of many biological oscillations is rather insensitive to environmental perturbations. Therefore, the function in block I will be chosen so as to minimize the influence of cP on T. Three different regulatory mechanisms will now be discussed. They are characterized by the way in which s(t) depends on k[c,,r - c(t - t,)]. s(t) = k[c,,,-c(t-&))I. 033) (9 (Bl) and (B3) then give c(t) = c,+k[c,,f-c(t-lt,)]. (B4) Forming the time average as in (B2) we get Z = cp + k[c,,‘-

Z],

or

It can be shown that k must be less than unity in order to avoid an everincreasing amplitude of c(t). Thus, the influence of c, on E is at least 0*5c, s(t) = k;

(ii)

[cref-c(t-t,,)].

w-3

(Bl) and (B6) then give c(t) = c,+kil

[cref-c(t-Q].

(B7)

Time averaging gives c = c,+k$(c,&-kx

lim 1 [c(t+z-t,,)--(t-t,,)]. t-m 7 Since c,,r is constant and c(t) has the period time r c = cp. Thus, 2 does not depend on c,,r at all but only on the perturbations (iii)

s(t) = kj [cref-c(f’-

t,)] dt’.

W9 cP. (B9)

(Bl) and (B9) then give c(t) = cp+ k d [cc

c(t’ - t,)] dt’.

@lo)

FEEDBACK

MODEL

FOR

BIOLOGICAL

RHYTHMS,

173

I

By forming first the time derivatives of (BlO) and then the time averages in steady state we get

$ = k[Crer-C(f-to)], and lim 1_[c(t+5)-c(t)] = k[c,,,-E]. t-rco 7 But c(t) has the period length z giving 0 = k[C,,f4], or c = C,,f. Wl) Thus, I! is fully determined by c,,, and not by cp. Therefore, the third block in Fig. 1 is chosen as an integrator corresponding to case (iii) above. Parallel combinations of the alternatives above are possible but have so far been avoided for the sake of simplicity. Appendix C The integral (4) of the present paper must be computed repeatedly to study the time course of the oscillation. The integral is of the type Z(t) = f h(t’)g(t - t’) dr’, 0

(Cl)

where t’ < t, W)

H, h, tl and t2 in (C2) are constants and t in (Cl) denotes real time. h(t’) is illustrated in Fig. 5. g(t--t’) in (Cl) is at time t known for t’ > 0. The properties of h(t’) can be used to facilitate the computation of Z(t) for t = 0, At, 2At, etc. The following notations will be used: Z’(t) = f i’ (t’-t,)g(t-t’) t1 Q(t) = i

I

tr (t2 -tJ

dt’,

e-h(t’--l*) g(t - t’) dt’,

= 7 g(t - t’) dt’. 11

cc31 w (C5)

174

A.

JOHNSSON

AND

H.

G.

KARLSSON

Then I(r) can be written as Z(f) = P(t)+Q(t). ((35) Now P(t), Q(t) and R(t) are supposed to be known for t = t,, where a stands for actual time, and shall be computed for I = t,+At. It can be shown rather easily that P(t,+Af)

= P(t,) + $ R(t,) + ;

i f,

-; Q(t,+At) R(t,+At)

= eFhA’Q(fa) + h = R(t,) +

r f,-At

(t’+At-rr)g(t,--r’) ,,%, (t’+Ar-fi)g(t,-t’)

1 (tz-t,) f2 Ar

g(r,-r’)

dt’-

At

dt’ -

e-h(t’+At-f2) g(t.-t’) ff’ g(t,-t’) tz-At

dt’.

dt’, dr’,

(C7) ((3 w

From (C7), (C8) and (C9) it follows that P(t), Q(r) and R(t) can be computed recursively. In the present papers all computations have been performed on a digital computer and the recursive technique above has been used to calculate integral (4). It is also possible to use analogue computers to simulate the system where integral (4) is involved, provided a device capable of delaying electrical signals is available. Such a device has been built (Karlsson, Borjesson & Johnsson, 1971) and analogue simulations are planned for the future.