A framework to understand the role of biological time in responses to fluctuating climate drivers

Understanding biological responses to environmental fluctuations (e.g. heatwaves) is a critical goal in ecology. Biological responses (e.g. survival) are usually measured with respect to different time reference frames, i.e. at specific chronological times (e.g. at specific dates) or biological times (e.g. at reproduction). Measuring responses on the biological frame is central to understand how environmental fluctuation modifies fitness and population persistence. We use a framework, based on partial differential equations (PDEs) to explore how responses to the time scale and magnitude of fluctuations in environmental variables (= drivers) depend on the choice of reference frame. The PDEs and simulations enabled us to identify different components, responsible for the phenological and eco-physiological effects of each driver on the response. The PDEs also highlight the conditions when the choice of reference frame affects the sensitivity of the response to a driver and the type of join effect of two drivers (additive or interactive) on the response. Experiments highlighted the importance of studying how environmental fluctuations affect biological time keeping mechanisms, to develop mechanistic models. Our main result, that the effect of the environmental fluctuations on the response depends on the scale used to measure time, applies to both field and laboratory conditions. In addition, our approach, applied to experimental conditions, can helps us quantify how biological time mediates the response of organisms to environmental fluctuations.

One of the biggest challenges faced by humanity is climate change [1][2][3] . A key characteristic of current climate change is the presence of extreme climatic events, i.e. strong fluctuating environmental conditions, manifested as storms, hurricanes and heatwaves. Hence, a challenge for ecologists consists in quantifying and predicting the responses of ecological systems (from populations to ecosystems) to such fluctuations [4][5][6][7] . Despite the emphasis in understanding ecology and evolution in fluctuating environments [8][9][10][11] most experiments concerning climate driven environmental variables focus on responses to constant conditions [11][12][13] . The currently growing body of work tackling responses to fluctuating environments highlights the complications in attempting to incorporate the role of variation in such environmental drivers [13][14][15][16][17] especially in dealing with their combined actions 18,19 . For instance, research on extreme events has identified five primary traits characterising heatwaves 20 , which should be added to the effect produced by the average condition experienced during that event.
An additional layer of complexity is given by the biological time (e.g. time to metamorphosis, to reproduction, generation times), characterising biological systems 21 . Biological time is governed by the interaction between the environmental and physiological processes and plays a critical role in driving fitness, population dynamics and community structure [21][22][23][24][25] . Biological time is critical because of three reasons. First, responses (e.g. survival rates) to environmental conditions are driven by a number of developmental processes operating at different biological time scales [26][27][28][29][30][31] , ranging from short (e.g. physiological acclimation) to medium (developmental plasticity) and long term (transgenerational plasticity, changes in gene frequencies). Second, if we characterise a fluctuation by its time scale, then the response will depend on the biological time characterising that species. From the perspective of organisms (e.g. with time scales ranging from those of from bacteria to trees), whether a fluctuation is long or short, is not determined by the chronological (= clock) time, but instead by its characteristic biological time 31 . Third, responses may be measured after a predetermined chronological time scale (e.g. at a given time in the year) or biological time scale (e.g. at maturity). The expression of biological responses in chronological www.nature.com/scientificreports/ time is obviously important to understand long term changes in seasonal habitats. However, measurements on biological time scales are critical for understanding population dynamics and evolution because such number drives individual fitness and the population growth rate. This is especially important in species that experience habitat shifts (e.g. most bottom living marine invertebrates, insects, migratory fish, birds and mammals). In such species, the biological time is "reset" at critical stages (e.g. metamorphosis, onset of migration) because, after the habitat shift, organisms experience the environmental conditions of a new habitat and the conditions in the old habitat might become irrelevant. Overall, understanding and quantifying the actual response to climate driven fluctuations requires that we also understand how responses are modified with a change in the time reference frame (from chronological to biological). We propose a framework to understand and quantify responses to fluctuations in one or more climate-driven environmental variables, considering biological time. We consider environmental fluctuations characterised by a magnitude (e.g. the amplitude) and a timescale (e.g. the period of the fluctuation or the time of exposure to a given magnitude: Fig. 1). The response is quantified at least once after the environmental drivers are experienced, with respect to chronological time or at a given life history event (e.g. at maturity). We use an experimental case and simulations to understand how biological responses to fluctuating environmental drivers are modified by the clock and biological time scale used to study the response. We structure this article as follows: First, we introduce definitions and a system of equations describing the biological responses. Second, the equations are explored using four specific cases. Third, in the methods section, we describe the experiments carried out to illustrate case 3.

