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 (m1 and m2) 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 (t1, t2) 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 t1, t2 and t\(^*\) we obtain τi = ti/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 time32,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 ti we have L = dτ/dti.
The response is expressed as a function of the predictors defined above, as R(m1, m2, t1, t2, t\(^*\)) = r[m1, m2, τ1 (t1), τ2(t2), τ\(^*\)]. The contribution of each predictor to the response is better understood by the partial derivatives 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.
$$\left[\beginarrayc\fracdRdm_1\\ \fracdRdm_2\\ \fracdRdt_1\\ \fracdRdt_2\\ \fracdRdt^^*\endarray\right]=\left[\beginarrayccccc1& \fracdm_2dm_1& \fracd\tau _1dm_1& \fracd\tau _2dm_1& \fracd\tau ^^*dm_1\\ \fracdm_1dm_2& 1& \fracd\tau _1dm_2& \fracd\tau _2dm_2& \fracd\tau ^^*dm_2\\ \fracdm_1dt_1& \fracdm_2dt_1& \fracd\tau _1dt_1& \fracd\tau _2dt_1& 0\\ \fracdm_1dt_2& \fracdm_2dt_2& \fracd\tau _1dt_2& \fracd\tau _2dt_2& 0\\ 0& 0& 0& 0& \fracd\tau ^^*dt^^*\endarray\right]\cdot \left[\beginarrayc\fracdrdm_1\\ \fracdrdm_2\\ \fracdrd\tau _1\\ \fracdrd\tau _2\\ \fracdrd\tau ^^*\endarray\right]$$
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 dmi/dmj for any i, j = 1, 2. (2) Terms accounting for cases where the magnitudes and time scales of stressors are related: dmi/dtj and dmi/dti. (3) Terms where biological time depends on the magnitude or time scale of the environmental fluctuation dτi/dtj and dτi/dmj. The terms of groups (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 = ti · L (m1, m2) we obtain dτi/dtj = dτi/dti = L (m1, m2). (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/dti = L + ti dL/dti.
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.
$$\fracdRdm=\frac\partial r\partial m+\frac\partial r\partial \tau ^^*\cdot \fracd\tau ^^*dm$$
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 rates32 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:
$$\fracdRdm=\frac\partial r\partial \mu \cdot \fracd\mu dm+\frac\partial r\partial \tau ^^*\cdot \fracd\tau ^^*dm$$
Equation (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:
$$B(m,t)=e^-\mu \left(m\right)\cdot t^^*=e^-\mu \left(m\right)\cdot \tau ^^*\left(m,t^^*\right)$$
By differentiation in both sides of Eq. (5) we get:
$$\fracdBdm=-e^-\mu \left(m\right)\cdot \tau ^^*\left(m,t^^*\right)\left\\tau ^^*\cdot \fracd\mu dm+\mu \cdot \fracd\tau ^^*dm\right\$$
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.
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τ\(^*\)/dmj, dτi/dmj or dτi/dtj 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):
$$\left[\beginarraycccccE& 0& P& P& P\\ 0& E& P& P& P\\ 0& 0& P& P& 0\\ 0& 0& P& P& 0\\ 0& 0& 0& 0& P\endarray\right]$$
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(mP, t\(^*\)). The driver controlling the biological time is temperature (mT) while the parasitic load (mP) controls survival. In such case, dτ\(^*\)/∂mP = 0, dR/dmP ≠ 0 and dR/dmT = 0. Although by definition the response in clock time does not depend on mT , it will do so in biological time. This is because, applying the matrix multiplication in Eq. (2), we obtain:
$$\frac\partial R\partial m_T=\frac\partial r\partial m_T+\frac\partial r\partial \tau *\cdot \fracd\tau *dm_T$$
$$0=\frac\partial r\partial m_T+\frac\partial r\partial \tau *\cdot \fracd\tau *dm_T$$
$$\frac\partial r\partial m_T=-\frac\partial r\partial \tau *\cdot \fracd\tau *dm_T$$
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 = R0 − mP·t\(^*\) and notice that, for a fixed clock time (t\(^*\)c) the effect of the magnitude of parasitism is constant (dR/dmP = −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 = R0 − mP·τ\(^*\)/L(mT). Hence, for a fixed biological time (τ\(^*\)c) we obtain ∂r/∂mP = −τ\(^*\)c/L(mT); 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 = R0 − (a1·mP + a2· mT)·t\(^*\) (a1, a2 are constants), we obtain dR/dmP = −a1t\(^*\); however, ∂r/∂mP = −a1τ\(^*\)c/L(mT). 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.
Case 3: role of clock and biological time scale of fluctuation
Previous examples did not consider, the time scale of the fluctuations as drivers of the response. Here we explore how a biological variable (= survival rate) responds to different levels of magnitude of a driver (= temperature) and to simultaneously changing the time scale of a fluctuation (from clock to biological time) of a second driver (= food limitation). As model, we use larval stages of a crab because there is sufficient information on the effect of temperature and food levels on survival and the timing of moulting33,34.
We performed the so-called point-of-reserve-saturation experiment (PRS35), 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 Europe36,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 survival38,39. In addition, because there is a single level of food limitation (= starvation), the magnitude of food limitation (mF) 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 pattern35, characterised by a parameter, PRS50. 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 t1 = t, τ1 = τ for simplicity), characterised by the starvation period; hence, R = R(t) = r[τ(t)] 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:
$$\fracdRdt=\frac\partial r\partial \tau \cdot \mathcalL(m_2)$$
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) = aTb (Fig. 4A, Supplementary Table 1, Supplementary Fig. 2), in consistence with previous studies36,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 (R2 < 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/(aTb).
Survival showed an S-shape pattern consistent with results found for other species35. When the starvation time was expressed in clock time (PRS50 = t50 in days) there was a dilation/contraction effect of the response curve, quantified by the PRS50 and driven by the effect of temperature on biological time (Fig. 4B, Supplementary Table 2). When time was expressed in biological time units (PRS50 = τ50), a single response curve explained 91% of total variation (Fig. 4C, Supplementary Table 3), irrespective of temperature. The estimate of parameters by temperature showed PRS50 in the range of 0.47–0.58% of moulting time with a slight decrease towards higher temperatures (Supplementary 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 PRS50 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 PRS50 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.
Case 4: biological time depends on the time scale of the fluctuation
Here, we generalise the above cases by considering situations where both the magnitude and time scale of an environmental fluctuation drive biological time. For simplicity, we consider a single driver. In such case, L = L(m,t) and the PDEs reduce to:
$$\frac\partial R\partial t=\frac\partial r\partial \tau \cdot \mathcalL+\frac\partial r\partial \tau \cdot t\cdot \fracd\mathcalLdt$$
$$\frac\partial R\partial m=\frac\partial r\partial m+\frac\partial r\partial \tau \cdot t\cdot \fracd\mathcalLdm+\frac\partial r\partial \tau *\cdot t^^*\cdot \fracd\mathcalLdm$$
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 = [k1/m + k2/t], while multiplicative responses will be illustrated L = k3/mt, with k1, k2 and k3 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 k1/m:
$$\frac\partial R\partial t=\frac\partial r\partial \tau \cdot \frack_1m$$
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 = k2/t , with ω = constant. We obtain R = exp(−ω·m·k2) 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:
$$\fracdRdm=\frac\partial r\partial \mu \cdot \fracd\mu dm+\left\\frac\partial r\partial \tau \cdot t\cdot \frac\partial \mathcalL\partial m+\frac\partial r\partial \tau ^^*\cdot t^^*\cdot \frac\partial \mathcalL\partial m\right\$$
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:
$$\fracdRdm=\frac\partial r\partial \mu \cdot \fracd\mu dm-\frack_1m^2\left\\frac\partial r\partial \tau t+\frac\partial r\partial \tau ^^*t^^*\right\$$
Case 4D. When the response is multiplicative, the phenological component responds non-linearly to the time scale of the fluctuation.
$$\fracdRdm=\frac\partial r\partial \mu \cdot \fracd\mu dm-\frack_3m^2\left\\frac\partial r\partial \tau +\frac\partial r\partial \tau ^^*\cdot \fract^^*t\right\$$
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 = k1/m or L = k3/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.