Balance equations can buffer noisy and sustained environmental perturbations of circadian clocks

We present a new approach to understanding how regulatory networks such as circadian clocks might evolve robustness to environmental fluctuations. The approach is in terms of new balance equations that we derive. We use it to describe how an entrained clock can buffer the effects of daily fluctuations in light and temperature levels. We also use it to study a different approach to temperature compensation where instead of considering a free-running clock, we study temperature buffering of the phases in a light-entrained clock, which we believe is a more physiological setting.


INTRODUCTION
Circadian oscillators are entrained by the daily cycles of light and temperature. It is therefore important that a clock is sensitive to their daily periodicity. On the other hand, there are very substantial stochastic dayto-day fluctuations in these environmental cycles. This can, for example, be seen in the time series for light intensity and temperature shown in figures 1 and 2. The daily fluctuations in both are substantial: the fluctuations in light have a coefficient of variation of approximately 36 per cent, while those of England's maximum and minimum temperatures in degrees centigrade over just the single month of September have a coefficient of variation 15 and 29 per cent, respectively.
The presence of such substantial noisy perturbations raises a number of questions. Our theoretical understanding suggests that in a typical model these perturbations will produce substantial fluctuations in the protein levels (e.g. figure 4). This seems to be in conflict with the need for the clock to provide robust signals to the genes that it is controlling. It therefore raises the question of whether the clock can be designed so that the daily protein variation is robustly entrained but at the same time the system manages to effectively buffer the stochastic variations.
In this paper, we explain how it is relatively easy for evolution to adjust a clock network so as to carry out such buffering of stochastic environmental fluctuations. We argue that this is potentially an important reason why an oscillator, rather than a direct measurer of the entraining signal (light or temperature) is used. The mathematical argument that we employ is very general and can be applied to a much broader class of regulatory and signalling systems and other environmental factors. It uses a combination of ideas behind balance equations (used previously to explain temperature compensation [1 -3]) and the principal component aspects of global sensitivity analysis [4 -8]. We consider fluctuations in the light level and temperature fluctuations.
To give examples of the effectiveness of the buffering mechanism for light fluctuations, we use a published model of the Arabidopsis thaliana circadian clock [7], and for temperature fluctuations, we introduce a temperature-dependent version of this model.
A related issue was addressed in a recent paper [8]. In a study involving artificial in silico evolution of clock networks, it was shown that a combination of the need to cope with multiple photoperiods and stochastic variation in the timing of dawn and dusk favoured the evolution of more loops and light inputs and greater complexity in the networks.
After considering short-term fluctuations, we turn to temperature compensation. Temperature compensation refers to the striking and defining feature of circadian clocks whereby their period only varies by a small amount over a physiological range of temperatures [9,10]. However, the exact value of the free-running period in constant conditions does not appear to have a direct selective value in the natural environment, as the clock will normally be entrained to diurnal day/ night cycles. One may therefore ask why temperature compensation has arisen during evolution. We address this question here and develop a theory to address the question of how evolution might act on forced entrained oscillators. We show that for typical systems, the phases of the proteins in the clock will vary significantly with temperature but that one can tune the clock to satisfy certain balance equations so that the changes in phase are buffered. We propose that the buffering of the free-running period characteristic of classically temperature-compensated systems is a consequence of this phase buffering. The idea that there is a connection between the free-running period and the entrained phase is not a new one. The way in which phase changes as one crosses an Arnold tongue is well understood in dynamical systems theory [5,11] and in circadian rhythms [12,13]. This relation between the period and the entrainment phase has been observed experimentally in physiologically relevant situations [14].
Throughout, we assume that our circadian clock is modelled by a set of ordinary differential equations where t is time, a vector x ¼ (x 1 , . . . , x n ) represents the state variables (namely, mRNA and protein levels) and k ¼ (k 1 , . . . , k s ) is a vector of parameter values. We assume that equation (1.1) has an attracting periodic solution x ¼ g(t,k) of period T. We study the properties of this solution. We will illustrate our results by using a model of the Arabidopsis circadian clock [7].