Results
Mathematical theory. We consider a biological response (e.g. body size, survival, biodiversity) to two environmental drivers (i.e. any abiotic or biotic factor) but the same idea may be applied to a larger number of drivers. The response depends of a set of predictors consisting in the magnitudes (m 1 and m 2 ) and time scales of fluctuation of two drivers (i = 1, 2); in addition, the response is quantified at least once after the fluctuations have been experienced (Fig. 1a).
Time is defined using two different frames; chronological (= clock) time (measured by clocks) and biological time. For the "clock" time scales of the fluctuations (t 1 , t 2 ) there are associated biological times (τ 1 , τ 2 ). Likewise, for the clock time at which the response is quantified ( t * ) there is an associated biological time (τ * ). Biological time is the proportion of (clock) time needed to reach a life history event (e.g. moulting, maturity). Hence, for t 1 , t 2 and t * we obtain τ i = t i /D and τ * = t * /D, (D = clock time needed to reach such life history event). We express the τ i and τ * in terms of a function L = 1/D. For instance, for t * we obtain: where L = L(ω) = D −1 (ω) characterises the timing of a life history event (with units as the inverse of clock time units). L depends on the set of predictors ω associated to the fluctuations; an important set of predictors will be defined by thermal fluctuations (the amplitude and time scales), which in ectotherm species have a strong influence on developmental time 32,33 . We find by differentiation that L provides the transform function between clock and biological time frames; for instance, if L does not depend on any t i we have L = dτ/dt i .
The response is expressed as a function of the predictors defined above, as R(m 1, m 2, t 1 , t 2, t * ) = r[m 1, m 2, τ 1 (t 1 ), τ 2 (t 2 ) , τ * ]. The contribution of each predictor to the response is better understood by the partial derivatives  g. an organism) is exposed to two fluctuating environmental drivers (E 1 , E 2 ). The fluctuations are characterised by predictors, i.e. the amplitudes (m1, m2) and time scales (t 1 , t 2 ). Measurement of the response, R, are made at different times ( t * 1 ) after the fluctuations occurred (black circle); additional measurements may be carried out at other fixed times after t * 1 (grey circles). (B) In a factorial experiment the process would be repeated so that observations are made for a minimum of two levels per predictor giving 16  www.nature.com/scientificreports/ with respect to each predictor; this defines a system of partial differential equations (PDE; Supplementary note 1) which expressed in matrix form give the following matrix equation.
In the PDE (Eq. 2), the left-hand side is a vector column of the derivatives of the response in clock time (R), with respect to each predictor; the right-hand side is the standard (= inner) product of a matrix (M) by a vector of the derivatives of the response in biological time (r), i.e. R = Mr. The matrix contains the derivatives of the predictors with respect to each other, with time both expressed in clock or biological scales; one can think of M as an object containing coefficients that transform r into R in the same way as a constant (= 1000) would transform kilometres into meters of distance. The large number of terms in M highlights the considerable diversity and the challenges in quantifying responses to multivariate environmental fluctuations. We show below how to use Eq. (2) to quantify the effect of fluctuating environmental drivers on biological responses, as mediated by biological time.
First, we note that M contains three groups of terms: (1) Terms accounting for situations where the magnitude of a driver affects the magnitude of the second driver (e.g. temperature drives oxygen concentration in aquatic habitats): these are dm i /dm j for any i, j = 1, 2.  (1) and (2) are zero when they are mutually independent, such as in a factorial experiment with orthogonal manipulation. We will set those to zero in the rest of this analysis.
Second, we note that for group (3) there are three scenarios: (3a) biological time does not depend on any environmental driver. This is the trivial case where biological time is proportional to clock time, not considered here; M is simplified to a diagonal matrix, i.e. with constants in the diagonal, and zero's otherwise leading to a single constant term per equation (3b). Biological time depends on the magnitudes of any or both drivers. In such case, τ 1 τ 2 , and τ * will be driven by the same equation: if τ i = t i · L (m 1 , m 2 ) we obtain dτ i /dt j = dτ i /dt i = L (m 1 , m 2 ). (3c) Biological time depends on the time scale of the fluctuations: in such case, differentiating Eq. (1) with respect to time, we obtain dτ i /dt i = L + t i dL/dt i .
Here, we explore four special cases where the equations are simplified to highlight the importance of biological time in modifying the responses as compared to clock time. We start with the simplest case where there is a single environmental variable and then we consider cases with two variables. We focus on cases representing the most frequent experiments carried out on multiple driver research, i.e. factorial manipulations where terms of the groups 1 and 2 are zero. Case 1: responses to the magnitude of a single variable. We start with the simplest case i.e. where the response is driven by the magnitude of a single driver, e.g. temperature (= m). Examples of this case are laboratory experiments quantifying the effect of temperature on body mass or survival of a given species, or mesocosm experiments quantifying effects of temperature on species richness where thermal treatments are kept constant over time. Here, the response is quantified at different times, both in the clock and biological frames. In such case we have R(m, t * ) = r[m, τ * (m, t * )] and the PDEs simplify to.
From Eq. (3), and because dR/dm ≠ dr/dm, we see that the response to the magnitude of the driver depends on a component quantifying the effect biological time: as long as dτ * /dm ≠ 0 the time reference frame affects the observed effect of m on the response. The simulation illustrated in Fig. 2 shows a case where there are differences between the observed responses at clock vs biological times. In the simulated experiment, there is a strong effect of the magnitude of the driver on the response at clock time, but such effect is much less pronounced at biological time. By contrast, there is no effect when the response is measured in the biological time frame.
Equation (3) (details in Supplementary code 1) captures an obvious but important feature of experiments manipulating temperature over the development of ectotherms, for instance, from birth to metamorphosis; namely that there is no consistent definition of a simultaneous event across the different time frames. Experiments are usually stopped at different clock times because organisms must be sampled at the same biological time. All points located in the horizonal line in Fig. 3 represent simultaneous events, as defined in clock time occurring at different temperatures (e.g. whether an animal is dead or alive); however, simultaneous events occurring in biological time are represented by the points on the curve. Hence, Fig. 2 gives a geometric representation of such fact. Temperature as a driver of developmental rates 32 is a central candidate to produce responses that differ at clock vs biological time.
We explore further this case with an example where the response is expressed as a function of time and an instantaneous rate μ(m) quantifying for instance mortality, growth or biomass loss. For this example, we obtain R(m, t * ) = r[μ(m), τ * (m, t * )]. By differentiating in both sides, we get:  (4) shows that m affects the response through two components: the instantaneous rate (dμ/dm) and the biological time (dτ * /dm). We call the first component "eco-physiological" and the second component "phenological" (m drives the timing of a biological event, e.g. time to maturation). Those components are not evident if the response is expressed in clock time; otherwise we would obtain dR/dm = ∂R/∂μ · dμ/dm.
In order to better understand Eq. (4), consider an example where the response is biomass loss experienced by an organism during the process of migration (e.g. towards a feeding or reproductive ground); when the access to food during migration is very limited the result should be a decrease in body mass reserves through time. Let biomass (B) be modelled as an exponential decaying function of time and an instantaneous rate of biomass loss μ; let μ depend on temperature (= m) such that, μ = μ(m). In such case we obtain: By differentiation in both sides of Eq. (5) we get: Equation (6) shows the eco-physiological (dμ/dm) and phenological components (dτ * /dm) within the brackets. If μ responds linearly to temperature, then dμ/dm would be represented by a constant quantifying the thermal sensitivity of biomass loss; the value of such constant would depend on physiological processes associated to use of reserves to sustain movement and the basal metabolic rate. Likewise, if τ * responds linearly to temperature, the dτ * /dm would be driven by a constant controlling the sensitivity of developmental time to temperature. www.nature.com/scientificreports/ Because biomass is a trait that is central to fitness, Eq. (6) gives the indirect contribution of phenological and physiological responses to fitness. Assuming that fitness should be maximised, adaptive responses should involve the mitigation of negative effect of m on both components of Eq. (5), represented by the partial derivative of the right-hand term. For instance, organisms with the ability to minimise the eco-physiological effect (through e.g. a compensatory physiological mechanisms) or the phenological effect (e.g. shortening the exposure time) would complete the migration minimal loss of reserves.
By generalization, Eqs. (4-6) help us to provide biological meaning to the terms of the matrix M: any term of the form dτ * /dm j , dτ i /dm j or dτ i /dt j represents the effect of an environmental driver on the timing of a phenological event; hence, they are phenological components. Terms that contain the effect of an environmental variable on an instantaneous rate are eco-physiological components. By substitution we find that the terms of the matrix in Eq. (2) can be classified in two categories according to whether the component is eco-physiological (E) or phenological (P): Case 2: multiple driver responses. Here we expand the previous case by looking at a response to the magnitude of two different drivers; i.e. keeping the levels of each driver constant over the duration of the experiment. Examples of this case are experiments quantifying the effect of temperature and nutrient load on body mass (e.g. in a rearing containers) or species richness (e.g. in mesocosms). This case is represented by the terms of first two rows of the matrix and the vectors of Eq. (2), with the terms of the remaining rows set to zero. Here, there are different scenarios, but we focus on the one highlighting the importance of biological time.
Consider a case where biological time depends on the magnitude of the first driver while the response is explicitly driven by the magnitude of the second driver (Fig. 4). For instance, the response may be the survival rate of a host organism exposed to different temperature and parasitic load. The response in clock time is described as R(m P , t * ). The driver controlling the biological time is temperature (m T ) while the parasitic load (m P ) controls survival. In such case, dτ * /∂m P = 0, dR/dm P ≠ 0 and dR/dm T = 0. Although by definition the response in clock time does not depend on m T , it will do so in biological time. This is because, applying the matrix multiplication in Eq. (2), we obtain: The second right-hand term in Eq. (8a) quantifies the effect of temperature on the response mediated by biological time. In order to better understand the responses, consider a simple linear response: R = R 0 − m P ·t * and notice that, for a fixed clock time (t * c ) the effect of the magnitude of parasitism is constant (dR/dm P = −t * c ); hence, the response can be understood, geometrically, as a flat surface with slope not depending on temperature. Now, note that under the specific conditions of our example, r = R 0 − m P ·τ * /L(m T ). Hence, for a fixed biological time (τ * c ) we obtain ∂r/∂m P = −τ * c /L(m T ); i.e. the importance of the parasitic effect depends now on temperature. In addition, this example is valid for the case of additive effects of any two environmental drivers: assuming R = R 0 − (a 1 ·m P + a 2 · m T )·t * (a 1 , a 2 are constants), we obtain dR/dm P = −a 1 t * ; however, ∂r/∂m P = −a 1 τ * c /L(m T ). In words, additive effects observed in clock time become interactive in biological time. This is illustrated in the simulation (Supplementary code 2) depicted in Fig. 4: the response in clock time depends on a single driver (parasite load); however, the response in biological time is interactive, i.e. the effect of parasite load depends on temperature. We performed the so-called point-of-reserve-saturation experiment (PRS 35 ), i.e. exposing groups of recently hatched larvae of the crab Hemigrapsus sanguineus to different initial feeding periods (= our time scale of fluctuation), after which they were starved until they either died or moulted to the second larval stage ( Supplementary  Fig. 1). H. sanguineus is originated from East Asia but has invaded the Atlantic shores of North America and North Europe 36,37 . This experiment was carried out at 4 temperature levels (15-21 °C), within the range of thermal tolerance of larvae of this species, i.e. where the magnitude of temperature does not affect survival 38,39 . In addition, because there is a single level of food limitation (= starvation), the magnitude of food limitation (m F ) is not considered as a variable in the example.
The response variable was the proportion of first stage larvae surviving the moulting event to the second stage, set to biological time τ * = 1. In response to different starvation periods (preceded by feeding), the survival shows a sigmoidal pattern 35 , characterised by a parameter, PRS 50 . This is the point of development where larval reserves are "saturated"; i.e. enough reserves have been accumulated during the previous feeding period to ensure survival and moulting to the next stage.
Under the conditions of the experiment, the survival proportion (= R) is driven only by the time scale of a fluctuation (here t 1 = t, τ 1 = τ for simplicity), characterised by the starvation period; hence, given that there is a single time of observation fixed to τ * = 1. Because biological time does not depend t, we get L = dτ/dt and: Equation (9) is represented in the PDE by the terms of row 3 and column 4 of M multiplied by the term of row 3 of the column vector r; dτ/dt = L(m), m represents the magnitude of temperature.
The relationship between biological time and temperature was best explained by a power function D(T) = aT b (Fig. 4A, Supplementary Table 1, Supplementary Fig. 2), in consistence with previous studies 36,40 . The interaction between starvation time and temperature was weak ( Supplementary Fig. 3); best models retained starvation time only at 21 °C where the percentage of explained variance was still low (R 2 < 0.2). The full range of starvation times resulted in a variation of developmental time of < 2 days, while the full range of temperature used resulted in variations of 8 days (range 5-14 days); hence, we approximated the model as L depending on temperature as L = 1/(aT b ).
Survival showed an S-shape pattern consistent with results found for other species 35 . When the starvation time was expressed in clock time (PRS 50 = t 50 in days) there was a dilation/contraction effect of the response   Table 4); the range of percent values found here is also consistent with findings in other species (40-60%) 35 . There were therefore important differences in the effect of temperature on estimates of PRS 50 depending on the choice between time scale (Supplementary Table 5). In synthesis, in the biological time scale we found a simple function showing that the PRS 50 was less responsive to a change in temperature than in clock time; we will address this point in the discussion in the context of physiological time keeping mechanisms. We now consider a response interpreted as a decay in performance of an organism, where longer time scales of the fluctuation increase the biological time, as expected for cases where organisms are exposed to suboptimal conditions (e.g. food limitation experienced over a given time scale). There are obviously many possible scenarios but for better understanding, we consider Cases 4A-D, where dL/dt is negative, reducing performance and where the functions linking m and t with L act additively or multiplicatively. Additive responses are will be illustrated with L = [k 1 /m + k 2 /t], while multiplicative responses will be illustrated L = k 3 /mt, with k 1 , k 2 and k 3 as constants.
The first two cases focus on Eq. (10a), which may be considered as an extension of Case 3. Case 4A: additive response: in such case dL/dt depends only on t and we get that dR/dt and dr/dτ differ by a factor k 1 /m: In addition, if L only depends on t, the result is that the response in clock time does not depend on the time scale of fluctuation (dR/dt = 0) because the terms of Eq. (10a) associated to L cancel out. However, the response in biological time does not need to be zero (dr/dτ may not be zero). For instance, assume that r(m, τ) = exp(−ω·m·τ), and L = k 2 /t , with ω = constant. We obtain R = exp(−ω·m·k 2 ) and dR/dt = 0, while dr/dτ = −ω·m·exp(−ω·m·τ).
Case 4B: Multiplicative response of m and t in L. In such case, dR/dt = 0 because the terms associated to L cancel out, while dr/dτ may not be zero.
The next two cases focus on Eq. (10b), which is an extension of cases 1 and 2. The effect of the time scale of the fluctuation depends again on how it relates, through L, to the effect of the magnitude of the fluctuation. We interpret the response as a decay in performance (or fitness), contributed by the "ecophysiological" and "phenological" components: In Eq. (11), the phenological component (within the brackets) is driven by two terms. The expression t·dL/ dm shows that the time scale of the fluctuation act as increasing the exposure to the suboptimal conditions and contributes to a further reduction in performance. Equation (11) takes different forms depending on whether the functions linking t and m with L are additive or multiplicative.
Case 4C: when the response is additive, the phenological contribution is proportional to the time scale of the fluctuation which then contributes to a further reduction of the performance: Case 4D. When the response is multiplicative, the phenological component responds non-linearly to the time scale of the fluctuation.
In both 4C and 4D the effect of the magnitude, m, on the response depends on whether the fluctuation is quantified in terms of clock or biological time. For example, when L = k 1 /m or L = k 3 /mt, dR/dm does not depend on t while dr/dm depends on τ. Overall, Cases 4A-D further extend the relevance of cases 1 and 2 in understanding the effects of fluctuating environmental drivers on the responses.

