MODELING SOME PROPERTIES OF CIRCADIAN RHYTHMS

Mathematical models have been very useful in biological research. From the interaction of biology and mathematics, new problems have emerged that have generated advances in the theory, suggested further experimental work and motivated plausible conjectures. From our perspective, it is absolutely necessary to incorporate modeling tools in the study of circadian rhythms and that without a solid mathematical framework a real understanding of them will not be possible. Our interest is to study the main process underlying the synchronization in the pacemaker of a circadian system: these mechanisms should be conserved in all living beings. Indeed, from an evolutionary perspective, it seems reasonable to assume that either they have a common origin or that they emerge from similar selection circumstances. We propose a general framework to understand the emergence of synchronization as a robust characteristic of some cooperative systems of non-linear coupled oscillators. In a first approximation to the problem we vary the topology of the network and the strength of the interactions among oscillators. In order to study the emergent dynamics, we carried out some numerical computations. The results are consistent with experiments reported in the literature. Finally, we proposed a theoretical framework to study the phenomenon of synchronization in the context of circadian rhythms: the dissipative synchronization of nonautonomous dynamical systems. 2010 Mathematics Subject Classification. Primary: 92B25, 34C15; Secondary: 68U20.

1. Introduction.It is a well experimental fact that under constant environmental conditions, the majority of the physiological variables of living organisms display endogenous, sustained oscillations with a period close to, but different from, 24 h, known as circadian rhythms.These rhythms allow organisms to predict, instead of merely react to cyclic changes in the environment.They also provide different species the possibility of a temporal organization.Because these circadian cycles can be synchronized with external signals, organisms are also able to adjust themselves to exogenous rhythms such as the day/night cycles.See, for example [24], [21].Daily rhythms are controlled by a network of circadian clocks that are able to reset themselves into a new phase when appropriate input signals are present.This occurs in order to impose the corresponding circadian signal on effector organs, mainly peripheral tissues.In mammals, the organization of the circadian network appeared to be hierarchical, with the master clock localized in the suprachiasmatic nuclei at the top of the hierarchy [6].Recent findings, however, have challenged the dominant role of these structures in physiological regulation, and there is increasing evidence that a close interaction between central and peripheral clocks is necessary to maintain robust circadian rhythms in functions and metabolism [3].In the circadian pacemaker system of birds, a hypothalamic region has also been detected as a possible equivalent to the mammalian suprachiasmatic nuclei [13].On the other hand, circadian biology of invertebrates has led to a different hypothesis: a circadian system of pacemakers is distributed in which no master clock has been identified.The model of distributed clockwork involves oscillators such as retinular cells, neurosecretory systems in the optic lobes, and putative brain pacemakers [33], [28], [1].Independently of the detailed organization of the circadian system of different phyla and even for the differents species, it is well accepted that pacemaker´s cells must oscillate in synchrony to produce a reliable global rhythm (See, for example, [14], [37], [31], [36], [2], and [11]).In these works a pacemaker is regarded as a network of coupled oscillators.It seems that the proposed models differ from each other in several respects: some play special attention in the geometric organization of the network, others in the molecular data or in the neurotransmitter involved, and yet others in the hierarchical organization of the oscillators.Because synchronization is a universal phenomenon [25], and thus independent of the nature of each implicated oscillator, we argue that besides modeling the details of the constituents of the network, it is also important to understand the essential and common features, universal mechanisms, of resynchronization in order to put them in the context of circadian rhythms.If we are able to achieve this goal, we could also deduce some well-known phenomena of the circadian rhythms such as the entrainment or the anatomical (spacial) organization of the network.In this sense it is our interest to study the main process underlying synchronization in the pacemaker of circadian systems.We propose a general framework to understand the emergence of synchronization as a robust characteristic of some cooperative systems.The coupling could be via intrinsic factors, for instance cell-to-cell communication.In this case all components are similar, there is no preferred element.Although we give a more general setting and refer to rigorous results after providing an abstract framework, we only consider simple examples in which explicit computations can be made and where the cooperative structure of the system gives rise to synchronization.
2. Mathematical model.In general, rhythms and particularly biological rhythms, arise as a temporal organization that emerges beyond a critical point of instability of a nonequilibrium steady state.This type of nonequilibrium self-organization represents dissipative structures that can be mantained only by the energy dissipation associated with the exchange of matter between the biological system and its environment.When a periodic phenomenon occurs in constant environmental conditions, as in the case of circadian systems, we have the clearest sign that the biological system operates beyond a point of instability.Endogenous rhythms, produced by a system and not by its environment, are indeed a sign of instability [12].Sustained oscillations due to a stable limit cycle can be viewed as temporal dissipative structures.That is why throughout this work we are going to use models exhibiting this behavior.
In order to fix ideas we will consider a system of n nonlinear coupled oscillators.We will assume that the oscillators are identical in the sense that they follow the same dynamics and that the corresponding parameters are taken to be equal.This is of course a reasonable assumption, but which should be modified later.For instance, in order to investigate what happens when the parameters have slight variations around an average value, which would be more realistic.In mathematical terms we have a system of n second order nonlinear differential equations: where x i (i = 1, • • • , n) represents the state of the i−th oscillator.In absence of any coupling, the dynamics of each oscillator is determined by the same function f .The dynamics of f must exhibit a stable limit cycle behavior.In principle it is not necessary to impose any specific functional form for f but, in order to show a specific example, we consider it to be a van der Pol oscillator.Thus: In the particular case of circadian systems the coupling mechanisms of the cells are complex and subtle.In fact, at present, this subject constitutes a very important line of investigation [37], [36], [2].In our model the coupling constants a ij represent the influence of the j−th oscillator on the i−th one.Since we are assuming cooperative behaviour, it is natural to take them all positive.The relative values of the a ij , as we can see in 2.1, can reflect different geometric configurations and different types of interactions.In fact, not only planar configurations, but also three dimensional arrays could be considered, e.g. a cube.However, for the sake of simplicity we consider in this work only two dimensional systems and leave more complex cases for future work.Regarding the relative magnitude of the matrix coefficients the basic assumption is made that they are inversely proportional to the square of distance between the oscillators.This is natural in some contexts when the oscillators interactions are mediated by physical fields, for instance the luminous pulse of fireflies or in general interactions in which dissipation can be neglected.In this case, standard conservation of energy considerations lead to this power decay of the intensity of the signal.However, in other biological contexts in which other signaling mechanisms are responsible for the interaction, e.g.cell to cell communication in tissues and organisms.For such systems, a different dependence of the coefficients would be appropriate.Studying different geometric configurations such as a linear or a square array might be carried out by setting specific values.This point, as we will see in section 2.2, is important for undersanding the emergence of the circadian temporal organization during the ontogeny.All this allows us to study the emergence of different spatio-temporal patterns of synchronization and to interpret them in the context of circadian rhythms.The µ is a positive constant that provides a measure of the overall strength of the coupling.
Rewriting the n coupled van der Pol oscillators (considering 1 and 2) in a matrix form we obtain: ..,n is the non negative R n×n matrix representing the oscillators networks structure, and finally I is the identity matrix in R n .
2.1.Effect of network architecture on synchronization.We only consider two very simple examples in which explicit computations can be made and the cooperative structure of the system gives rise to synchronization.As it was pointed out before, other more complex situations both regarding the geometry and the dependence of the coefficientes reflecting the nature of the interactions among oscillators can be considered.In many cases it is also natural to assume that there is no influence of the oscillator on itself, i.e a ii = 0 and that the influence of the j−th oscillator on the i−th is equal to the influence of the i−th oscillator on the j−th one.That is because matrix A of equation 3 is symmetric with vanishing entries on the principal diagonal.We consider four oscillators organized in two different configurations: a single line (all equally spaced) and in a square (Figure 1).2.1.1.Single file.Supposing that the distance between adjacent oscillators is equal to 1 and considering 3 we obtain: with n = 4 The parameters γ and ω affect the characteristics of each van der Pol oscillator; if γ is small (γ = 0.2) we have a quasilinear oscillators whereas if it is large then we have a relaxation oscillator (γ = 5); ω is related with the frequency of each isolated oscillator.In Figure 2 we show the temporal evolution of x 1 from van der Pol equation (see 2) when γ = 2, ω = 3, and using (0, 0) as initial condition.Near limit cycle frequency of oscillation is 0.25.These parameters, as all others in the manuscript, were chosen for simplicity in the numerical calculations and are only qualitative.In fact, the next step regarding the connection with real biological systems is the quantitative tuning of realistic parameters in a specific experimental situation.  .We used γ = 2, ω = 3, and (0, 0) as initial conditions.The frequency of oscillation near the limit cycle is 0.25.
In Figure 3 we show the simulation for equation 4 for two different initial conditions.We chose µ = 1 (in section 2.2 we are going to vary this parameter).
Figure 3. Temporal evolution of four coupled oscillators when they are arranged on a row.In the top part we used the initial condition (0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1) and in the bottom part (0, 0, 0, 0.5, 0, 0, 0, 0).The emergent dynamics depends on the initial conditions.In the top part of the figure it can be observed a synchronized behavior between oscillators 1 and 2, and 3 and 4 but not among all.In the bottom part of the figure it can be observed that the four oscillators are synchronized (although not with the same phase).Note the change in frequency of the emergent dynamics with respect both cases presented here, but also with respect the isolated oscillator (Figure 2).2.1.2.Square.We consider a square with side of length equal to 1.We obtain the following equation: In Figure 4 we show numerical solutions of equation 5 for two different initial conditions.We used the same parameters as for the linear arrangement.In all cases simulations were performed with Mathematica 7. In the top part we used the initial condition (0, 0, 0, 1, 0, 0, 0, 0) and in the bottom part (0, 0, 0, 0.5, 0, 0, 0, 0).As in the single line case.the emergent dynamics depends on the initial conditions.In the top part of the figure it can be observed a synchronization among all oscillators.In the bottom part of the figure can be observed that the four oscillators get in synchrony (although not in the same phase).Note the change in frequency of the emergent dynamic with respect both cases presented here, but also with respect the isolated oscillator (Figure 2).