GLOBAL SENSITIVITY ANALYSIS AND ITS PRINCIPAL COMPONENTS
If g(t) is the solution of equation (1.1) mentioned above, then the change dg(t) in g caused by a change dk ¼ (dk 1 , . . . , dk s ) in the parameter vector k is where the linear map M is given by We can regard M as a map from the parameter space R s to a space of time series. For our purposes, the appropriate space of time series is the Hilbert space H defined in appendix A, and we also consider the subspace H 0 of H spanned by the functions @g/@k j (t) on 0 t T, j ¼ 1, . . . , s. In Rand [6], it is shown that there are a set of numbers s 1 ! s 2 ! . . . ! s s , a set of orthonormal vectors V 1 , . . . ,V s of the parameter space R s and a set of orthonormal vectors U 1 , . . . ,U s in H such that MV i ¼ s i U i , M*U i ¼ s i V i with the following optimality property: for all k ! 1, the average error given by is minimized over all orthonormal bases of H 0 . At this minimal value, e k 2 ¼ cs k 2 , where c is an absolute constant. The s i are uniquely determined and the V i and U i are, respectively, eigenvectors of MM* and M*M. Thus, the s i are the eigenvalues of M*M. If they are simple eigenvectors, then the U i and V i are uniquely determined.
M* is the adjoint to M and is given by It follows that the ijth element of M*M is given by k @g @k i ; @g @k j l  Since this is self-adjoint, it has real positive eigenvalues and these are s 1 ! s 2 ! . . . ! s s .
It follows from the above discussion that a change dk to the parameters leads to a change dg in the solution g such that and that the s i decay as rapidly as is possible for any such representation. Here, the W ij are the elements of the inverse matrix W ¼ V 21 . The speed with which the s i decay for the Arabidopsis clock model [7] is shown in figure 3. Inspection of figure 3 shows that this decay is exponential for the case for the clock model of Locke et al. [7], and this is typical [4,6,15] as also shown in the electronic supplementary material, figure S1. Note that when applying this theory, one can either apply it directly to the model parameters k j or one can apply it to the parameters h j ¼ log k j . It is often more relevant to use the second approach because, in regulatory and signalling systems, the values of two parameters may differ by an order of magnitude or more. Therefore, it is more appropriate to consider relative changes in the parameters k j than the absolute changes. The only change to the theory when considering relative changes instead of absolute ones is to replace the linear operator M above by M . D, where D ¼ diag(k) is the diagonal matrix whose diagonal is made up of the parameter values. If we denote quantities for the latter relative case using a tilde, we have Ṽ ¼ D . V, W ij ¼ k j W ij and s i ¼ s i . This follows from the fact that for small changes dk to the model's parameters, dh j ¼ dk j /k j . The scaled changes dh j also have the advantage of being non-dimensional.