Discussion and conclusions
We have introduced a mathematical framework to better understand and quantify the responses of biological systems to fluctuations in climate driven stressors. Central to that framework is the need to consider biological time as playing a role in driving ecological and evolutionary processes 20 . We did not consider the average magnitude of the fluctuation because responses to fluctuations can differ from responses to the average in two different ways. First, under extreme fluctuations critical tolerance levels may be surpassed leading to the collapse or irreversible shift of the biological system under study. The most obvious case is when the temperature surpasses critical tolerance levels leading to death 42 but experiencing less extreme levels may lead to negative carry-over effects or acclimatory responses 15 , not present when average conditions are experienced. Second, when the relationships between environmental drivers (e.g. temperature) and biological responses (e.g. physiological performance) is non-linear, average conditions do not predict well the expected response to a fluctuating environment 13,17 .
We started with the simple Case 1 of a single variable and no fluctuation. This case may be considered as trivial but helped as a step to understand more complex cases. For instance, through Case 1 we noted that when the matrix M of the PDE is not the identity matrix, simultaneous events in biological time are not so in clock time; this known observation in experimental research was then identified geometrically as the difference between a curve and straight lines when plotting the response in the space defined the magnitude of an environmental variable and the time scale of observation (Fig. 2). With further analysis of Case 1, we noted that the M contained two types of terms interpreted as eco-physiological and phenological effects. The specific example of the migrating organism provided further biological interpretation to the components of the differential equation. We focused on the negative effects of temperature (biomass loss) and then noted that adaptations should minimise either one or both the phenological and the physiological components if biomass loss were to be minimised. This is for instance the case of the evolution of early life histories of marine invertebrates to habitats characterised by limited food availability 43 . Where food is available, marine invertebrate tend to develop through feeding larval stages; however, where food is too limiting, most species develop through non-feeding larvae with abbreviated larval phase. In such case, the allocation of maternal reserves into eggs contributes to minimise both the ecophysiological a phenological components with respect to survival, because both the mortality rates and developmental time are independent of food availability.
In Case 2, we introduced the magnitude of a second variable, to explore more complex scenarios studied through factorial experiments carried out under constant conditions; i.e. not yet considering the time scale of a fluctuation. Here, the presence of an interactive response depended on the scale used to measure time. This finding is central to climate change biology, given the interest on interactive effects of multiple environmental variables on biological systems 12,18,19,44 . A critical question is whether effects are additive or whether they are antagonistic or synergistic. Additive effects refer to situation where the response can be modelled from the isolated effect of each single variable; many biological responses are however synergistic or antagonisitic 44 . Synergistic responses occur when the combined effects are larger than the expected contribution of each separate variable. Synergistic responses are critical when they are negative, for instance the combined effect of habitat loss and an environmental stressor, as they can drive ecosystem collapse. By contrast, antagonistic responses imply a mitigation effect. For management, it is essential to get the response right because resources for actions are limited and disrupting synergies may be considered a priority 51 . In such context our findings suggest that management depends on using the correct time frame to measure the response. Interactive effects also change when responses at low levels of organization (e.g. consumer resource functional responses) are used to predict those at higher levels (e.g. population dynamics 45 ), because of the non-linear nature of the function mapping the response across levels. In our case there is a non-linearity in that the components of M are partial derivatives which depend on the predictors (i.e. the original functions are non-linear).
We used Case 3 to explore responses including the time scale of fluctuation and to better understand how experimental results are interpreted in the light of the PDE. We studied responses in larval stages because of the relevance of marine larvae in driving climate-change effects on marine organisms: most marine organisms (e.g. mussels, crabs, fish) develop through a pelagic larval stage, and larval dynamics affect species range 38,39 , population dynamics 46 , connectivity 47 and community structure 48 . Warming provides a new context where larvae need to cope with fluctuations in e.g. food abundance (or other variables) in a scenario of increased metabolic demands due to higher temperatures 49 .
Case 3 highlighted the importance of understanding how temperature drives time keeping mechanisms in biological systems, understood as those responsible for setting the pace and regulating the timing of life history events 50 . For Case 3, the phenological component of the PDEs captured the effect of temperature on PRS 50 which instead reflects hormonal control of the so called "D 0 -threshold". This threshold is surpassed when moulting hormones are triggered and the premoult period starts 51 ; after D 0 , development proceeds at a rate that is independent of food levels and larvae will moult. Case 3 therefore highlights the importance of understanding how temperature drives the hormonal regulation of moulting, for the formulation of mechanistic models predicting survival. By extension, knowledge of role of hormones and other signalling mechanisms 50,52,53 should help the formulation of models in other species. In cases where organisms undergo acclimation, an important question is how the time scale of acclimation relates to time keeping mechanisms, including hormonal processes and metabolic rates 31 . Acclimation speed correlates negatively with body size, most likely driven by the positive effect of metabolic rates on acclimation speed 54 . Perhaps acclimation time, as a fraction of developmental time, varies little with temperature or alternatively, acclimation time and the developmental processes setting phenological events have different sensitivities to temperature.
In Case 4, we introduced effects of the time scale of a fluctuation on the function L and found equations that may be considered extensions of the previous cases. A critical question is to determine situations when responses may follow Cases 1-3, 4 or be further simplified into the trivial case. An important point is therefore to determine  Supplementary Fig. 4). However, whether the developmental time is driven by time scale of the fluctuation depends on the timing of the fluctuation in relation to size thresholds reached as organisms grow [55][56][57] . Our work suggest that we need to understand how such thresholds relate to time-keeping mechanisms. We used the PDEs to understand the importance of the choice of time scale, within an experimental setting. In addition, Case 3 suggest the set of conditions where simple models predict responses to environmental fluctuations. In the example, the response may be approximated by a function quantifying the relationship between temperature and biological time and a second function controlling the timing of the starvation period. Notice that under the range of temperatures considered, the transformation from clock to biological time led to a simple model with high predicting capacity (> 90% of explained variation) although it ignores the significant (but small) effects of food limitation on developmental time (Supplementary Fig. 3). In similar cases, data of the response from a narrow temperature range would give an approximation of ∂r/∂τ when scaled in biological time.
In that case, additional data on the effect of temperature on biological time may be used predict the response. An important point is the set of conditions where equation-8 may be used with safety. For our experimental system, we hypothesise that equation-2 had a good fit because the developmental time varied little with the time of starvation and the range of temperatures was within the so called "pejus range" 58 i.e. where survival was high irrespective of the magnitude of temperature (dR/dm T ≈0). However, assumptions are not valid if the response fall within Case 4 or at temperatures beyond the pejus threshold where a change in temperature have strong effects on the response.
Under the conditions of equation-8, one may combine mechanistic sub-models as modules (e.g. for each of the partial derivatives). In Case 3, the first module is given by L which may be modelled from metabolic theories (MTE 32,33 ). For instance, in the MTE, the effect of temperature on biological time, is represented by the Arrhenius equation, which instead will determine L. For the second module (represented by ∂r/∂τ), we can associate the response to hormonal control of development which drive the timing the switch of the sigmoid survival function. More in general, sigmoid responses are characteristic of populations or ecosystems exhibiting regime or phase shifts 59,60 . At the population level, phase shifts reflect an unstable equilibrium (saddle points) point driven by thresholds associated to density-dependent changes in mating and reproduction. Hence, the mechanisms associated to such phase shifts would be captured in the response function expressed on biological time scale (e.g. the generation time for population level responses).
There are two important points concerning relevance of our approach to characterizing biological responses to environmental fluctuations in the field. First, our main finding, i.e. that responses to environmental fluctuations depends on the scale used to measure time, is valid for both field and experimental conditions. The use of the experimental set up only facilitates teasing apart the independent contributions of each of the predictors (i.e. the magnitudes and time scales of fluctuations). However, whether one can directly apply the equations straight away to field conditions, depends on meeting the assumptions used to formulate the equations; in this sense our approach is not different from any other experimental approach and the usual recommendations apply 18 . There are three assumptions: (1) the right few predictors are identified (e.g. magnitudes and times scales of temperature and any other factor). (2) No covariation among predictors; (3) fluctuations characterised by well-defined values of predictors. A first challenge in field applications is the complexity shown by natural environmental fluctuations. For instance, real fluctuations (e.g. heatwaves) consist of a sequence of oscillations; in addition, fluctuations may be characterised by descriptors other than the period and amplitude (e.g. the rate of daily temperature increase in tropical habitats) 20 . A second challenge is that environmental variables covary in the field 12 . This include situations, not considered here, where different environmental variables fluctuate sequentially 12,31 , and an additional time scale must be included in the PDEs (i.e. the one separating the fluctuations). Third, in the field, fluctuations may be characterised by means and variances of the predictor values and attempts to model the average biological response need to consider issues associated to non-linear responses 13,17 . However, the most likely scenario is that field observations inform the design of future experiments. For example, field studies can identify the main environmental variables, the most important traits characterising the fluctuations, whether fluctuations of different variables occur sequentially or simultaneously 12,18 .
Another important question is what time frame should we choose. The selection of the appropriate reference frame will depend on the question asked by the researcher. In experiments aimed at determining the effects of environmental fluctuations on body size or survival at e.g. maturation, biological time will be the obvious choice. Biological time will be a choice in situations where organisms experience habitat shifts through the life cycle, for example in species where the larval habitat is aquatic and the adult habitat is terrestrial. In such case, once the aquatic larvae metamorphose to a terrestrial juvenile stage, the importance of the larval habitat conditions for the survival of the juvenile is likely to be low; hence, what matters are the conditions experienced as larvae up to the time of metamorphosis. There are however cases where the decision is less clear. For example, where larvae and adults coexist or where key environmental variables (e.g. weather conditions) can affect both the larval and adult habitat. Another example is the one given by, mesocosm experiments, used to study the effect of warming on populations and communities 18 ; under warming, it is likely that populations fluctuate over more generations that in the absence of warming. If the number of generations is important, it will be helpful to analyse the responses considering both clock and biological time (i.e. after a fixed number of generations). The same experiment may be designed to understand the importance of a fluctuation occurring at a characteristic clock time scale, associated to e.g. seasonal fluctuations or 5-day long heatwaves; in such case, the scale of the fluctuation should be kept in clock time. Overall, we will profit from analysing the response in both time frames, as they will provide different pieces of information.