Discussion and further remarks.
Little is known about the connectivity and topological properties of the pacemaker cellular network, but, undoubtedly, it is a central problem [14].Computational biology starts with the complex interactions we already know about, and uses computer-aided mathematics to explore the consequences of those interactions.We intend to create models in which networks are specified in terms of elements and interactions, the network "topology", but the numerical values that quantify those interactions (the parameters) are deliberately varied over wide ranges.As a result, the study of such networks must focus not on the exact values of outputs, but rather on the qualitative behavior.By investigating how such behavior changes for different parameter sets, when "exploring the parameter space", one starts to assemble a comprehensive picture of all the possible qualitative dynamics that a network can produce.If one of such behavior seems useful (to the organism), it becomes a candidate for explaining why the network itself was selected, i.e., it is seen as a potential purpose for the network.If experiments subsequently support assignments of actual parameter values to the range of parameter space that produces such behavior, then the potential functionality becomes viable [17].Characterization of anatomical and functional connectivity in other regions of the brain of a mammal, different of the pacemaker, e.g.cortex, revealed small-world network properties.For example see [20], [38], [4].Small-world network is a type of mathematical graph in which most nodes are not neighbors of one another, but most nodes can be reached from every other by a small number of steps.Specifically, a small-world network is definided to be a network where typical distance L between two randomly chosen nodes (the number of steps required) grows proportionally to the logarithm of the number of nodes N in the network, that is L ∝ log N [8].The speculation that this architecture is dynamically advantageous, for example more synchronizable or error-tolerant, need to be sharpened and then confirmed or refused mathematically [34].
In this section, we only vary the topology of the network but not the connectivity.To make a study including a more detailed mathematical presentation of the different qualitative behavior of the network or the ability to be synchronized constitutes part of our work in progress.Also we have to consider the connectivity of the network as another important element.

