Abstract
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 freerunning clock, we study temperature buffering of the phases in a lightentrained clock, which we believe is a more physiological setting.
1. 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 daytoday 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 temperaturedependent 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 shortterm 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 freerunning 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 freerunning period characteristic of classically temperaturecompensated systems is a consequence of this phase buffering. The idea that there is a connection between the freerunning 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 1.1where 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].
2. Global sensitivity analysis and its principal components
If g(t) is the solution of equation (1.1) mentioned above, then the change δg(t) in g caused by a change δk = (δk_{1}, … , δk_{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 ℝ^{s} to a space of time series. For our purposes, the appropriate space of time series is the Hilbert space ℋ defined in appendix A, and we also consider the subspace ℋ_{0} of ℋ 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 σ_{1} ≥ σ_{2} ≥ … ≥ σ_{s}, a set of orthonormal vectors V_{1}, … ,V_{s} of the parameter space ℝ^{s} and a set of orthonormal vectors U_{1}, … ,U_{s} in ℋ such that MV_{i} = σ_{i}U_{i}, M*U_{i} = σ_{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} = cσ_{k}^{2}, where c is an absolute constant. The σ_{i} are uniquely determined and the V_{i} and U_{i} are, respectively, eigenvectors of MM* and M*M. Thus, the σ_{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 where It follows that the ijth element of M*M is given by Since this is selfadjoint, it has real positive eigenvalues and these are σ_{1} ≥ σ_{2} ≥ … ≥ σ_{s}.
It follows from the above discussion that a change δk to the parameters leads to a change δg in the solution g such that and that the σ_{i} decay as rapidly as is possible for any such representation. Here, the W_{ij} are the elements of the inverse matrix W = V^{−1}. The speed with which the σ_{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 η_{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 · Δ, where Δ = 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 Ṽ = Δ · V, W̃_{ij} = k_{j}W_{ij} and σ̃_{i} = σ_{i}. This follows from the fact that for small changes δk to the model's parameters, δη_{j} = δk_{j}/k_{j}. The scaled changes δη_{j} also have the advantage of being nondimensional.
3. Balance equations for sustained and daily fluctuations
Now suppose that a subset {k_{j1}, … , k_{jq}} of the parameters depend upon a parameter p, i.e. k_{jm} = k_{jm}(p). Then, 3.1where the latter sum is over j in 𝒮 = {j_{1}, … , j_{q}}. Thus, if we want dg/dp = 0, we require 3.2for i = 1, … , s. The s equations in (3.2) are called the balance equations and the sums ∑_{j∈𝒮} 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 ∑_{j∈S} W_{ij}dk_{j}/dp is multiplied by σ_{i} and therefore if the σ_{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 dg/dp. 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 φ_{m} of the mth component as the time that it reaches its maximum value, then [6] 3.3A derivation of equation (3.3) is provided in appendix A. Thus, the balance equation to obtain dφ_{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.
4. 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 timedependent light intensity by θ(t). For example, θ(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 θ(t) by αθ(t), where α fluctuates daily around its mean value of α = 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 αθ(t). In general, the light input will occur in multiple terms in the form k_{j}θ(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}θ(t) is replaced by k_{j}αθ(t). Let 𝒮 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}αθ(t), j ∈ 𝒮, is replaced by k_{j}(c_{j}α + d_{j})θ(t). It is also natural to assume that c_{j} + d_{j} = 1 since this implies that k_{j}(c_{j}α + d_{j}) fluctuates around k_{j}. One can always reduce to this case by scaling the k_{j} in advance. We call this the lightmodified 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 4.1where the W̃_{ij} are as above for the approach where one uses relative changes in parameters instead of absolute ones so that one uses η_{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 ∈ 𝒮, are listed in table 1. The light parameters that have the greatest effect on the balance equations (by ranking σ_{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 α are normally distributed with mean μ and standard deviation σ, which we denote by α ∼𝒩(μ,σ). 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 lightmodified model was constructed as described for which the lefthand 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 α ∼ 𝒩(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 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).
5. 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 ∈ 𝒮, 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(−E_{j}/RT), where A_{j} is a constant specific to the individual parameter and R is the gas constant (8.314472 × 10^{−3} kJ mol^{−1}K^{−1}). In our case, we need to deal with the fact that we have different night and day temperatures and therefore we assume 5.1where t is in the daytime if t_{l} ≤ t mod τ < t_{d}, where [t_{l}, t_{d}] denotes the range of day hours, and otherwise t is in the night. Moreover, ε and η denote the fluctuations in day and night temperatures, T_{D} and T_{N}, respectively. We consider the following firstorder Taylor series expansion: 5.2where k_{j,0}^{D} = A_{j} exp(−E_{j}/RT_{D}) and k_{j,1}^{D} = − 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 𝒪(0.1) (i.e. the order of a large number of Locke parameters), with sensible activation energy E ≈ 50 kJ mol^{−1}, the coefficient of the secondorder term is of order 𝒪(0.001).
From the observation in equation (3.1), day temperature variations ε will have the following effects on changes to the solution g: 5.3and the following changes to the phases: 5.4
Similar expressions can be derived for night temperature fluctuations. Together, these give rise to two sets of balance equations, 5.5Or alternatively, 5.6since, by the above, k_{j,1}^{D} = −k_{j,0}^{D}E_{j}/RT_{D}^{2}. 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} = σ_{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=1−4} σ_{i}W̃_{ij}. 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=1−4}log_{10}σ_{i}W̃_{ij}. Only seven other parameters have sensitivity higher than 30 per cent of the maximum. These eight parameters will be temperature sensitive.
Note that the key parameters to control the buffering were parameters linked to PRR7/9 components, namely LHYdependent 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 lightindependent 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}θ + k_{j}^{N}(1 − θ), where morning and evening terms were adjusted to be a fixed percentage higher or lower than the Locke model values. Note that the term θ 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 temperaturesensitive 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=1−4}log_{10}σ_{i} W̃_{ij}. 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 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^{−1} 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 temperaturesensitive 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).
6. Application to temperature compensation
Instead of temperature compensation on a freerunning clock (i.e. dp/dT ≈ 0), we are interested in temperature compensation of the lightentrained clock, which we believe is a more physiological setting. The aim is to minimize protein and phase changes in the context of changing temperature, dφ/dT ≈ 0. We therefore consider models that are entrained to the day–night cycle rather than the freerunning 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 ∈ 𝒮, are temperature dependent and describe them by Arrhenius equations, k_{j} = A_{j} exp(−E_{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^{−1} ≤ E_{j} ≤ 150 kJ mol^{−1}. 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.
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 6.1for i = 1, …, s. When the σ_{i} decrease rapidly, we need only consider these balance equations for the first few i in order to get dg/dT or dφ_{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 temperaturesensitive 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).
7. 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 lefthand 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 nearconstancy of the freerunning period of the circadian clock under changing temperature. However, it is not clear how evolution acts on the freerunning 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 nearconstancy of the freerunning 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 same sign. For example, for temperature compensation, we require that the relevant quantities W_{ij} do not all have the same sign. This puts constraints on the network structure, and this leads to a prediction about what network structures can be expected.
8. Methods
All the computations were carried out using Matlab and XPPAUT. In particular, the global sensitivity calculations were done using the Matlabbased 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.
Acknowledgements
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.
Footnotes

One contribution of 16 to a Theme Issue ‘Advancing systems medicine and therapeutics through biosimulation’.
 Received October 7, 2010.
 Accepted November 4, 2010.
 This Journal is © 2010 The Royal Society
This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Appendix A
A.1. Definition of Hilbert space ℋ
This is the L^{2} Hilbert space of ℝ^{n}valued functions U(t) = (U_{1}(t), … ,U_{n}(t)), U^{′}(t) = (U^{′}_{1}(t), … ,U^{′}_{n}(t)), 0 ≤ t ≤ T, with inner product and norm given by U^{2}_{L2} = 〈U,U 〉_{L2}.
A.2. Derivation of equation (3.3)
This was presented in Rand [6]. Let t = ϕ_{m}(k) be the time when the concentration of the mth component, g_{m}, is at a maximum or a minimum value. Hence, ϕ_{m}(k) satisfies ġ_{m}(ϕ_{m}(k),k) = 0. Differentiating both sides with respect to k_{j} and rewriting,