BALANCE EQUATIONS FOR SUSTAINED AND DAILY FLUCTUATIONS
Now suppose that a subset fk j 1 , . . . , k j q g of the parameters depend upon a parameter p, i.e. k j m ¼ k j m ( p). Then, where the latter sum is over j in S ¼ f j 1 , . . . , j q g. Thus, if we want dg/dp ¼ 0, we require . . , s. The s equations in (3.2) are called the balance equations and the sums P j[S W ij dk j /dp are called balance sums.
Note that the ith balance equation can only be solved if, for this value of i, the W ij dk j /dp do not all have the same sign. This places a constraint on the W ij and hence on the system as a whole.
Now we come to the reason why using the above principal components U i is important. Note that in equation (3.1), each term P j[S W ij dk j /dp is multiplied by s i and therefore if the s i decrease rapidly, as usually is the case for the sort of systems that we consider, then the importance of the balance equations in ensuring dg/dp % 0 decreases rapidly as i increases, and obtaining the balance equations for just a few low i will substantially decrease jdg/dpj. Inspection of figure 3 shows that this is certainly the case for the clock model of Locke et al. [7] and this is typical [4,6,15].
For circadian oscillators, we are often particularly interested in the changes in phase of the various components of the clock. If we define the phase w m of the mth component as the time that it reaches its maximum value, then [6] dw m dp ¼ À A derivation of equation (3.3) is provided in appendix A. Thus, the balance equation to obtain dw m /dp % 0 agrees with equation (3.2). In practice, the balance sums are never zero, and the aim of the balancing is to reduce them substantially in the sense that the ratio of the sum after balancing to that before is substantially less than one.

APPLICATION TO DAILY LIGHT FLUCTUATIONS
Our aim here is to demonstrate a simple mechanism that can enable the clock to filter out the sort of substantial variation in daytime light intensity that is observed in the Harvard Forest data [17], figure 1. We denote the time-dependent light intensity by u(t). For example, u(t) might be the function that is 1 between dawn t l and dusk t d and zero elsewhere or it might be a slightly smoothed version of this discontinuous function. We model fluctuations in the light by replacing u(t) by au(t), where a fluctuates daily around its mean value of a ¼ 1.
We then consider a model where the effect of light is to modulate a subset of the terms in the differential equation so that they are functions of au(t). In general, the light input will occur in multiple terms in the form k j u(t), where k j is one of the parameters of the system. If the light level is fluctuating as above, then the term k j u(t) is replaced by k j au(t). Let S be the set of parameter indices for the parameters that occur in this way.
Suppose that this original model is not buffered against variation in daytime light intensity. To buffer it, we can assume that there is some simple regulation of the light inputs, so each of the terms k j au(t), j [ S, is replaced by k j (c j a þ d j )u(t). It is also natural to assume that c j þ d j ¼ 1 since this implies that k j (c j a þ d j ) fluctuates around k j . One can always reduce to this case by scaling the k j in advance. We call this the light-modified model.
If c j and d j are fixed (and not regarded as parameters), then the new value for W ij in this new system is c j W ij , where the latter W ij is the value in the original system. Therefore, the balance equation for absolute changes in the parameters is where the W ij are as above for the approach where one uses relative changes in parameters instead of absolute ones so that one uses h j ¼ log k j . This is a very general formulation and is, for example, directly applicable to the model of Locke et al. [7] for the Arabidopsis clock that introduces light in the way described. The light parameters k j , j [ S, are listed in table 1. The light parameters that have the greatest effect on the balance equations (by ranking s i W ij ) are LHY transcription (q 1 ), light and TOC1mediated induction (q 2 and n 4 ) of Y transcription and accumulation of protein P ( p 5 ).
We assume that the variations a are normally distributed with mean m and standard deviation s, which we denote by a N(m,s). As shown in figure 4 (and in the electronic supplementary material, figure  S2), the daily variation in the amplitude of the limit cycle solution of this model is substantial and quantitatively reflects the variation in the light amplitude. A light-modified model was constructed as described for which the left-hand sides of equation (4.1) are substantially smaller, as is given in table 1.
To implement this, the values of six of the c j were chosen fairly arbitrarily but relatively close to one and then the remaining three were calculated by solving the linear balance equations corresponding to the first three principal components. Our balanced model shows less variation in the output protein and gene levels ( figure 5 for a N(1,0.2)). For all genes and proteins, the balanced model shows less variation in the concentration levels than the unbalanced Locke model (electronic supplementary material, figure S3 Table 1. Left: parameters for the balanced model. Each light parameter k j has light intensity of the form a j L(t). For the Locke model, a j ¼ a ¼ 1, while for the balanced model, a j ¼ c j a þ d j , where c j are listed in the table and d j are d j ¼ 12c j . Right: the sum P j W ij dk j /da for i ¼ 1,. . ., 8 evaluated for the Locke model (dk j /da ¼ 1) and the balanced model (dk j /da ¼ c j ). The corresponding singular values s i are plotted in figure 3.  and table S1). The phases of the balanced model show less variation for all components, except the LHY gene and proteins, which are of same order as those of the Locke model (electronic supplementary material, table S2).

APPLICATION TO DAILY TEMPERATURE FLUCTUATIONS
We consider a model with fixed but differing day and night temperatures, T D and T N . Below, we derive a set of balance equations that eliminate the effects of temperature fluctuations. We assume that the temperature enters the model through a set of parameters k j , j [ S, which are temperature sensitive. A standard assumption for the temperature dependence of each of model parameters k j is that it is similar to that for rate constants of chemical reactions and is described by the Arrhenius equation. This expresses the dependence of the rate constant k j on the temperature T and activation energy E j as k j ¼ A j exp(2E j /RT), where A j is a constant specific to the individual parameter and R is the gas constant (8.314472 Â 10 23 kJ mol 21 K 21 ). In our case, we need to deal with the fact that we have different night and day temperatures and therefore we assume where t is in the daytime if t l t mod t , t d , where [t l , t d ] denotes the range of day hours, and otherwise t is in the night. Moreover, e and h denote the fluctuations in day and night temperatures, T D and T N , respectively. We consider the following first-order Taylor series expansion: where k j,0 D ¼ A j exp(2E j /RT D ) and k j,1 D ¼ 2 k j,0 D E j /RT D and similarly for the night parameters. Note that this is a very good approximation and higher order terms can be neglected, since for a parameter of order O(0.1) (i.e. the order of a large number of Locke parameters), with sensible activation energy E % 50 kJ mol 21 , the coefficient of the second-order term is of order O(0.001).
From the observation in equation (3.1), day temperature variations e will have the following effects on changes to the solution g: and the following changes to the phases: Similar expressions can be derived for night temperature fluctuations. Together, these give rise to two sets of balance equations, Or alternatively, since, by the above, k j,1 In order to make a model temperature dependent, we have to define the dependence of its parameters on temperature. It is likely that all parameters in a regulatory system are actually temperature dependent, but this temperature dependence will have little effect for a parameter k j if all the sensitivities S ij ¼ s i W ij are small for that value of j. Therefore, in order to keep the model reasonably tractable, we will only introduce temperature into those variables with a significant sensitivity.
Moreover, to determine the relative importance of parameters, since values of some parameters differ by an order of magnitude or more, it is more appropriate to compare relative changes of parameters and to use the log parameters and the W ij as described in the last paragraph of §2. We selected parameters that ranked highest when ordered by max i¼124 js i W ij j. The plot of top 20 is shown in figure 6. In this selection, we include only parameters that come linearly in the model and exclude Hill coefficients as these are not rate constants. Parameter n 6 has the highest max i¼124 log 10 js i W ij j. Only seven other parameters have These eight parameters will be temperature sensitive.
Note that the key parameters to control the buffering were parameters linked to PRR7/9 components, namely LHY-dependent transcription (n 6 ), mRNA degradation (m 16 ), translation ( p 6 ) and cytoplasmic protein degradation (m 17 ) as well as the nuclearcytoplasmic transport (r 12 ). Aside from these, the list also includes TOC1 light-independent transcription (n 2 ), mRNA degradation (m 4 ) and LHY mRNA degradation (m 1 ). Each of the parameters k j was separated into a morning and an evening term, k j ¼ k j D u þ k j N (1 2 u), where morning and evening terms were adjusted to be a fixed percentage higher or lower than the Locke model values. Note that the term u is the light function. For most parameters, we chose this to be 10 per cent, but for m 17 , we had to choose a significantly smaller value of 1 per cent, so that gene and protein waveforms of the temperature-sensitive model would match those of the Locke model as closely as possible. We confirmed that the parameters had a strong effect on the model by verifying that several of them appeared in the top 10 per cent of parameters ordered by their magnitude of max i¼124 log 10 js i W ij j. We chose night and day temperatures to be T ¼ 285.15 K and T ¼ 290.15 K, to be close to the temperature data from UK Met Office (figure 2), with mean maximum and minimum temperatures T ¼ 292.79 K and T ¼ 283.92 K.
We calculated the energies of the new model and checked the balance equations. We tried several different combinations of energies, then recalculated the morning and evening parameters and then computed the balance equations for the model. From the starting Locke model, we evolved two models with differing values of energies (table 2). We labelled them a . Power sensitivity spectrum of the Locke model. Each group of bars corresponds to the values of log 10 js i W ij j for a parameter k j . These are only plotted for those i for which log 10 js i W ij j is significant (here that is i ¼ 1 -4). The parameters k j are ordered by max i¼124 log 10 js i W ij j and only 20 most sensitive ones are plotted. Not taking into account Hill terms g 2 and g 7 (for reasons outlined in the main text), parameter p 6 has the highest max i¼124 log 10 js i W ij j. Only seven other parameters ( parameters up to and including m 1 ) have sensitivity higher than 30 per cent of the maximum, and they will be candidates for temperature-sensitive parameters. balanced and an unbalanced model according to their fit to the balance equations (table 3). The balanced model was chosen so that originally it had activation energies of about 30 kJ mol 21 for almost every component except for that of the parameter m 17 whose activation energy was chosen to be significantly lower to make balancing easier. Since later we adjusted the temperature range to fit with the data seen in figure 2, these values appeared slightly smaller. The unbalanced model was made by changing four energies from the list of the balanced model, in order to make a worse fit to the balance equations.
In fact, our initial temperature-sensitive model gave the best balance equations, so it was chosen as the balanced model. The balanced model shows less variation in peak concentrations and phase variations than the unbalanced model (figures 7 and 8 and electronic supplementary material, tables S3 and S4).

APPLICATION TO TEMPERATURE COMPENSATION
Instead of temperature compensation on a free-running clock (i.e. dp/dT % 0), we are interested in temperature compensation of the light-entrained clock, which we believe is a more physiological setting. The aim is to minimize protein and phase changes in the context of changing temperature, dw/dT % 0. We therefore consider models that are entrained to the day-night cycle rather than the free-running clock. This hypothesis [16] that a balance of opposing reactions could allow temperature compensation was first tested by Ruoff using a simple model for an oscillatory feedback loop [1 -3]. He used an Arrhenius representation for temperature dependence and deduced a balance equation for the local period slope dp/dT in terms of the activation energies E j and control coefficients for each of the parameters.
We also assume that parameters k j , j [ S, are temperature dependent and describe them by Arrhenius equations, k j ¼ A j exp(2E j /RT), as described above. The temperature t is in kelvins in the range 285.15 K T 300.15 K. Activation energies, E j , must be in the range 1 kJ mol 21 E j 150 kJ mol 21 . We insist that some of the activation energies are substantial because otherwise the parameters only have weak dependence upon temperature. In fact, for a given balanced system, since these energies enter linearly into the balance equation, scaling them by a factor just scales the divergence from perfect balance by that factor. Table 3. The balance sums P j W ij k 0,j D E j and P j W ij k 0,j N E j and the corresponding singular values s i for the unbalanced model (UB) and the balanced model (B). To get the true sums, divide each column by 1/RT 2 , where T D ¼290.15 K and T N ¼285.15 K.  We now apply equation (3.2) where p is replaced by temperature T. From the relation dk j /dT ¼ k j E j /RT 2 , we deduce that the balance equations are . ., s. When the s i decrease rapidly, we need only consider these balance equations for the first few i in order to get dg/dT or dw m /dT small. Note that the ith equation can only be solved if, for this value of i, the W ij do not all have the same sign. This is because k j E j is always positive. This places a constraint on the W ij and hence on the system as a whole. Only certain networks can be balanced.
We assume that the model parameters from Locke et al. [7] correspond to a model at T ¼ 288.15 K. Since this model does not include temperature, we selected the temperature-sensitive parameters at each temperature end by checking the parameter sensitivity spectrum, as described in §5.
We consider three models in which temperature dependence is inserted into the Locke model. In two of them (models 1 and 2), the parameter values at T ¼ 288.15 K are as in the original model but although they have the same activation coefficients A j , they have different activation energies E j (table 4).
The energies for model 1 were selected so that they were substantial and roughly satisfied the balance equation (table 6). We then iteratively calculated the W ij 's for the new model and more exactly rebalanced the equation. The iterative procedure converged reasonably well. For model 2, the energies were modified so that the balance equations were not well satisfied (table 6). Some energies were decreased as well as increased (table 4). We could easily do this as the two models share W ij 's at T ¼ 288.15 K.
The temperature compensation in model 2 is substantially worse than that in model 1, demonstrating the importance of balancing. We find that model 1 is better at local temperature compensation than model 2, with less variation in its protein concentrations and phases from the original published model [7], cf. TOC1 protein times series in figure 9, and also it fares better at global temperature compensation (electronic supplementary material, tables S8 and S9). Moreover, model 1 can temperature compensate in constant light conditions (figure 10), with model 2 becoming arrhythmic outside a narrow temperature range. Model 3 shares activation energies of model 1 but has different activation coefficients A j . Thus, its parameter values do not agree exactly with the Locke model at any temperature (table 5). This model is better balanced than either of the other two (table 6) and is significantly better compensated for phase than the other models (tables 7 and 8 and electronic supplementary material, table S10). Moreover, in continuous light conditions, its period is much better compensated than model 2 and slightly better than model 1 (figure 10).

DISCUSSION
We have attempted to show how a combination of ideas behind balance equations and the principal component aspects of global sensitivity analysis gives a new  approach to understanding how to design regulatory networks that are buffered against certain fluctuating environmental perturbations. As well as presenting some concrete examples of applications to circadian clocks, we have described the general theory behind this. From this, it is clear that it can be applied to a much broader class of regulatory and signalling systems and environmental factors other than light and temperature.
If buffering of a particular environmental fluctuation has significant selective pressure for the organism in question, then we can interpret this selective pressure as acting on the balance equation. There will then be selective pressure on the quantities that make up the left-hand side of the equation (e.g. system parameters, activation energies) to make them balance to zero. Understanding this makes it much clearer how evolution can act to achieve what appear to be quite complex tasks.
Temperature compensation has been one of the driving dogmas of circadian biology and has been interpreted in terms of the constancy or near-constancy of the freerunning period of the circadian clock under changing temperature. However, it is not clear how evolution acts on the free-running period since in physiological conditions the clock is entrained to the day-night cycle and has a period of 24 h. We show that through certain balance equations it is possible to buffer the changes in phase over relatively large ranges of temperature. Although we do this for a specific example, the mathematical approach suggests that this buffering should be possible for an extremely broad range of clock models. We suggest that the observed near-constancy of the free-running period is a consequence of the nearconstancy of the phases of the entrained clock.
To balance a balance equation, it is necessary that the terms making up the equation do not all have the   peak times trough times

METHODS
All the computations were carried out using Matlab and XPPAUT. In particular, the global sensitivity calculations were done using the Matlab-based Time Series Sensitivity Analysis Package available from http:// www2.warwick.ac.uk/fac/sci/systemsbiology/software/, and the period calculations were performed using XPPAUT, available from http://www.math.pitt.edu/ bard/xpp/xpp.html.
We are particularly grateful to Paul Brown who contributed extensively to the software used by us in preparing this paper. We are also grateful to the ROBuST team for many useful discussions on this topic and particularly to Andrew Millar. The comments of two anonymous referees were very helpful. This research was funded by BBSRC SABR grant BB/F005261/1 (ROBuST Project) and EU BIOSIM Network Contract 005137. D.A.R. is also funded by EPSRC Senior Research Fellowship EP/C544587/1.

A.2. Derivation of equation (3.3)
This was presented in Rand [6]. Let t ¼ f m (k) be the time when the concentration of the mth component, g m , is at a maximum or a minimum value. Hence, f m (k) satisfies g˙m(f m (k),k) ¼ 0. Differentiating both sides with respect to k j and rewriting, @f m @k j ¼ À ð@ _ g m =@k j Þðf m ðkÞ; kÞ € g m ðf m Þ ¼ À ð@=@tÞj t¼f m ð@g m =@k j Þðf m ðkÞ; kÞ € g m ðf m Þ ¼ À P i s i W ij _ U i;m ðf m Þ € g m ðf m Þ :