2.2.
Effect of the strenghth of interaction on the network.As we anticipate in a previous section, we are going to check the effect of modifying the value of µ in the dynamics of equation 1 when we use equation 2. We recall that, µ is a positive constant providing a measure of the overall strength of the coupling; when µ = 0 we have n uncoupled oscillators, and when µ → ∞ we have the maximum coupling among the oscillators of the network.
To show the qualitative behavior of the network for differents values of µ we use a particular arrangement: the square; parameter values are the same of those of the section 2.1.2.In Figure 5 we show the numerical simulation for µ = 0.2 and µ = 1.1 (as in the past cases γ = 2, and ω = 3).We can observe several things: 1.After a transient time, oscillators display a synchronized behavior.
2. The frequency of oscillation of the emergent dynamic is different from the natural frequency of each oscillator.3. When µ rises the frequency of oscillation of the emergent dynamic diminishes.2.2.1.Other considerations.One approach to the problem of the origin of circadian oscillations involves the study of the ontogeny of the corresponding rhythm.For the expression of a circadian pattern it is necessary that the anatomical substrate reaches maturity and establishes the necessary structural and physiological relationships between its parts.It is expected that during development, different structures and functions begin to show some temporal organization that eventually will acquire circadian characteristics.This implies that underlying the sense of time of the organism, there are present a number of changes in its anatomical and functional organization.Indeed, during all the ontogeny process, the organism must exhibit successive changes in the period, relative amplitude and activity level of their rhythms.It can be assumed that variation in these parameters involves changes in the main elements of a circadian system, namely the structures implicated in the generation, expression and synchronization of the rhythm.In both cases we used the same initial condition: (0, 0, 0, 1, 0, 0, 0, 0).Note the change in frequency of the emergent dynamic.
In order to understand the mechanisms underlying circadian systems, in our team work the crayfish has been used as the biological model for many years.
When an adult crayfish with no other experimental modification is kept under complete darkness (except for test light pulses) and the ambiental temperature is controlled and fixed to 20 • C, a clear circadian rhythm is obtained by measuring the amplitude of the visual photoreceptors' electrical response to light (electroretinogram, ERG).The ERG was recorded by means of a probing electrode implanted into the cornea to lead the photoreceptor's electrical response.This was elicited by a 200 lux flash of 1 ms duration, produced at 20 min intervals.The period of oscillation of amplitude of the visual photorecepto's electrical response is about 23.4 h and shows an amplitude near to 100 µV .The ERG is just the sum of the activity of visual photoreceptors and its amplitude depends on both the probability that a photon reaches the retina and the photoreceptor cells excitability.
If we follow the same experimental protocol to measure the ERG rhythm but in very young crayfish (immediately after hatching), we found that the ERG amplitude is very low (∼ 4 µV ) and shows clear high frequency (ultradian) periodical fluctuations with a period ranging from 15 min to 4 h.A four week old crayfish expreses a higher ERG amplitude (∼ 50 µV ) and, for the first time, the presence of a circadian oscillation.The analysis of these recordings also reveals an ultradian rhythm similar to those recorded from the youngest crayfish, but now superimposed on an oscillation with circadian (although somewhat irregular) periods.Older animals (around 5 months after hatching) show a progressive increment in ERG amplitude and period length, as well as a progressive disappearance of the high frequency cycles.From this moment on, amplitude and oscillation period show the characteristic values for this species (Figure 6).One is tempted to propose that during ontogeny, before the complete maturation of the neuroendocrine system, there are only ultradian frequency oscillators, which progressively vanish until they disappear under the influence of some circadian neurosecretions released principally from the sinus gland (in the eyestalk).The appearance of a circadian rhythm coincides with the maturation of the neuroendocrine system.[9].Moreover, the 4 h pattern found in the youngest crayfish is similar to the pattern found in the ERG amplitude recorded from the isolated eyestalk of an adult crayfish [29], and similar to the pattern obtained from long term recordings of the ERG from crayfish deprived of one or both sinus glands show a very irregular circadian rhythm with superimposed ultradian cycles [21], [22].Another series of experiments reinforces the idea of ultradian oscillations underlying circadian oscillation are reported in [10]: long term recordings of the ERG were made in the presence of different doses of deuterium oxide (D 2 O); the results showed a lengthening of both circadian and ultradian periods.Moreover, there was a direct relationship between the D 2 O dose and the length of period in both circadian and ultradian cycles.
Unfortunately, we don't have a conclusive explanation of the physiological role of the ultradian cycles in the ERG circadian rhythm.One is tempted to propose that in ontogeny, before the complete maturation of the neuroendocrine system, there are only ultradian oscillators.In our model this corresponds to the case when µ has a small value, i.e. when the coupling among oscillators is weak; in adult animals resurgence of ultradian oscillation coincides with exogenous perturbations (in our model when we modify the state vector; data not shown).We propose that circadian oscillation is an emergent behavior from the interaction of several ultradian oscillators, and that the interaction among these oscillators is achieved by means of a neurohormone, specifically the pigment dispersing hormone (PDH): a minimum concentration of this hormone is required for the crayfish to show circadian oscillations.In our model this corresponds with a enough large value of µ.
3. Mathematical analysis.In this section we present for the sake of completness some remarks on rigorous results regarding synchronization and quote one of them in order to exemplify how precise statements can be made.The result we quote is not directly applicable to our case, but it is still simple to formulate.There are however, more general results and we comment on them later.
Suppose that we have two autonomous ordinary differential equations in R d with parameter λ > 0 also satisfies a one-sided dissipative Lipschitz conditions and, hence, it also has a unique equilibrium (x λ , ȳλ ), which is globally asymptotically stable.This phenomenon is known as synchronization for the coupled dissipative deterministic system [19].The problem of synchronized behaviour with strong dissipation was studied in [16], [5] and references therein.The theorem they establish is: Theorem 1.Let x λ (t), y λ (t) be the solution of the coupled system