Methods
Experiments were carried out in automated fully programmable incubators (RUMED-EcoLine R ). At each combination of temperature and starvation period, three replicate groups (10 larvae each) were kept in well oxygenated and filtered natural seawater (mesh = 1 μm, salinity = 32.5 PSU) in glass vials of 100 ml. Water and food were changed every day (except at 21 °C where food was checked every ~ 12 h, at 7:00 and 18:00 h); at such times larvae were checked for moulting or mortality (dead larvae were removed from cultures). Larvae were fed freshly hatched Artemia sp nauplii, provided at libitum (density 5 nauplii per ml) during the feeding periods.
The effect of temperature (T) or feeding period on developmental time (D) were evaluated through model selection 61,62 using general least squares for model fitting and Akaike information criterium (AIC) for model comparison. Analyses were carried out in R using the package nlme 61 . Four models were compared, i.e. linear ( D = −a · T+b), exponential ( D = a · e −b·T ), power ( D = a · T −b ) and Arrhenius ( D = a · e b (T+273 ), where a and b are constants. Models were fitted after appropriate transformations, log(D) for exponential, log(D)-log(T) for power and log(D) vs 1/(T + 273) for Arrhenius model.
Effects of initial feeding periods on survival were evaluated using the sigmoidal dose response function: where f m and f M are the asymptotic minima and maxima respectively, k is the slope parameter and PRS 50 is the timing of the inflection point where f(PRS 50 ) = f(x M )/2. Model fitting was carried out by non-linear regression (in GraphPad Prism software), with feeding period expressed in both chronological and biological time scales. When the biological time scale was used, a single model was fitted to data from all temperatures. In that case, there were sufficient degrees of freedom to enable appropriate estimation of the four model parameters. In addition, separate models were fitted by temperature using at biological and chronological time. For those models, the estimation of some parameters was unreliable (e.g. extremely large confidence intervals); hence, we focused on estimating PRS 50 , which is the parameter determining the time of the inflexion point in the curve. Therefore, we set f m = 0 and f M = M, with M being the average survival of the last three points. In addition, for 19 and 21 °C we set k to a constant value obtained after our initial attempt to estimate k.

Data availability
Data will be available in the portal PANGEA.