Figure 1 .
Figure 1.Two different geometric configurations of four nonlinear oscillators.In the top part it is represented the straight line configuration and in the bottom part the square configuration.In all cases arrows indicate interaction between oscillators (each one numbered).The strength of the interaction is the inverse of the square of the distance.

Figure 4 .
Figure 4. Temporal evolution of four coupled oscillators when they are arranged in square.In the top part we used the initial condition (0, 0, 0, 1, 0, 0, 0, 0) and in the bottom part (0, 0, 0, 0.5, 0, 0, 0, 0).As in the single line case.the emergent dynamics depends on the initial conditions.In the top part of the figure it can be observed a synchronization among all oscillators.In the bottom part of the figure can be observed that the four oscillators get in synchrony (although not in the same phase).Note the change in frequency of the emergent dynamic with respect both cases presented here, but also with respect the isolated oscillator (Figure2).

Figure 5 .
Figure 5. Temporal evolution of four coupled oscillators when they are arranged in square for two values of the parameter µ.In the top of the figure we use µ = 0.2 and in the bottom part µ = 1.1.In both cases we used the same initial condition: (0, 0, 0, 1, 0, 0, 0, 0).Note the change in frequency of the emergent dynamic.

Figure 6 .
Figure 6.Ontogeny of the ERG circadian rhythm.See text for explanation.