Beyond Arrhenius: Nonlinear and negative temperature scaling of biological rates from multi-step mechanisms

Beyond Arrhenius: Nonlinear and negative temperature scaling of biological rates from multi-step mechanisms | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Beyond Arrhenius: Nonlinear and negative temperature scaling of biological rates from multi-step mechanisms Simen Jacobs , Federico Vazquez , View ORCID Profile Nikita Frolov , View ORCID Profile Lendert Gelens doi: https://doi.org/10.1101/2025.09.01.673554 Simen Jacobs 1 Laboratory of Dynamics in Biological Systems, Department of Cellular and Molecular Medicine, KU Leuven , Herestraat, 49, Leuven, Belgium Find this author on Google Scholar Find this author on PubMed Search for this author on this site Federico Vazquez 2 Instituto de Cálculo, FCEyN, Universidad de Buenos Aires and Conicet , Intendente Guiraldes 2160, Cero + Infinito, Buenos Aires C1428EGA, Argentina Find this author on Google Scholar Find this author on PubMed Search for this author on this site Nikita Frolov 1 Laboratory of Dynamics in Biological Systems, Department of Cellular and Molecular Medicine, KU Leuven , Herestraat, 49, Leuven, Belgium Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Nikita Frolov Lendert Gelens 1 Laboratory of Dynamics in Biological Systems, Department of Cellular and Molecular Medicine, KU Leuven , Herestraat, 49, Leuven, Belgium Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Lendert Gelens For correspondence: lendert.gelens{at}kuleuven.be Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Temperature shapes all biological processes, particularly during the early development of ectothermic organisms. A widely used framework for describing temperature dependence is the Arrhenius equation, which predicts an exponential increase in rates with temperature. However, biological rates often deviate from this prediction when measured across broader temperature ranges. While negative apparent activation energies are often attributed to protein denaturation, this cannot explain similar behavior observed at temperatures where enzymes remain stable. These broader scaling patterns remain mechanistically unexplained. Here we present a general Markov chain framework for modeling biological timing as cascades of reversible, temperature-dependent steps. Applying this model to 121 published datasets spanning diverse species, biological timescales, and temperature ranges, we find that a consistent three-zone scaling pattern emerges: Arrhenius-like behavior at low and high temperatures, and a quadratic exponential regime at intermediate temperatures. We show that this pattern arises naturally from differences in activation energies across steps in the network. The quadratic exponential regime is an emergent feature of averaging across many steps and is robust to variation across network realizations. In contrast, Arrhenius-like scaling at the extremes tends to be more variable and originates from smaller sub-networks. Apparent negative activation energies can emerge naturally from the dynamics of multi-step networks, even in the absence of protein denaturation. Our framework provides a unified mechanistic explanation for diverse temperature-scaling behaviors in biology and may help predict how developmental and physiological processes respond to environmental change. Although we focus on development, the model is broadly applicable to biological systems governed by multi-step reaction networks. Introduction Temperature influences all living organisms, shaping their physiology, behavior, and overall fitness. This is particularly critical during early embryonic development, where temperature acts as a key regulator of developmental rates, impacting biochemical reactions, metabolic activities, and organismal performance. While endothermic organisms maintain stable internal temperatures, ectotherms are especially vulnerable to environmental temperature changes, making them sensitive to the effects of climate change. Understanding how temperature governs biological processes is thus essential for predicting the consequences of global warming on biodiversity, particularly for species whose habitats and developmental success are tightly linked to thermal conditions ( 1 – 3 ). Biological processes exhibit complex temperature dependencies across scales, from enzymatic catalysis to organismal growth. These dependencies are often represented by thermal performance curves, which display optimal performance at specific temperatures and declines at extremes. Such curves provide insights into how temperature governs life processes, ranging from enzymatic activity to embryonic development ( 4 , 5 ). The Arrhenius equation ( 6 – 8 ) is widely used in biology to describe the temperature dependence of reaction rates. Initially derived for chemical reactions, it takes the form: where κ ( T ) is the rate constant, A is a pre-exponential factor, E the activation energy, R the universal gas constant, and T the temperature in Kelvin. This response function is typically plotted as ln κ ( T ) versus 1 /T , or using the associated time τ ( T ) = 1 /κ ( T ): which yields a simple linear relationship. The Arrhenius equation predicts a monotonic increase in rate, which holds reasonably well over moderate temperature ranges ( 9 , 10 ). If the Arrhenius equation were universally valid and processes had similar activation energies, a few degrees of warming would speed up rates in a predictable way. For example, processes would speed up by ≈23% (assuming Q 10 = 2, that is, the rate roughly doubles per 10 ° C). In reality, activation energies typically differ, so even under ideal Arrhenius kinetics different steps accelerate by differ-ent amounts and relative timings drift with temperature, potentially disrupting coordination 1 . Moreover, empirical observations show that even modest temperature increases often lead to pronounced deviations from Arrhenius behavior, especially at the upper end of the physiological temperature range. For instance, a temperature change of just a few degrees can trigger high fever in humans or completely arrest development in ectothermic embryos. In these non-Arrhenius regimes, biological rates frequently plateau, decline, or cease altogether. Such breakdowns arise from a range of factors, including protein denaturation, the activation of stress response pathways, and the nonlinear dynamics of complex biochemical networks ( 12 , 16 – 18 ). These phenomena underscore the limitations of the Arrhenius framework in fully capturing biological temperature responses, particularly near physiological boundaries or across broad thermal ranges. They also highlight the need for mechanistic models that can explain not only species-level differences, but also variability across different realizations of biochemical networks. To understand how such deviations appear throughout development, we examine temperature scaling across both early and later stages. Early embryonic cleavages are often driven by well-defined biochemical oscillators with minimal regulation, while later stages involve more complex and tightly regulated processes. Yet, as we show below, non-Arrhenius behavior occurs across this entire developmental spectrum, suggesting that it arises from general features of biological timing rather than stage-specific mechanisms. Figures 1A and B illustrate this distinction. Panel A shows a schematic of the early, rapid cleavage divisions in the X. laevis frog embryo, regulated by a periodic cell cycle oscillator ( 12 ). These divisions occur without gap phases, checkpoint controls, or transcription ( 19 , 20 ), allowing fast, regular cycles of proliferation. In contrast, Figure 1B sketches the slower and more variable timing of later developmental transitions, which involve diverse regulatory pathways, including gene expression, signaling, and morphogenesis. Download figure Open in new tab Fig. 1. Developmental timing scales non-Arrhenius. (A) Schematic of the first cleavages in the early X. laevis embryo (adapted from ( 11 )), driven by periodic cell cycle oscillations ( 12 ). (B) Schematic of selected later developmental stages during X. laevis development (adapted from ( 11 )), showing varying and longer durations between stages. (C) Period of the early embryonic cell cycle in X. laevis ( 12 ), plotted on a logarithmic scale as a function of inverse absolute temperature. The best fits using the quadratic exponential (QE, solid line) and Arrhenius (dotted line) models are shown. (D) Periods of the early embryonic cell cycle for five species— X. laevis, X. tropicalis, D. rerio ( 12 ), and C. briggsae, C. elegans ( 13 ), with best fits using either the double exponential (DE, dashed line) or QE (solid line) models. (E) Timings of four developmental stages in D. melanogaster : from the 14th cell cleavage to germband retraction, germband retraction to first breath ( 14 ), total embryonic duration, and total pupal duration ( 15 ). Since both DE and QE models fit the data equally well, only the simpler QE fit (solid line) is shown. Temperature scaling behavior across these developmental stages is illustrated in Figures 1C–E . In panel C, the early cell cycle period of Xenopus laevis is plotted on a logarithmic scale against inverse absolute temperature. The data appear approximately linear at lower temperatures, consistent with Arrhenius scaling, but deviate significantly at higher temper-atures. This upward curvature is more accurately described by the quadratic exponential (QE) model than by a single Ar-rhenius fit ( Table 1 ). Panel D extends this analysis to five species: X. laevis, X. tropicalis, D. rerio, C. briggsae , and C. elegans , showing the temperature dependence of early cleavage timings. Most species display Arrhenius-like scaling at low temperatures, but deviate at higher temperatures. These deviations are abrupt in C. elegans, C. briggsae , and D. rerio , and more gradual in Xenopus . In some cases, the slope even becomes negative at high temperatures, implying negative local activation energies. The best-fitting model varies by species: the double exponential (DE) function fits the more abrupt deviations better, while the simpler QE law fits gradual ones equally well ( Table 1 ). Panel E shows four later developmental intervals in D. melanogaster , including embryonic and pupal durations. These processes are more heterogeneous than early cleavages, yet they also deviate from Arrhenius scaling, especially over broader temperature ranges. Here too, the QE equation provides an excellent fit. View this table: View inline View popup Download powerpoint Table 1. Proposed modifications of the Arrhenius law for timings τ . Note that E 2 can be negative. Together, these findings show that non-Arrhenius scaling is a general feature of developmental timing, appearing in both periodic and non-periodic processes, across diverse taxa and regulatory network architectures. The variability in scaling patterns between species and developmental stages highlights the need for a general mechanistic framework that can account for such differences. To describe these deviations more accurately, several extensions to the Arrhenius model have been proposed. These include the QE equation ( 12 , 14 , 18 , 21 ), the DE equation (13) , the polynomial exponential (PE) equation (22) , and others ( 23 , 24 ). These models introduce additional temperature-dependent terms to improve data fits (see Supplementary Note 8). As shown in Figures 1C–E , both the QE and DE equations do a better job at describing cleavage and developmental timing than the Arrhenius equation. In species such as C. elegans, C. briggsae , and D. rerio , the DE law fits better due to abrupt deviations and negative local activation energies. In contrast, for X. laevis, X. tropicalis , and later stages of D. melanogaster , the QE law fits well with fewer parameters. These results suggest that different empirical models are suited to different biological contexts, depending on species, developmental stage, and underlying regulatory complexity. Many of the models described above are empirical, aiming to fit the experimental data without necessarily explaining the underlying mechanisms. The Arrhenius equation itself was originally proposed as an empirical description of the exponential increase in reaction rates with temperature. Only later was a theoretical basis developed by Kramers ( 25 , 26 ) and Eyring ( 27 ), who linked reaction rates to transitions over energy barriers in chemical reactions. More recently, we provided a mechanistic explanation for non-Arrhenius scaling in the early embryonic cell cycle of Xenopus laevis , using both embryos and egg extracts ( 12 ). By modeling the cell cycle as a relaxation oscillator, we showed that different activation energies among temperature-sensitive biochemical processes can account for the observed deviations from simple Arrhenius behavior. However, that framework is limited to oscillatory systems and does not extend to the temperature scaling of more complex, non-periodic biological processes. In complementary work, Voits and Schwarz ( 18 ) developed a statistical model that explains quadratic exponential scaling using a graph-based interpretation of mean first-passage times from a biochemical master equation. Such mechanistic approaches, drawing from dynamical systems and statistical physics, can help explain the physiological origins of temperature dependence in specific biological processes. In this study, we present a general mechanistic framework based on Markov chains to explain non-Arrhenius temperature scaling in biological timing. We first map the ordinary differential equation (ODE) model for the cell cycle oscillator ( 12 ) to an equivalent Markov chain model, recovering the period as a mean first-passage time. We then show that temperature-dependent multi-step cascades produce a mean first-passage time with Arrhenius flanks and a central QE regime, the latter consistent with prior graph-theoretic work ( 18 ). Applying this framework to 121 published datasets spanning diverse rates, durations, and species, we consistently observe the Arrhenius–QE–Arrhenius pattern, with dataset-specific boundaries. The QE regime is a robust consequence of averaging across many steps, whereas the flanking Arrhenius regimes are governed by smaller subnetworks and vary more across realizations. Simple motifs (cycles, competing paths) also account for apparent negative activation energies. Together, the ODE-to-Markov bridge and cascade theory provide a unified, testable account of how network structure shapes thermal responses in development and related processes. Results The early embryonic cell cycle oscillator’s period is controlled by reaction-specific temperature sensitivities A model based on coupled ordinary differential equations The early embryonic cell cycle of the frog Xenopus laevis (and other fish and amphibians) functions as an autonomous oscillator, with oscillations driven by the interplay between Cyclin B synthesis, its regulated degradation, and the activation of Cyclin-dependent kinase 1 (Cdk1). These oscillations are governed by nonlinear interactions within a regulatory network, consisting of Cyclin synthesis ( 29 ), a bistable switch for Cdk1 activation ( 30 , 31 ), and Cyclin degradation ( 32 , 33 ) ( Figure 2A ). The simplicity of this system allows for mathematical modeling of the cell cycle using ODEs parameterized from experimental data ( 28 , 34 ) (Supplementary note 1). These models successfully reproduce experimentally observed oscillations ( Figure 2C ), which form a closed trajectory, a limit cycle, in the (Cyc, Cdk1) phase plane ( Figure 2D , red/blue). The limit cycle orbits around the system’s nullclines (brown), illustrating the dynamics of the cell cycle. Download figure Open in new tab Fig. 2. The temperature scaling of represent rising (+) and falling (−) Hill functions, which depend on the variable x . (B) Markov chain approximation of one cycle of the 2-ODE embryonic cell cycle model. The MFP time E [𝒯], calculated between the initial state (00) and the final state (00) (passing through all intermediate states in an anti-clockwise direction), is identified with the period of the cycle. (Left: one-step approximation) (Right: n-step approximation.) (C) A numerical simulation of the 2-ODE cell-cycle model, demonstrating a limit cycle as a time series, for a specific parameter set. (D) A numerical simulation of the 2-ODE cell-cycle model, demonstrating a limit cycle of the relaxation-oscillator type in phase space, for a specific parameter set. The nullclines are also drawn. (E) Illustration of how the shifting of the nullclines of the 2-ODE cell-cycle model with temperature can create new fixed points that cause the model to get stuck and destroy the limit cycle. (Left: increasing T) (Right: decreasing T). (F) Comparison of the period of the Markov model and the ODE model. Different step sizes and the continuous limit with “infinitely” small step sizes are shown. The ODE model, uses the same specific parameter set as the other figures. (G) Transition times and of the Markov model in function of temperature. The temperatures where the model gets stuck in a high Cyc and Cdk1 state, or in a low Cyc and Cdk1 state, correspond to asymptotes (dashed black lines) of the continuous limit in (F). As the timescale of Cdk1 activation/inactivation is much faster than that of Cyclin synthesis and degradation, relaxation oscillations occur ( 35 – 41 ). During interphase, the system slowly ascends along the low-Cdk1-activity portion of the S-shaped nullcline, then abruptly jumps to the high-Cdk1-activity portion of the same nullcline ( Figure 2D , red → blue curve). Subsequently, the system descends along the high-Cdk1-activity branch due to active APC/C-mediated Cyclin degradation, before falling abruptly back to the low-Cdk1-activity branch to complete the cycle ( Figure 2D , blue →red curve). These dynamics produce characteristic sawtooth-shaped oscillations in Cyclin B levels ( Figure 2C , brown). This model has previously provided us with valuable information on the influence of temperature on the system ( 12 ). By introducing Arrhenius scaling to the rate constants ( k s , k d , k a , and k i ) (Supplementary note 1), the model allows to study how the oscillation period changes with temperature. When all activation energies ( E a ) are equal, the period exhibits Arrhenius-like scaling with a single E a . However, when activation energies differ across reactions, the model predicts deviations from Arrhenius scaling, particularly at high temperatures. These deviations, including negative effective activation energies, arise from the nonlinear interactions within the network and the differing temperature sensitivities of individual reactions. Two distinct mechanisms were identified as drivers of non-Arrhenius scaling and thermal limits in embryonic cell cycles. First, thermal limits can occur when one or more reactions deviate from Arrhenius scaling, exhibiting a thermal optimum. In such cases, the system’s behavior is predominantly shaped by the dynamics of these biphasic reactions ( 12 ). Second, even when all reactions follow Arrhenius-like scaling, differences in activation energies among reactions lead to emergent thermal limits and deviations from linear Arrhenius plots. For example, if the activation energy for Cyclin synthesis is higher than that for Cyclin degradation, the ratio k s /k d increases with increasing temperature. This causes the Cyclin nullcline to shift upward until it no longer intersects the middle portion of the S-shaped Cdk1 nullcline. At this point, oscillations cease, leaving the system in a M phase-like steady state with high Cdk1 activity ( Figure 2E , left). Conversely, if temperatures decrease, the Cyclin nullcline shifts downward, stabilizing an interphase-like steady state with low Cdk1 activity ( Figure 2E , right). Strikingly, experiments using cycling Xenopus extracts revealed that both mechanisms contribute to overall temperature scaling. Cyclin synthesis exhibited a biphasic temperature response, and an imbalance in activation energies was also observed: Cyclin synthesis was more temperature-sensitive than degradation ( 12 ). A model based on Markov jump processes The ODE analysis shows that temperature scaling can be a network property. However, extending this analysis beyond relaxation oscillators remains a significant challenge. Boolean models motivate a discrete state space for regulatory networks ( 42 ), but their discrete time hinders analysis of temperature-dependent periods. We therefore map the two-ODE oscillator to a continuous-time Markov jump process on a discrete state graph with protein concentrations having binary states: active (1) or inactive (0). We introduce a continuous-time framework into the state graph by assigning transition rates κ ij to transitions between states ( i ) → ( j ). For small time intervals dt , the probability of such a transition is given by: This describes a Markov jump process, a stochastic, continuous-time process on a discrete state space with minimal memory ( 43 ). The future state of the system depends only on its current state, making these processes both analytically and numerically tractable. In a Markov process, the time required to transition from state i to state j is called the first passage time 𝒯 ij . We are especially interested in the first passage time from the initial state 1 to the final state n , denoted as 𝒯 1 n = 𝒯. As 𝒯 varies between realizations, it is treated as a random variable. However, the mean first passage (MFP) time 𝔼 [𝒯] a fixed value for the process, can be calculated and used to study biological transition timings. For simplicity, we assume that the higher order moments of 𝒯 are negligible compared to the mean, a reasonable approximation for large-scale biological systems evolved for specific purposes. Temperature dependence is introduced into the Markov jump process by assigning Arrhenius scaling to the transition rates. Specifically, the inverse rate is modeled as: where A ij is the pre-exponential factor, E ij is the activation energy, R is the universal gas constant, and T is the absolute temperature in Kelvin. To approximate the two-ODE model, we developed a simplified Markov chain in which Cyclin B and active Cdk1 are limited to binary states: ON (1) or OFF (0) (Supplementary note 3). The ON state represents a high concentration of Cyclins and high Cdk1 activity, while the OFF state corresponds to a low concentration of Cyclins and low Cdk1 activity. This yields four possible system states: (0, 0), (1, 0), (1, 1), and (0, 1). Transition times τ ij were assigned based on the corresponding rate constants in the two-ODE model. For example, τ s represents the time for Cyclins to be synthesized [transitions (0, 0) → (1, 0) and (0, 1) → (1, 1)], while τ d governs Cyclin degradation [transitions (1, 0) → (0, 0) and (1, 1) → (0, 1)]. Similarly, ϵτ a and ϵτ i were used for the faster activation and inhibition of Cdk1, leading to transitions (1, 0) → (1, 1) and (0, 1) → (0, 0), respectively. To replicate the two-rate dynamics of Cyclin degradation, we in-troduced two pre-exponential factors, (low degradation when APC is OFF) and (high degradation when APC is ON). These factors were scaled to match the ratio of extremal values in the Hill function of the two-ODE model. The resulting Markov process is represented as a weighted directed graph ( Figure 2B , left). While backward reactions are possible, they were omitted in the transitions for Cdk1 activation and inhibition for simplicity. At a given temperature, the MFP time for the system to traverse all states in an counterclockwise direction, starting from (0, 0) and returning to (0, 0), was calculated and identified as the period of the cycle. The results were then compared to the temperature-dependent period of the two-ODE model. To improve the match between the Markov chain and two-ODE models, intermediate states were added for Cyclin B. Instead of binary states, Cyclin B concentrations were divided into n states, allowing transitions in steps of 1 /n . The right-hand side of Figure 2B shows the resulting state graph. As n increases, the MFP time of the Markov chain converges to the period of the ODE model ( Figure 2F ). In the continuous limit ( n →∞), the system exhibits vertical asymptotes at extreme temperatures (dashed lines), corresponding to thermal limits where oscillations cease. The Markov model provides an intuitive explanation for the breakdown of Arrhenius scaling when E s ≠ E d . Because transitions are stochastic in both directions, net forward drift is observed when the mean waiting times satisfy: When E d and E s differ, their scaling with temperature diverges, causing one of these inequalities to fail at extreme temperatures (see Figure 2G ). This results in the system becoming trapped in backward cycles during either Cyclin synthesis or degradation, halting oscillations (when n → ∞). Physically, this corresponds to one process dominating the other, preventing the system from completing a full cell cycle. In short, the Markov formulation reproduces the ODE period and the loss of oscillations at thermal extremes, while making a forward–backward bias explicit. This, in turn, motivates the three-regime scaling we derive next for linear cascades. Cascades of reversible transitions explain temperature scaling in developmental systems Motivated by the ODE to Markov mechanism above, we now examine linear cascades of reversible steps as a simplified model of developmental progression. This shifts the focus from periodic cycles to one-way progressions, where a system transitions from an initial to a final state through a series of potentially reversible steps. Many later developmental events can be conceptualized this way, as sequential, partly reversible processes rather than a dedicated oscillatory mechanism. To generalize this framework, we allow the transition times between states to differ. The resulting state graph, a linear cascade of cycles, is shown schematically in Figure 3A . It consists of n states, where each state i connects to the next state i + 1 with a forward transition time τ i,f and to the previous state i − 1 with a backward transition time τ i,b . 2 The process begins in state 1 and terminates in state n . For this linear cascade, the MFP time from 1 to n over all realizations of the Markov process can be calculated analytically: Download figure Open in new tab Fig. 3. The mean first passage time of a linear cascade of reversible processes exhibits three distinct temperature-scaling regimes. (A) Schematic of a linear cascade of reversible Markov jump processes. (B) Ingredients of the transitions times’ temperature dependence: (i) natural logarithms of forward (blue) and backward (red) transition times are distributed normally at the reference temperature T 0 , with the backward distribution shifted by Δ with respect to the forward one; (ii) activation energies of forward and backward transitions are drawn from normal distribution with the same mean and variance; (iii) resulting Arrhenius scaling of forward (blue) and backward (red) transition times (mean±s.d.). Gray vertival line indicates inversed reference temperature (1 /T 0 ). (C) Simulation of the temperature dependence of the MPF time for system (A): (i) Arrhenius scaling is fitted at high and low temperatures (Regions I and III, yellow lines). In Region II, the QE approximation from Eq. (7 ) provides a good fit. Boundaries of the QE region are determined from Eq. (8 ). (ii) MPF time plotted alongside with effective passage times (family of thin black curves) and the longest effective passage time (bold black curve). (D) Variation of MPF time temperature scaling across multiple random realizations of system (A): (i) MPF time from 100 realizations are plotted alongside with the QE approximation; (ii) Relative variation of MPF time (s.d. / mean) across realizations versus temperature for different sizes of a linear cascade (A). (E) Examples of individual transition time scaling when E i,f > E i,b (i) and E i,f < E i,b (ii). In both panels, blue and red lines show scaling of individual forward and backward transition line, respectively; blue and red shadings show the areas where either forward or backward transition from state i is more preferable, as shown by the respective pictograms. (F) Forward-to-backward transition time distribution illustrating the fraction of forward-dominated ( r 1, red bars) transitions, at different temperatures symmetric to the reference one. where When forward transitions dominate (i.e., r i = τ i,f /τ i,b ≪1), the MFP time simplifies to the sum of all forward transition times, 𝔼 [𝒯] ≈ τ i,f . However, when backward transitions become significant τ i,b = 𝒪 ( τ i,f ) or τ i,b < τ i,f , 𝒯 increases, as each effective passage time is amplified by a factor that accounts for the ratios of forward and backward transitions in prior cycles. This factor is always greater than 1, meaning faster backward transitions increase the MFP time considerably. To model large-scale biological processes, we fix a reference temperature T 0 and generate transition statistics in three steps (see also Supplementary Note 5). Baseline times at T 0 : for each step i , we draw for-ward transition rates ln τ i,f ( T 0 ) from and back-ward transition rates ln τ i,b ( T 0 ) from , where sets a forward bias at T 0 (forward faster than backward), see Fig. 3B(i) . Activation energies: independently for each directed transition, we draw activation energies E i,f and E i,b from , with biologically plausible values (e.g., µ E ≃ 50 kJ mol −1 , σ E ≃ 10 kJ mol −1 ; Ref. ( 44 ); see Fig. 3B(ii) ). Arrhenius extrapolation from T 0 : for any temperature T , the transition times are given by which yields the forward/backward (= f/b ) time distributions shown in Fig. 3B(iii) . In the limit of many steps ( n → ∞) this leads to a QE scaling law ( Eq. (70 ) of the Supplementary note 5): where A is a prefactor determined by the baseline statistics at T 0 (and depends on n, µ E , and σ E ). Although illustrated here with normally distributed activation energies, the QE form is generic, arising for broad classes of energy distributions beyond normality. The QE scaling law is valid within a specific temperature range: where forward transitions dominate, meaning ⟨ r ⟩ = ⟨ τ i,f /τ i,b ⟩ < 1. Here, T c and T h are the temperature limits of the QE scaling at cold and warm temperatures: Note that the above expressions for the temperature limits T h and T c only have a physical meaning for , i.e., when the distributions of forward and backward transition times at the reference temperature are well-separated. Within this range, the MFP time follows a QE scaling law as described in Eq. (7 ) ( Figure 3C ). This behavior arises from the collective contribution of many forward-dominated transitions with comparable timescales in the sum given by Eq. (5 ). As a result, the system displays strongly reduced variability across realizations, particularly as the number of steps increases ( Figure 3D , Figure S3 ). This regime is therefore both more predictable and more robust, reflecting the emergent regularity of large ensembles of similarly activated steps. The validity of the QE regime largely depends on the relative width of the activation energy distribution and mean directional bias |Δ|. In contrast, outside this range, the distribution of effective transition times includes large tails, outliers dominate the total MFP time, and the law of large numbers is no longer valid ( Figure 3C ). In these extreme temperature regions, the MFP time scales in an Arrhenius-like manner with either a positive or negative constant activation energy, depending on the temperature limit: positive activation energy in the cold regime and a negative activation energy in the warm regime. In both cases, the dynamics are effectively dominated by the longest rate-limiting effective passage time (panel (ii) in Figure 3C ). This leads to strong sensitivity to stochas-tic variation: despite identical distributions of activation energies across realizations, substantial variability emerges in the observed slopes, transition sharpness, and crossover temperatures. This variability persists even for large systems (large n ), indicating that the statistical averaging in long cascades does not apply in these regimes ( Figure 3D , Figure S3 ). To better understand where the deviations from the QE scaling law originate, let us zoom into the scaling of individual reactions. Since the activation energies of forward and backward transitions of individual reaction are generally unequal, two cases are possible: (1) E i,f > E i,b and (2) E i,f < E i,b . In both cases, forward transitions are preferable around the reference temperature T 0 , favoring the QE approximation. However, backward transitions in case (1) dominate at lower temperatures, while the opposite holds for the case (2) ( Figure 3C , Figure S2 , and Supplementary note 5). Thus, moving away from T 0 in either direction facilitates the backward-dominated transitions with r i = τ i,f /τ i,b > 1 through either case-(1) reactions on the colder side or case-(2) reactions on the warmer side ( Figure 3F ). Such backward-dominated transitions trap the cascade in the local cycle motifs, thereby extending overall MFP time beyond the QE prediction at the temperature extremes. Together, these results show that complex biological processes can transition between distinct scaling regimes ( Figure 4A ). Quadratic exponential scaling arises in intermediate thermal regimes due to cumulative contribution of many steps in large cascades, while single-exponential (Arrhenius) behavior dominates at thermal extremes due to the disproportionate influence of a small number of transitions. This framework provides a mechanistic basis for understanding both the robustness and variability observed in temperature-dependent biological timing. Download figure Open in new tab Fig. 4. General scaling motifs underlying temperature-dependent biological timing. (A) Schematic overview showing how biological systems combine quadratic exponential scaling at intermediate temperatures — driven by the collective averaging of many similar transitions — with Arrhenius-like single-exponential scaling at low and high temperature extremes, where one or few rate-limiting steps dominate. (B) Sequential two-step motif: at low temperatures, the slower high-barrier reaction sets the overall rate; at high temperatures, the faster low-barrier reaction dominates, but both regimes maintain positive effective activation energies. (C) Cycle motif: at high temperatures, reverse transitions become dominant, trapping the system and producing an emergent negative macroscopic activation energy despite all microscopic steps following Arrhenius scaling. Here, E 23 > E 12 and E 12 + E 23 < E 21 . (D) Feedforward motif: two parallel forward pathways contribute, with temperature shifting dominance between them, producing a negative activation energy within a narrow intermediate temperature range. Simple network motifs explain the mechanisms behind non-Arrhenius scaling and negative activation energies Figure 4A illustrates how single-exponential Arrhenius behavior, whether with positive or negative activation energies, dominates at thermal extremes, due to the disproportionate influence of a small number of rate-limiting transitions. We focus here on three minimal network motifs that play a critical role in generating deviations from simple Arrhenius responses at these temperature extremes. The simplest case is two sequential reactions with different activation energies ( Figure 4B ). At low temperatures, the system’s overall rate scales with the reaction having the largest activation energy, while at high temperatures, it scales with the smallest. This produces the widely used double-exponential response curve. However, at both extremes, the system still exhibits scaling with a positive activation energy. So how can minimal network motifs, where each microscopic transition obeys the Arrhenius law with a fixed positive activation energy, give rise to a system where the MFP time scales with a negative activation energy? We identify two minimal state-graph motifs that illustrate how negative activation energies can emerge from systems with only positive microscopic activation energies (see Supplementary Note 7). The first is a simple cycle ( Figure 4C ), where the system can transiently loop between states 1 and 2 before reaching the endpoint 3. At high temperatures, when reverse transitions become significantly faster than forward ones, the MFP time scales with an effective negative activa-tion energy. The second motif involves two parallel forward-directed paths ( Figure 4D ), where temperature determines the likelihood of taking one path over the other. Within a narrow intermediate temperature range, the shift in path dominance produces a negative activation energy. These motifs demonstrate how macroscopic negative activation energies can arise from underlying microscopic processes that individually follow Arrhenius scaling, emphasizing the importance of network topology and transition probabilities in shaping overall temperature dependence. Crucially, both motifs feature multiple possible routes to the endpoint, with the most probable path shifting with temperature. In contrast, systems with only a single possible path, such as the two-step sequential motif, cannot generate negative macroscopic activation energies in their MFP scaling. The temperature-dependent scaling in the cycle and feed-forward motifs is qualitatively distinct. For the cycle, the macroscopic activation energy can become negative in the high-temperature limit, as the system becomes increasingly trapped in reverse loops. For the feedforward motif, negative activation energies emerge only in an intermediate transition region, where the dominant contribution to the MFP time shifts between two competing forward paths that control behavior at opposite thermal extremes. The Markov framework captures temperature scaling across diverse species, biological timescales, and temperature ranges To evaluate the applicability of the Markov cascade model, we assembled a diverse set of experimental datasets densely sampled in temperature, i.e. many measurements per curve across a broad thermal span, covering growth and developmental processes, including early embryonic cell cycles and stage transitions across multiple species ( 12 – 15 , 45 – 88 ). In total, we analyzed 121 datasets from 48 sources, comprising 2174 measurements with ≈18 temperature points per curve (and minimum 5), with a typical temperature span of 30 ± 17 K (minimum 10K). Each temperature response curve was fit with a composite model: a central QE regime [Eq.( 7 )], flanked by single-exponential regimes at low and high temperatures (See Supplementary note 8). This approach enables a unified description of thermal scaling across biologically relevant ranges, with transitions between regimes governed by the underlying distribution of activation energies and the balance of forward and backward transitions in the Markov chain. Figure 5A shows representative examples, highlighting the diversity of observed scaling behaviors, from processes following QE across the full range to cases with deviations at either thermal extreme. Across all datasets, the model provided ex-cellent fits (median R 2 = 0.98, 90% CI [0.89, 0.99]), accurately capturing both gradual and sharp transitions in scal-ing. The majority of experimental temperature responses, approximately 59%, were fully described by the QE law, while the rest displayed deviations at temperature extremes: 32% only at warm temperatures, 6% only at cold temperatures, and 3% at both ends ( Figure 5B ). This pattern suggests that while many biological systems operate within a robust QE regime near their thermal midpoint, transitions to Arrhenius-like scaling are common at extreme temperatures. At the same time, datasets exhibiting deviations at colder temperatures were less represented among the analyzed sources. It may reflect the fact that the processes slow down faster at the colder side, making their experimental measurement more challenging. Download figure Open in new tab Fig. 5. The Markov cascade model explains temperature scaling across diverse biological systems. (A) Example fits of experimental data using a composite model with a central quadratic exponential (QE) regime flanked by exponential regimes at low and high temperatures. (B) Proportions of N = 121 datasets classified as fully QE or showing deviations at warm, cold, or both extremes (left chart) and as positive or negative scaling (right chart). (C) Distributions of fitted QE parameters, including mean activation energy µ E , variance and E c . (D) Scatterplot of datasets in space, with markers indicating taxonomic groups (see Legend) and colors indicating scaling regime (QE or deviations, similar as in B). (E) Scatterplot similar as in (D) but in ( µ E , T 0 ) space. Kernel density estimation compares reference temperature T 0 in datasets with positive and negative scaling ( p < .01 via the Mann-Whitney test; solid lines indicate sample medians). (F) Temperature response curves normalized to their value at reference temperature T 0 , overlaid and colored by µ E , showing conserved yet distinct scaling behaviors across datasets. Figure 5C presents the distribution of key fitted parameters. The mean activation energy of the QE regime, µ E ( Eq. 7 ), typically ranged from 40–80 kJ/mol, consistent with systems composed of many elementary steps, each drawing from a moderate activation energy distribution (0–200 kJ/mol). These observed mean activation energy values also overlap those reported by Dell et al . ( 10 ). The width of the distri-bution, , varied across systems and determined the size of the QE window as well as the scaling curvature. For temper-atures outside this window, we extracted effective activation energies from the flanking exponential regimes, confirming that distinct scaling laws govern different parts of the thermal response. We further visualized these relationships in the 2D parame-ter space , shown in Figure 5D . Here, we marked different major clades (Archaea, Bacteria, Fungi, and Plan-tae) with different symbols, and indicate data from frogs (including frog egg extracts), Insecta, Nematoda, and all other sampled groups separately. Most datasets cluster around a positive mean activation energy µ E with moderate variance . Two anomalous groups emerge: (i) systems with exceptionally large variance and (ii) systems with a negative mean activation energy of the QE regime ( µ E < 0). These outliers are enriched for phytoplankton, algae, bacteria, viruses, thermophilic fungi, and heat-tolerant plants, taxa that routinely experience severe thermal conditions. This pattern suggests distinct temperature-compensation strategies: a large is consistent with stronger compensation across the operative thermal window, whereas the biological interpretation of µ E < 0 is less direct. The plot of the reference temperature T 0 against µ E ( Figure 5E ) shows that the µ E 0, consistent with a change of the operating range towards warmer conditions. Consistently, responses normalized by T 0 and colorcoded by µ E ( Figure 5F ) reveal that µ E sets the local slope at T 0 , effectively sliding systems along the QE curve: positive µ E aligns with the cooler-temperature side, whereas µ E < 0 aligns with the warmer-temperature side (see schematic in Figure 5A ). Negative quadratic exponential temperature scaling can be recapitulated in the Markov cascades Data reveal two distinct settings in which the apparent activation energy can be negative: (i) at the hot extreme , where reverse loops in cyclic motifs dominate and the macroscopic slope becomes negative; and (ii) within the intermediate QE range , where a local negative slope is observed over a finite window. The hot-edge case has a simple mechanistic account via reverse-loop trapping in cycles (cf. Fig. 4C ). In contrast, the mechanism in the QE window case is unclear. We therefore explore thermodynamically consistent network options, where all microscopic steps obey Arrhenius with positive activation energies, showing how negative apparent E a can appear within the QE window. The most direct way to reproduce negative scaling is to implement a linear Markov cascade in which the activation energy distribution is mirrored around zero ( Figure 6A–C ). Drawing microscopic one-step transition energies from this distribution generates negative linear scaling of individual reactions, producing a negatively scaled temperature response curve that qualitatively resembles those observed in the experimental datasets ( Figure 5A, F ). However, interpreting a negative mean activation energy simply as the average over negative microscopic transition energies can be misleading, since such a perspective may contradict fundamental thermodynamic principles. Download figure Open in new tab Fig. 6. Negative quadratic-exponential temperature scaling can be recapitulated in the Markov cascade of cycles. (A)-(C) Negative MFP time temperature scaling results from the negative scaling of individual reactions in a linear Markov cascade. (A) Schematic of a linear cascade, where each reaction has negative activation energy ( E { f,b } < 0). (B) Temperature scaling of individual reactions: (i) normal distribution from which activation energies E { f,b } are drawn; (ii) “negative” scaling of individual forward and backward reactions in blue and red, respectively (mean±s.d.). (C) Resulting temperature scaling of MFP time in the linear cascade (blue). The bold red curve shows QE approximation, the thin black curves display effective passage times, and the shading indicates the range where MFP time scaling follows the QE approximation. (D)-(F) Cascade of cycles with an apparent negative scaling branch composed of positively scaled reactions can also exhibit negative MFP time scaling. (D) Schematic of a linear cascade of cycles, where each reaction has positive activation energy ( E { f,b } > 0), yielding negative apparent activation energy of cycles ( E cycle < 0). (E) Temperature scaling of individual reactions: (i) normal distributions from which activation energies E { f,b } are drawn; (ii) “positive” scaling of individual forward and backward reactions in blue and red, respectively (mean±s.d.). (F) Resulting temperature scaling of MFP time in the linear cascade of cycles (blue). The bold red curve shows QE approximation, the thin black curves display individual cycle transition times, and the shading indicates the range where MFP time scaling follows the QE approximation. We propose instead that systems with negative µ E may reflect a higher degree of regulation, in which transitions between states are multi-step and can thus give rise to negative apparent activation energies. To illustrate this, we turn to minimal network motifs composed solely of positive-energy transitions, for which we have previously demonstrated the possibility of negative scaling ( Figure 4B–D ). Without loss of generality, we focus on a cycle motif and consider a Markov cascade of such cycles ( Figure 6D ). In this setup, forward and backward one-step transitions, each positively scaled but dif-fering in mean activation energies ( Figure 6E ), together produce negative scaling of individual cycles at higher temperatures and positive scaling at lower temperatures ( Figure 6F ). As expected, the MFP time of a cascade of such cycles then displays negative QE scaling within a defined temperature range under physically realistic conditions. The simple cycle examined here serves merely as one example, as many other multi-step motifs could plausibly underlie negatively scaled reactions. Conclusion Temperature fundamentally shapes biological timing. While the Arrhenius equation captures how individual activated steps accelerate with temperature, its limits become clear across broader thermal ranges and complex multi-step processes. Here, we introduced a Markov cascade framework that provides a mechanistic account of these departures. Across linear cascades of reversible steps, the MFP time exhibits a consistent three-zone pattern: a QE regime at intermediate temperatures, with Arrhenius-like behavior at cold and hot extremes. The QE regime arises from statistical averaging over many temperature-sensitive steps and is robust to distributional details, yielding predictable and low-variability timing. At thermal extremes, averaging breaks down and a few effective steps dominate, restoring single-exponential (Arrhenius) scaling and increasing variability. We derived analytical boundaries between these regimes in terms of activation energy statistics and forward-backward bias, and we used an ODE → Markov mapping of the embryonic cell-cycle oscillator to reproduce the observed thermal limits of oscillations. Apparent negative activation energies need not imply negative microscopic barriers or enzyme denaturation. In our framework they emerge naturally from network topology. Simple motifs, such as cycles or feedforward parallel pathways, are sufficient to generate such effects, emphasizing how global thermal responses can emerge from local, reversible steps. The theory agrees with data: applied to more than 100 densely sampled datasets spanning taxa, rates, and developmental intervals, the composite QE–Arrhenius picture provides excellent fits, explains when coordination is preserved, and clarifies why breakdowns can occur near physiological limits. These results offer practical guidance for measurement (need for dense temperature sampling) and interpretation (treat warm-edge slowdowns and negative slopes as network effects unless denaturation is independently demonstrated), and they suggest how thermal robustness may be organized by network architecture in natural systems. Our QE prediction is consistent with graph-theoretic first-passage analyses of large networks ( 18 ). We specify when QE holds for finite linear reversible cascades, explain the reemergence of Arrhenius scaling at the extremes and motif-driven negative apparent activation energies, and validate these predictions across a broad empirical dataset. More broadly, the Markov cascade model provides a unified, testable framework for forecasting how biological timing shifts with temperature. By linking microscopic activation statistics to macroscopic timescales, it identifies thresholds of developmental failure, guides targeted perturbations (e.g., altering network topology or biases), and can be extended beyond development to other multi-step physiological and ecological processes in a warming world. Discussion Robustness of biological timing to temperature changes Development Embryonic development is a cascade of coordinated biochemical events, from molecular regulation to large-scale morphogenesis. Despite large changes in stage duration with temperature, the final outcome often remains remarkably consistent. A striking example is Drosophila , where the relative durations of embryonic and pupal stages are preserved across a broad thermal range ( 14 , 70 ) ( Figure 7A ). Similar invariance is seen in X. laevis embryos ( 12 , 14 ) ( Figure 7B ) and in Medaka fish somitogenesis, where somite spacing remains constant even though the period varies more than threefold ( 89 ). Related results show that transcriptional output and cell cycle timing scale proportionally in early Drosophila embryos ( 90 ), preserving gene expression levels across temperatures. Download figure Open in new tab Fig. 7. Robustness of biological timing to temperature change. (A) Mean durations of five stages of early embryogenesis in X. laevis (Data from Ref. ( 12 )). (B) Mean durations of four developmental stages in D. melanogaster (Similar to Figure 5E , data from Refs. ( 14 , 15 )). (C) Simulated mean durations of successive developmental stages using linear Markov cascade of increasing size (from 50 to 4000 steps) overlayed with individual realizations. In (A)–(C), left panels display scaling of biological timing per stage, while panels to the right show centered duration, i.e., duration with extracted mean over temperature range per stage, overlapped with the universal scaling curve. While these examples highlight striking robustness, recent large-scale single-cell profiling in zebrafish shows that this resilience is not absolute: certain cell types are disproportionately sensitive to thermal stress, leading to structural defects ( 91 ). Such findings illustrate that robustness at the whole embryo scale can coexist with local vulnerabilities, and that temperature sensitivity can propagate from molecular processes to developmental coordination ( 92 ). The preservation of relative timing despite variable absolute rates is often attributed to active temperature compensation, such as opposing enzymatic reactions with tuned activation energies ( 93 ). Our results suggest that it can also arise passively from network structure. In the QE regime, timing emerges from averaging over many similarly scaled, reversible steps, which suppresses variability and preserves stage-to-stage timing ratios even when overall rates change. Simulations ( Figure 7C ) show that this holds even if stages differ in step number and reaction rates, and remains robust to stochastic variability in individual rates. At thermal extremes, this robustness breaks down as averaging gives way to a few rate-limiting steps. Variability increases, and relative stage durations can shift, a pattern observed experimentally in some systems at high or low temperatures ( 93 ) and consistent with the idea that coordination is maintained only within the Arrhenius-like regime of physiological rates ( 13 ). In our framework, such deviations reflect a transition from QE behavior to Arrhenius-like regimes at the edges of the viable temperature range. Temperature compensation In chronobiology, the maintenance of relative timing across temperature is sometimes called temperature compensation of phase relationships. For example, the neuronal oscillator in the crab Cancer borealis preserves phase relationships across temperatures even though its period changes ( 94 , 95 ). A stricter definition is used for systems such as the circadian clock, where the pe-riod itself is temperature-independent ( Q 10 ≈1) ( 96 – 99 ). Mechanistic explanations for strict compensation include an-tagonistic balance—opposing processes with matched activation energies ( 100 , 101 ), and enzyme-limited schemes, where several reactions share a common catalytic bottleneck ( 102 , 103 ). Network-based models have also shown that temperature robustness can emerge from architecture and feedback alone ( 104 – 106 ). Our framework extends these ideas: compensation can emerge not only from fine-tuned opposing processes but also from distributed averaging across many steps with constrained activation energy distributions. In some cases, this averaging naturally produces similar scaling for opposing reactions, leading to antagonistic balance without explicit tuning. In others, large control networks in the QE regime can exhibit very low apparent activation energies, , balancing mean and variance in individual reaction activation energies to achieve robustness. Negative activation energies Negative apparent activation energies have long been reported in experimental studies of thermal scaling, particularly at warm temperatures, but their mechanistic basis has remained unclear. At temperature extremes, they may arise from processes that deviate from simple Arrhenius kinetics, for example, biphasic behavior such as protein denaturation ( 107 , 108 ), or nonlinear system dynamics approaching bifurcations ( 12 ). In Markov networks, cycles and feedforward motifs ( Figure 4C–D ) can likewise produce negative slopes when transitions become dominated by reverse steps or competing paths. Within the QE regime, negative activation energies can emerge more subtly as apparent network-level properties. Although each microscopic step obeys Arrhenius scaling with a positive energy, motifs such as asymmetric cycles generate local negative slopes, and cascades of such motifs can display an overall negative µ E ( Figure 6 ). These values therefore do not reflect physically negative microscopic energies, but rather emergent signatures of distributed, multi-step regulation. Biologically, such mechanisms may explain why negative µ E values are enriched in organisms adapted to extreme or fluctuating environments ( Figure 5D ), where they may represent adaptive strategies for maintaining coordination under thermal stress. Evolutionary and ecological implications Across approximately 100 datasets, we find that many biological systems predominantly operate within the QE regime, yet deviations, particularly at elevated temperatures, are frequently observed. This pattern suggests that biological systems may have evolved to function optimally within a buffered thermal window, outside of which developmental timing becomes increasingly erratic or prone to failure. Supporting this idea, recent work ( 12 ) showed that, across multiple species, the native environmental temperature closely aligns with the minimum of the QE scaling curve. In other words, these organisms appear to have evolved to maximize developmental speed while remaining within the more robust, low-variability region of the QE regime. These findings carry important implications for understanding how organisms adapt to thermal environments and for predicting the biological consequences of climate change. Consistent with this view, Begasse et al. ( 13 ) observed that the fertile temperature window in C. elegans and C. briggsae overlaps with the range in which cell division adheres to Arrhenius kinetics, highlighting a correlation between developmental robustness and fertility limits. Our results emphasize the importance of characterizing thermal response curves at high resolution, not just at a few “representative” temperatures, to uncover the underlying structure of thermal sensitivity and to identify critical tipping points in developmental systems. Connections to empirical biophysical models Our results also relate to the Sharpe–Schoolfield model, a biophysically inspired framework that combines temperature-dependent reaction kinetics with reversible enzyme inactivation at temperature extremes ( 23 , 24 , 109 – 111 ). Like our model, it captures three distinct thermal regimes: intermediate speeding up with quadratic-like curvature, a linear scaling at cold extremes, and a slowdown with negative activation energy at hot extremes. However, it attributes the deviations at both ends to enzyme inactivation. In contrast, our Markov framework shows that similar triphasic scaling can emerge from statistical averaging and network effects, even without assuming enzyme denaturation. This provides a mechanistic foundation for empirical patterns observed across taxa. Other commonly used models, such as the QE and DE, typically capture only one or two scaling regimes; while useful for fitting smooth or abrupt deviations, they often fail to describe the full thermal range seen in biological data. Toward predictive models of temperature scaling Together, these insights highlight the need for mechanistic models that explain not only what scaling behavior occurs, but why, and how it might shift under environmental change. By combining statistical physics with biological realism, the Markov cascade model provides a versatile tool for studying temperature-dependent timing. It unifies empirical observations across diverse taxa and stages, and offers mechanistic explanations for scaling deviations that previously lacked theoretical grounding. As global temperatures continue to rise, such generalizable models will be essential for identifying thresholds of developmental failure and for predicting species-specific vulnerabilities. Data availability statement The numerical codes to reproduce the figures in this study are openly available in G it L ab ( 112 ), and as an archived repository on RDR by KU Leuven [Upcoming]. The experimental datasets used and their analysis is also openly available on RDR [Upcoming]. Supplementary Note 1: The two-ODE embryonic cell cycle oscillator model The two-ODE cell cycle model is given by the following set of ordinary differential equations Here, and are rising and falling Hill functions, re-spectively, representing an increase or decrease in the activity of protein i in response to changes in x : The parameters a i , b i , n i , and K i denote the basal activity, maximum increase in activity, Hill coefficient, and activation threshold, respectively ( 12 ). A. Nullclines and dynamics To analyze the system’s behavior, we examine its nullclines—curves in the ( cyc, cdk1 ) phase space along which the time derivative of one variable is zero. The cyc -nullcline is defined by setting ,producing a curve in which the synthesis and degradation of cyc are balanced. Similarly, the cdk1 -nullcline is obtained by setting , describing where activation and inactivation of cdk1 are in equilibrium. The intersection points of the nullclines correspond to steady states of the system. Depending on the location and stability of these points, the system may exhibit monostability, bistability, or oscillatory behavior. In particular, the interaction of fast positive feedback and slower negative feedback can give rise to a stable limit cycle—an isolated, closed trajectory in phase space that corresponds to sustained oscillations in time. As shown in Figure 2C–D , the concentrations of cyc and active cdk1 oscillate periodically over time (panel C), and the system traces out a loop in the ( cyc, cdk1 ) phase plane (panel D). This loop follows a characteristic relaxation oscillator pattern: slow accumulation phases alternate with rapid transitions, producing sharp transitions between activity states (e.g., from low to high cdk1 activation). B. Temperature dependence The constants k i are rate constants which are assumed to have an Arrhenius dependence: This introduces a temperature dependence in the model and allows to study the temperature dependence of the limit cycle oscillations. An example in which the period has been plotted in function of the temperature is shown Figure 2F (ODE model). In numerical simulations, it has been observed that when a limit cycle is observed at intermediate temperatures, when E d ≠ E s it gets destroyed in the high- or low-temperature limit( 12 ). This can be explained, based on the relative positions of the nullclines in ( Cyc,Cdk1 ) phase space. The position of Cyc -nullcline is determined by the ratio If E d ≠ E s , then the Cyc -nullcline shifts up or down with changing temperature untill it interesects the cdk1 -nullcline in the upper or lower branc, which creates a fixed point and destroys the limit cycle. This is shown in Figure 2E for the case E d > E s . C. Simulations Figure 2 in the main text contains a number of simulations of this model and its temperature dependence (C: time series, D: phase space, E: nullclines, F: temperature dependence). For all of these simulations we used the parameters in Table S1 which are taken from ( 12 ). Supplementary Note 2: Markov jump processes and mean first passage times In this work, we approximate biological systems as Markov jump processes. Generally speaking, a Markov process is a stochastic process characterized by the key property of memorylessness: the future state of the system depends only on its present state, not on its past history. That is, for a random time-dependent variable X ( t ) undergoing a Markov process, the probability distribution of X ( t + δt ) depends only on the current state at time t . In Markov jump processes specifically, X ( t ) can only occupy one of a finite set of discrete states S = {1, 2, …, n }, while transitions between states occur at random times in continuous time t . The likelihood of a transition from state i to state j is governed by the transition rate κ ij . Informally, if the system is in state i , the probability of transitioning to state j within a small time interval δt is approximately P [ i → j ] ≈ κ ij δt . Since we are particularly interested in timing, it is useful to define a transition time: This transition time τ ij represents the average time to make the transition from state i to state j , assuming no other transitions from i occur during that period. A Markov process can be conveniently represented as a weighted, directed graph: nodes correspond to discrete states, and directed edges with weights κ ij represent transitions i → j 3 . Examples of such graphs are shown in Figures 2 – S1 . An equivalent representation is the generator matrix G , with entries G ij = κ ij for i ≠ j and G ii = ∑ j≠ i κ ij . We aim to study the temperature dependence of biological processes. To this end, we postulate that each transition time in the Markov jump process is temperature-activated according to an Arrhenius law: These transitions should be understood as small-scale thermodynamic events that are sufficiently well described by the Arrhenius equation within the relevant temperature range. Real biological systems are not accurately described by arbitrary state graphs. In regulatory developmental systems in particular, we know they have evolved to reliably reach a specific final state under appropriate conditions ( 113 ). We therefore restrict our analysis to simply connected, directed state graphs in which state n is the only possible terminal state. In the language of Markov chains, n is the unique absorbing state, with total escape rate . We also assume that the initial state is known and fixed, so that the system always starts in state 1: X ( t = 0) = 1. In this setting, a natural quantity to study is the time 𝒯 it takes for the process to reach the final state X = n for the first time, starting from X = 1. This is known as the first passage time. Since 𝒯 depends on the stochastic realization of the process, it is a random variable. Markov theory provides tools to compute its statistical moments ( 18 , 43 ). The most informative of these is the mean first passage (MFP) time, 𝔼 [𝒯], particularly because directed biological systems are expected to minimize higher-order moments through evolution. The MFP time from state i to the absorbing state n can be computed by solving the following system of equations: where T i = 𝔼 [𝒯 | X (0) = i ] is the MFP time from state i to n , and G is the generator matrix of the Markov process. Supplementary Note 3: Markov approximation of the two-ODE embryonic cell cycle oscillator The goal of this supplementary note is to give a detailed overview of how certain features of the two-ODE embryonic cell cycle oscillator can be captured in an approximating Markov jump process. The approximation becomes better as the number of discrete states in the process is increased. A. Construction of a four-state Markov chain To accurately capture the behavior of a temperature-dependent biological system, we need to describe its state and how it evolves over time and with temperature. In the ODE formulation of the embryonic cell cycle model, the state is specified by the continuous concentrations of the proteins Cyc and Cdk1. The system of differential equations in Eq. (12 ) governs their time evolution, and temperature dependence enters through the rate constants in Equation 14 . We aim to approximate the behavior of this ODE system using a Markov jump process. The first step is to replace the continuous state space with a discrete one. In first instance, we limit Cyc and Cdk1 to two binary states: ON (1) or OFF (0). The ON state represents high Cyclin concentration and high Cdk1 activity, while the OFF state corresponds to low Cyclin concentration and low Cdk1 activity. Combining the states of the two proteins yields a discrete state space with four possible configurations: (0,0), (1,0), (1,1), and (0,1) (see Figure 2B , left). The next step is to introduce time evolution that mimics the behavior of the ODEs. If the ODE model is in a low- Cyc state (0,X), it can transition to a high- Cyc state (1,X) via the synthesis rate constant k s . Conversely, it can return from (1,X) to (0,X) due to degradation governed by k d . We therefore introduce two corresponding transition times, τ s and τ d , to allow switching between low and high Cyc states in the discrete-state Markov chain. We note that amount of degradation in the ODE model is regulated by cdk1 via the Hill function .To reflect this, we distinguish between two regimes: for low cdk1 ((1,0) → (0,0)), and for high cdk1 ((1,1) → (0,1)). The fast switch from low to high cdk1 activity ((X,0) → (X,1)) in the ODE model is controlled by the activation term with rate constant k a , which becomes active only when Cyc concentration is high. For this reason, we introduce a fast transition time ϵτ a between (1,0) and (1,1) in the Markov chain. Finally, to complete the model, we add a fast transition time ϵτ i from (0,1) to (0,0), mimicking cdk1 inhibition controlled by the rate constant k i (see again Figure 2B , left). The final step in constructing the Markov jump process is to introduce temperature dependence. We want the transition times to reflect the same temperature dependence as the rate constants in the ODEs. This can be achieved by assuming an Arrhenius form where E i is taken to be the same as for the corresponding rate constant k i . The prefactors A i are undetermined at this stage. The only physically motivated constraint is on the ratio between the fast and slow degradation regimes at the reference temperature. This ratio is chosen to match the ratio of minimum to maximum activity of the Hill function There is no direct relationship between the remaining preexponential factors and the rate constants . Their numerical values can be determined by fitting to the period of the ODE model (see below). Once the model is constructed, we extract a period that can be compared to that of the original ODE model. This is done by computing the mean time it takes the Markov chain to return from state (0,0) to (0,0), passing through all other states in the process. Mathematically, this is equivalent to computing the MFP time to (0,0)* starting from (0,0) in a linear cascade of states such as in Figure S1A . From Eqs. (5 ) and ( 6 ), we obtain When we fill in all parameters and compare this temperature-dependent MFP time with the period of the ODE model, we find that the ODE model exhibits a stronger temperature dependence (see Figure 2F , one step vs. ODE model ). Specifically, in the Arrhenius plot, the ODE model shows a steeper negative slope in the high-temperature limit and a steeper positive slope in the low-temperature limit. This difference arises because four discrete states are insufficient to fully capture the complex temperature dependence of the continuous ODE model. B. A linear cascade of multiple reversible reactions generates a strong temperature dependence We hypothesize that adding more intermediate states could yield a better approximation of the continuous ODE model. Therefore, in this section, we study in more detail a linear cascade of n reversible transitions, with forward transition times τ f , backward transition times τ b , and a final transition time τ 0 , as illustrated in Figure S1D . We define the ratio If E f < E b , this ratio increases in the high-temperature limit; if E b < E f , it increases in the low-temperature limit. The MFP time from the initial to the final state is given by Evaluating the sums yields the closed-form expression: In the limit where r ≪ 1, backward transitions can be neglected, and the MFP time simplifies to the sum of all forward transition times and the final transition time: If r ~ 1, the MFP time begins to increase. Physically, this corresponds to the system becoming trapped in cycles of forward and backward transitions: when r ~ 1, the system frequently revisits previous states, and as the number of intermediate states grows, the number of potential loops increases, leading to a longer time to reach the end state. When r ≫ 1, the MFP time is approximately and increases rapidly compared to the r ≪ 1 case. Note that in this regime, the dependence on n becomes exponential. Since r is defined by Eq. (21 ), it has an exponential dependence on temperature. This implies that the transition from the r ≪ 1 to the r ≫ 1 regime occurs abruptly, producing a sharp, temperature-dependent increase in the MFP time in the thermal limits where r becomes large. Moreover, this effect is amplified for larger values of n . For the regime r < 1, we can also consider the large- n limit of Eq. (23 ). To obtain a finite expression in this limit, we rescale the forward and backward transition times as and . Physically, this reflects the idea that introducing an infinite number of intermediate states only makes sense if the transition times become infinitesimally small. Note that . Taking the limit n → ∞ while holding and fixed yields: For future reference, we note that this large- n limit does not yield a finite expression when r ≃ 1. C. Adding extra steps to the Markov chain reproduces the temperature dependence of the ODE model In the previous section, we showed that adding intermediate states to a linear cascade of reversible processes can produce a steeper temperature dependence of the MFP time in the thermal limits. This motivates us to extend the four-state Markov chain by adding additional steps, such that the MFP time of the chain more closely matches the period of the ODE-based cell cycle oscillator when plotted in an Arrhenius diagram. In numerical simulations, we observe that Cdk1 activation and inhibition occur rapidly compared to the synthesis and degradation (see Figure 2C ). Based on this, we introduce additional steps in the chain, specifically for the synthesis and degradation processes of Cyc. Instead of modeling each process as a single step, we approximate them as n discrete reversible steps, each with identical forward and backward transition times. The resulting Markov chain is illustrated in Figure 2B . Again, we can map one period of the chain to the a linear cascade of forward and backward reactions such as in Figure S1B . This chain can be seen as two consecutive cascades as in figure S1D . Therefore its MFP time can be calculated by applying Eq. (5 ) to each of the individual chains and summing the results. In the continuous limit n → ∞, the chain should approximate the dynamics of the original ODE model. From Eq. (26 ) applied to the two consecutive chains, we find: We used this expression to fit the MFP time of the Markov chain to the period of the ODE model as a function of temperature by tuning the unknown pre-exponential constants. A comparison is shown in Figure 2F . Once these constants are determined, we can compare the MFP time for finite values of n to the ODE model period. We find that n = 50 already gives a close approximation within the temperature range where the ODE model exhibits a limit cycle. Only near the thermal boundaries, where the limit cycle is lost, that is when and when (see Eq. (27 )), does the MFP time fail to fully capture the temperature dependence. This alternative formulation of the cell cycle model as a Markov chain offers new insights into the period increase and eventual breakdown of oscillatory dynamics in thermal limits. The MFP time for Cyc synthesis follows Eq. (5 ), with and similarly for degradation If E d ≠ E s , both ratios are exponentially temperature-dependent. Within the functional temperature range of the oscillator, both synthesis and degradation must be forward-directed, meaning r s , r d < 1 and As the temperature changes and r d → 1, the MFP time of degradation increases dramatically and the total MFP time diverges. Physically, the system becomes trapped in forward-backward loops during degradation, and the limit cycle is destroyed. Likewise, if r s → 1, backward cycling dominates during synthesis, again disrupting the periodic behavior. Since r s and r d have opposite temperature dependencies, one of these breakdowns will occur in the high-temperature limit and the other in the low-temperature limit. This behavior is illustrated in Figure 2G . Supplementary Note 4: Temperature scaling of a sum of random transition times: Arrhenius vs. quadratic exponential In many cases, the MFP time of a biological system can be approximated as a the sum of a large number, N , of (random) effective transition times which scale Arrhenius-like around some reference temperature T 0 . Examples given in this paper include Eq. (5 ) and Eq. (42 ). In this section we derive the temperature scaling of such a system and show that two possible forms of scaling emerge at different temperature ranges: QE scaling when we can average over a large number of effective transition times and Arrhenius scaling when one of the transition times dominates. A. Distribution from the central limit theorem Assuming forward propagation r ≪ 1 in the cascade and that preexponential factors A i and activation energies of forward reactions E i are random variables drawn from arbitrary distributions, we can rewrite Eq. (5 ) as In this expression, we took β = 1 / ( RT ) −1 / ( RT 0 ) such that A i is the effective transition time at the reference temperature. In the limit N → ∞, the central limit theorem assures that if all exponentials are independent and identically distributed random variables, then 𝔼 [𝒯] will be normally distributed with mean µ and variance σ 2 over different realizations, where E i and A i ( i = 1,.., N ) are randomly drawn from their respective distributions where where ⟨…⟩ denotes the mean with respect to the distributions of A i and E i ( i = 1,.., N ). Rewriting Eq. (34 ) as ,we express 𝔼 [ A i ] = µ A and use the cumulant form for : where κ 1 = µ E and κ 2 = Var[ E ]. Neglecting cumulants of order 3 and higher, one gets a closed-form expression for the MFP time: Note that this general form is valid for arbitrary distributions of A i and E i . Likewise, we find the variance of MFP time: At extremal temperatures the variance in A i at the reference temperature can usually be neglected which means that B. Relative variation We note that µ from Eq. (37 ) has a QE temperature dependence. Whether 𝔼 [𝒯] will follow this temperature scaling of µ for a typical realization drawn from the distributions of A i and E i , depends on the relative variation In the region where is a good approximation of 𝔼 [𝒯] for all realizations, i.e., simulations of a system each with a given set of ( A i , E i ), thus the system scales according to a QE law. If , which is the case for both hot and cold temperature extremes, energy variation Var contributes to a of the MFP time. The picture that emerges is that of two different regimes with a boundary at extremal temperatures that is located at σ/µ = 1 or at Supplementary Note 5: Linear cascades of reversible processes with random transition times In this section, we approximate biological systems as linear cascades of reversible processes with n distinct transition steps, as schematically illustrated in Figure S1E . A. The mean first passage time The MFP time for such a cascade is given by where the are the effective forward transition times, defined as Each represents the average time required to move from state 1 to state i accounting for all possible backward excursions along the cascade, which we refer to as an effective passage time. These effective times reduce the reversible cascade to an equivalent process with purely forward transitions. B. Forwards and backwards regime By analogy with Eq. (21 ), we define ⟨ r ⟩ = ⟨ τ i,f /τ i,b ⟩. If ⟨ r ⟩ ≪ 1, it suggests that that on average forward transitions dominate the propagation in the cascade, so we can approximate and from Eq. (42 ) we find In this forwards regime, we can ignore the backwards transitions in the computation of the MFP time and the system behaves as if there are only forward transitions. If⟨ r ⟩~ 1 the backward transitions will increase the MFP time. Physically, this happens because the system can get trapped in backward cycles before it exits towards the final state in the forward direction. If ⟨ r ⟩≫ 1, we say we are in the backwards regime. The boundary between the two regimes is determined by the temperature at which ⟨ r ⟩ = 1. C. Introducing temperature dependence Connecting the results in the previous subsections with the rest of the paper, we again assume all transition times have an Arrhenius dependence on temperature In this, T 0 is a reference temperature at which we suppose the cascade is working in the forward regime and β ( T ) = 1 / ( RT ) −1 / ( RT 0 ) is introduced to facilitate calculations. We now choose a distribution of transition times by imposing a normal probability distribution 𝒩 ( µ, σ ) with mean µ and standard deviation σ 2 on the activation energies and the the transition times at the reference temperature: Note that all activation energies are assumed to be identically distributed, but the forward and the backward transition times have a different mean at the reference temperature which is controlled by . In the forward working regime we expect that Δ < 0. D. Temperature-dependent distributions of forward and backward transition times Assuming that the activation energies as well as the logarithms of forward and backward transition times at the reference temperature T 0 are drawn from the normal distributions, Eqs. (46 )-( 48 ), distri-bution of τ { f,b } at temperature T is given by: where µ f = µ 0 and µ b = µ 0 − Δ are means of the forward and backward transition times at T 0 . Next, one can derive the distribution of the forward-to-backrward ratio r = τ f /τ b at a given temperature T as: Or, explicitly: From Eqs. (50 ) and ( 51 ), it follows that only the variance of the forward-to-backward ratio distribution is temperature-dependent, and, importantly, it increases symmetrically relative to the reference temperature T 0 ( Figure 3D ). In other words, as the system moves away from the reference temperature, the probability of getting r i > 1 increases on both cold and warm sides, amplifying the cycle motifs in the linear transition cascade. E. Temperature-dependent distributions of transition times in different populations Moving forward, we want to get more insight into the development of a heavy tail of the distribution of r away from T 0 . We take a closer look on the distribution of transition times depending on the balance between respective forward and backward activation energy. Taking into account that in the populations with E i,f > E i,b and with E i,f < E i,b transition times scale differently, i.e., backward transitions dominate at low temperatures in the former case and at the higher in the later case ( Figure 3C ), let us treat these two populations separately and find the distributions of conditioned variables τ { f,b } | E f > E b and τ { f,b } | E f E b and E { f,b } | E f < E b ( Figure. S2A ). To do so, consider a bivariate distribution ( E f , E b ), where each component individually is drawn from the normal distribution in Eq. (46 ). Next, define a vector: and set ρ = 0, assuming no correlation between E f and E b . To facilitate the analysis, introduce new variables: New variables ( U, V ) are obtained by linear transformation (rotation) of ( E f , E b ), and, thus, are also drawn from a bivariate normal distribution. From the latter it follows that: Now, without the loss of generality, consider the distribution of E f conditioned by E f > E b , or, equivalently, by U > 0. Defining E f through the new variables as E f = ( U + V ) / 2, we express sought conditioned variable as: where U | U > 0 is drawn from a truncated normal distribution with mean and variance: Substituting Eqs. (53 ), ( 55 ), ( 57 ), and ( 58 ) in Eq. (56 ), we obtain the distribution of forward activation energies E f in the population of transitions with E f > E b : Repeating the same for E b | E f > E b and for the E f E b and E { f,b } E f < E b , we find that sought transition times are also normally distributed: with respective means and variances: From the above equations, it follows that the variation of tran-sition times have the same temperature dependence in both populations, while the means scale differently: the forward and backward transition times get closer at lower temperatures and diverge at higher temperatures in the E f > E b population and the opposite occurs in the E f E b transitions amplify the cycles, while E f < E b transitions contribute to the cycles at warmer temperatures ( Figure S2B ). F. Temperature scaling Let us now take a look at the scaling of the MFP time. Based on the analysis of the forward and backward transition times in the previous section we distinguish between a forward regime where ⟨ r ( T ) ⟩ ≫1 and the MFP time is dominated by forward transitions, and a backward regime where ⟨ r ( T ) ⟩ ≫1 and the MFP time is dominated by the effect of backward transitions. Since for the lognormal distribution in Eq. (51 ) we have we find that the boundaries between the forward and the backward regime are located at Since the temperature dependence is exponential, the transition between the forwards and backwards regime will happen over a narrow temperature range. F.1. Temperature scaling in the forwards regime In S4 we derive the temperature scaling of a system which has a MFP time that can be written as the sum of exponentials such as in Eq. (44 ). The conclusion is that when the number of exponentials is large and σ β is small, the MFP time scales according to a QE temperature law that emerges from an averaging over all exponentials. Using the results in Eq. (34 ) and Eq. (41 ) and applying them to the forwards approximation in Eq. (44 ) where σ β = σ E and N = n , we find that in the forwards regime for all temperatures β for which At more extreme temperatures, individual transition times tend to dominate the sum in Eq. (44 ) and the system will scale Arrhenius-like, following the dominant forward transition time. However, if n is sufficiently large enough then and the MFP time scales quadratically within the entire forward regime. F.2. Temperature scaling in the backwards regime To understand what the temperature dependence of the MFP time in the temperature region ⟨ r ⟩ > 1 looks like, we can start from Eq. (43 ) and note that with Eq. (45 ) the effective transition times take again the form of a sum of exponentials with and A i,k a constant factor that is determined as the product and quotient of pre-exponential factors of individual transition times. The random activation energies E i,k can be positive or negative. The variance on these effective activation energies is much higher than that on the individual forward activation energies since This means that while the number of effective transition times is We conclude that and that if n → ∞. Thus, in the large-number limit, there will be no region of QE scaling where we can average over a large number of transition times and the MFP time always scales Arrhenius-like in the backwards regime. The exponential with will dominate the MFP time in the high- β or low-temperature limit and cause a positive linear slope in the Arrhenius diagram. The exponential with will dominate in the high-temperature limit and will typically cause a negative linear slope in the Arrhenius diagram. G. Numerical simulations To check these theoretical results we conducted a series of numerical simulations of linear cascades with transition times drawn from a random distribution. For a first round of simulations we used the parameter values in Table S2 (middle column) and n variable and used these parameters to assign temperature-dependent random transition times to each of the transitions in the linear cascade via Eqs. (46 ), ( Eq. (47 )), and ( 48 ). The motivation for the parameters of the activation energies comes from ( 44 ) which gives a survey of a large number of enzymatic reactions and their temperature dependence. Afterwards we calculated the MFP time of this cascade in function of temperature via Eq. (42 ) for a temperature range of · 1000/ K . In the QE temperature region we compared the MFP time with the theoretical prediction Eq. (70 ). Outside of this region we determined a best fit of a linear Arrhenius function with positive activation energy in the low-temperature and negative activation energy in the high-temperature region to the MFP time. This process was repeated 100 times for different values of n . The results are given and discussed in Figure 3F and the accompanying text. We also repeated this procedure for a uniform distribution of activation energies with the same mean and standard deviation as the normal distribution and reported no significant differences, confirming the reliability of the approximation via a normal distribution. In a second round of simulations we repeated the same procedure but now in a setup with a feedforward loop. The exact parameters that were used are given in Table S2 (right column). The network topology is the same as in the first row but contains an additional shortcut as shown in Figure 3B , now contains a feedforward loop with fixed transition times that has been engineered to become active at low temperatures. We did this by making E control ≪ µ E such that at low temperatures τ control ≪ τ i,f for an average τ i,f . Therefore, the shortcutted path becomes much more probable at low temperatures and will always be taken. The results for the simulations are given in Figure 3E and G and discussed in the accompanying text. Supplementary Note 6: Linear cascade of backward cycles with random transition times In this section we derive the temperature scaling of the network topology in Fig. 6D and show how it can lead to negative QE temperature scaling in a certain temperature range. A. The mean first passage time The MFP time is the sum of the MFP time of n − 1 individual cycles (see Eq. (92 )) B. Temperature scaling As before we call ,We can again distinguish between a forward regime when ⟨ r ⟩ ≫ 1 and a backwards regime when ⟨ r ⟩ ≪ 1. B.1. Forward regime In the forward regime, the contribution from the backward cycles to the temperature scaling can be neglected. The MFP time is the sum of 2( n − 1) forward transition times. We assume a normal distribution for the forward activation energies and a mean of µ 0, f for the forward pre-exponential factors. From Eq. (34 ) we find that in the large number limit For finite-sized systems we expect to see deviations from the QE in the form of Arrhenius scaling at cold temperatures corresponding to B.2. Backwards regime In the backward regime the MFP time will be dominated by the contributions from the backward cycles and we can neglect the individual forward reactions in Eq. (77 ) to find We assume a normal distribution and for the forward and backward activation ener-gies and a means µ 0, f and µ 0, b for the forward and bakcward pre-exponential factors. Again using, Eq. (34 ) we find that in the large number limit a QE scaling law emerges If µ E,b > 2 µ E,f then this will be a negative QE. For finite-sized systems we expect to see deviations from the QE in the form of Arrhenius scaling at hot temperatures corresponding to Typically, it will be the individual exponentials with large negative activation energies that will dominate the MFP time at these extremal temperatures. C. Total temperature scaling Suppose that we are dealing with a finite-sized system, that r ( β * ) = 1 and that µ E,b > 2 µ E,f . Based on the discussion above we can then discertain between four different regions of temperature scaling Arrhenius scaling at the hot side when Negative QE scaling at the hot side with a large curvature when QE scaling with a smaller curvature at the cold side when Arrhenius scaling at the cold side when D. Generalizations In the examplatory topology of this section, the difference between the forward and backward activation energies should quite high to get a negative QE in the bakcwards regime. However, one could also consider cascades of cycles with m forward and backward reactions where m > 1. In the backwards regime, the largest exponentials would dominate and the approximate MFP time would become In this case one would get a negative QE for which makes that the difference between the forward and backward activation energies need not be large to get a negative slope. Note that the curvature also increases for larger backward cycles. If multiple backward cycles with different lengths would be present in the same cascade, the cycles with the longest length would always dominate in the backward regime and the other ones could be ignored. Supplementary Note 7: Simple motifs with a negative activation energy One of the main paradoxes to be solved, is how the overall negative activation energy emerges from a series of small-scale transitions which all have positive activation energies. In this section, we discuss one example that illustrates when this is not possible, and two minimal examples of state graph topologies for which this is possible. The two last motifs highlight the two physically different mechanisms that might lead to an overall macroscopic negative activation energy in the temperature scaling of the MFP time. A. Cascade of sequential forward reactions Consider a simple cascade of forward-directed transitions such as depicted in Figure S1C . The MFP time is just the sum of all individual transition times If we calculate the temperature dependence of 𝔼 [𝒯] we find that We conclude that the MFP time will never show a negative slope in the Arrhenius diagram if there is only one possible path towards the final state. An example with two forward transitions is shown in the bottom half of Figure 4B . To make this plot we used the parameter in Table S3 (left). B. A cycle Consider the simplest state graph with one cycle such as depicted in Figure 4D 4C ??. Calculating the MFP time, using Equation Eq. (42) and Eq. (43) , results in Besides the sum of the two forward transition times, an extra term appears. It captures the possibility of trapping the system in the cycle (1) ↔( 2 ) a number of times before it exits towards state ( 3 ). Look at the particular case: E 21 > E 12 + E 23 . In the lowtemperature limit , systems scales Arrheniuslike with a positive activation energy E = max{ E 12 , E 23 }. But when the following scaling is observed However, does not correspond to physical temperatures. Therefore, it depends on the pre-exponential constants whether this limit is already observed for physical temperatures with T > 0. If this is the case, 𝔼 [𝒯] scales Arrhenius-like with a negative activation energy in the high-temperature limit. An example of this scaling behaviour can be seen in the bottom half of Figure 4C . To make this plot we used the parameter in Table S3 (middle). C. A feedforward loop Consider a minimal feedforward loop without cycles but with two forward directed paths such as in Figure 4E . The MFP time is given by Where P [1 → 2] is the chance that the path along state ( 2 ) is taken and P [1 → 3] is the chance that the path along state ( 3 ) is taken. Take E 12 < E 13 and τ 12 , τ 13 ≪ τ 24 , τ 34 . This setup makes that the followed path, is controlled by the first transition step, but the total MFP time is entirely dominated by the second step. In the low-temperature limit τ 12 ≪ τ 13 . Therefore the path along state ( 2 ) is taken and 𝔼 [𝒯] ≈ τ 24 . In the high-temperature limit, the reverse could happen which would result in 𝔼 [𝒯] ≈ τ 34 . Note that in both of these limits the local activation energy is constant and positive. However, if the transition between these two limits happens fast and τ 34 ≪ τ 24 the local activation energies in the transition between the two limits can be negative. An example for which this happens is shown in the bottom half of Figure 4D . To make this plot we used the parameters in Table S3 (right). Supplementary Note 8: Data analysis A. Dataset construction We conducted a literature search to find studies measuring the timing of various biological processes under temperature variation. We included recordings of durations of early and later embryonic development stages in different model organisms, as well as bacteria and phytoplankton growth rates. When available, we extracted data from public online repositories or article tables; otherwise, we used PlotDigitizer ( https://plotdigitizer.com ) to get data from the article figures. Overall, we collected 121 datasets from 48 source, resulting in a total amount of 2174 data points (18 points per set on average, minimum of 5 points per set). On average, temperature range spanned 30±17 K with the minimum of 10 K . B. Fitting procedure Fitting procedure aimed to select and fit the best model explaining temperature scaling of biological timing for individual curves. The candidate models were: Quadratic-exponential (QE) , which supposedly describe scaling of biological timing in a wide range of temperatures following Markov-chain formalism (See Eq. (7 )): with and ; Quadratic-linear-exponential (QLE) , capturing deviation of biological timing from the QE model at either cold or hot temperature extreme: Quadratic-double-linear-exponential (QDLE) , capturing deviation of biological timing from the QE model at both cold and hot temperature extremes: In both QLE and QDLE models, E { c,h } = RA { c,h } . To find and fit the best model to the data, we used a multi-step fitting procedure described below. Step 1: quadratic function was fitted on the plane (1000 /T , ln τ ( T )) in a whole temperature range to the all datasets. Step 2: datasets were labeled as QE based on the goodness-of-fit evaluation ( R 2 > 0.95) and additional visual inspection for less clean data. Step 3: datasets, exhibiting deviations from the QE function, were subjected to automated model fitting/selection process. Deviations were identified by linear regression on residuals of a quadratic fit (20% of datapoints at each side). If linear regression coefficient exceeded a predefined threshold, the deviation was identified. Based on the identified deviation (on the hot side, on the warm side, or both) either one-sided QLE or two-sided QDLE function was fitted. Step 4: similar to Step 2 , datasets were labeled as QLE or QDLE based on the goodness-of-fit evaluation ( R 2 > 0.95) and additional visual inspection for less clean data. Step 5: datasets, failing the automated fitting, were subjected to a similar model/fitting selection but with the manual setting of deviation boundaries. Parameter fitting was performed using the least squares algorithm implemented in the least_squares function from scipy.optimize package for Python. Goodness-of-fit was evaluated using the r2_score function from sklearn.metrics package for Python. Fitting code was optimized using Generative AI (ChatGPT-5o). Supplemental Tables View this table: View inline View popup Download powerpoint Table S1. All parameters which were used to simulate the 2 ODE cell-cycle model in Figure 2C,D,E,F . All parameters in the upper half are taken from ( 12 ) where they are shown to represent aspects of a realistic cell division cycle. The parameters are the values of the rate constants k i at T = T 0 . The activation energies in the bottom half are chosen by us to highlight specific aspects of the model. View this table: View inline View popup Download powerpoint Table S2. All parameters which were used to simulate the linear cascade of reversible processes in Figure 3 . The left column corresponds to panels D and F, the right column to panels E and G. View this table: View inline View popup Download powerpoint Table S3. All parameters which were used to simulate the three examples in Figure 4 . The left, middle and right column correspond to panels B, C and D. Supplemental Figures Download figure Open in new tab Fig. S1. Different linear cascades of reactions which were considered to explain the temperature dependence of biological networks. (A) One period of the minimal Markov approximation of the two-ODE cell cycle oscillator mapped to a linear cascade (Section A) (B) One period of the Markov approximation of the two-ODE cell cycle oscillator with additional states (Section C). (C) Cascade of forward reactions. This cascade has no macroscopic negative activation energy (Section B). (D) Cascade of fixed reversible reactions followed by a forward reaction. This cascade can produce a strong temperature dependence with negative activation energies in the high-temperature limit (Section C). (E) This cascade is the most general one and can produce negative activation energies in the high-temperature limit and a QE scaling at intermediate temperatures (section A) Download figure Open in new tab Fig. S2. Temperature-dependent distributions of forward and backward transition times. (A) Scatterplot showing bivariate distribution of E f versus E b . Histograms to the top and to the right show distributions of E f and E b conditioned by E f > E b (green) and E f E b (top row) and E f < E b (bottom row) at different temperatures. Download figure Open in new tab Fig. S3. The mean first passage time varies little in the intermediate temperature range. Simulation of the temperature dependence of the MFP time for system (A) and in Fig. 3 . Individual transition times are drawn from a distribution with µ E = 50 kJ/mol, σ E = 10 kJ/mol, µ 0 = −20, σ 0 = 0.2, µ r = 2. In the intermediate temperature regime, the QE prediction from Equation ( 7 ) provides a good fit. Boundaries of the QE regime are determined from Equation ( ?? ). (A) The MFP time as a function of temperature for 100 realizations of transition rate distributions. (B) The relative variation (standard deviation/mean) of MFP time across realizations as a function of temperature. ACKNOWLEDGEMENTS L.G. acknowledges funding by the KU Leuven Research Fund (grant number C14/23/130 and IDN/25/007) and the Research-Foundation Flanders (FWO, grant number G074321N). We thank Stephan Grill for providing us with the data from Ref. ( 13 ), and Daniel Cebrián-Lacasa for feedback on the manuscript. Funder Information Declared KU Leuven, https://ror.org/05f950310 , C14/23/130 , N/25/007 Research Foundation - Flanders, https://ror.org/03qtxy027 , 074321N Footnotes https://gitlab.kuleuven.be/gelenslab/publications/Jacobs_Beyond_Arrhenius.git ↵ 1 Here Q 10 ( T ) ≡ κ ( T + 10 ° C) /κ ( T ) is the local fold-change in rate per 10 ° C. For Arrhenius κ ( T ) = A exp(− E a /RT ), so Q 10 depends on both T and E a (e.g., at 20 ° C, E a = 50–80 kJ mol −1 gives Q 10 ≈ 2.0–2.8). For durations τ = 1 /κ , the corresponding factor is . In non-Arrhenius (curved) regimes, a single Q 10 is only a local approximation ( 12 ); for a general Δ T , the factor is . ↵ 2 Irreversible forward transitions can also be incorporated by setting the corresponding backward transition time to infinity. ↵ 3 Note that edges are often labeled by their transition time τ ij . This implicitly refers to a transition rate of . References 1. ↵ A Cossins . Temperature biology of animals . Springer , 1987 . 2. ↵ JM Sunday , AE Bates , MR Kearney , RK Colwell , NK Dulvy , JT Longino , and RB Huey . Thermal-safety margins and the necessity of thermoregulatory behavior across latitude and elevation . Proceedings of the National Academy of Sciences , 111 ( 15 ): 5610 – 5615 , 2014 . OpenUrl Abstract / FREE Full Text 3. ↵ Alexander C Ruane . Synthesis report of the ipcc sixth assessment report (ar6) . 2024 . 4. ↵ JF Gillooly , JH Brown , GB West , M Van Savage , and EL Charnov . Effects of size and temperature on metabolic rate . Science , 293 ( 5538 ): 2248 – 2251 , 2001 . OpenUrl Abstract / FREE Full Text 5. ↵ S Gilbert . Developmental biology . 2000 . 6. ↵ Svante Arrhenius . Über die reaktionsgeschwindigkeit bei der inversion von rohrzucker durch säuren . Zeitschrift für physikalische Chemie , 4 ( 1 ): 226 – 248 , 1889 . OpenUrl 7. ↵ SR Logan . The origin and status of the arrhenius equation . Journal of Chemical Education , 59 ( 4 ): 279 , 1982 . OpenUrl 8. ↵ KJ Laidler and MC King . Development of transition-state theory . The Journal of Physical Chemistry , 87 ( 15 ): 2657 – 2664 , 1983 . OpenUrl CrossRef Web of Science 9. ↵ JF Gillooly , E. Charnov, GB West , M Van Savage, and JH Brown . Effects of size and temperature on developmental time . Nature , 417 ( 6884 ): 70 – 73 , 2002 . OpenUrl CrossRef PubMed Web of Science 10. ↵ AI Dell , S Pawar, and VM Savage . Systematic variation in the temperature dependence of physiological and ecological traits . Proceedings of the National Academy of Sciences of the United States of America , 108 ( 26 ): 10591 – 10596 , 2011 . OpenUrl Abstract / FREE Full Text 11. ↵ J Faber and PD Nieuwkoop . Normal table of Xenopus laevis (Daudin): A systematical and chronological survey of the development from the fertilized egg till the end of metamorphosis. Garland Science , 1994 . 12. ↵ Jan Rombouts , Franco Tavella , Alexandra Vandervelde 1, Connie Phong , James E. Ferrell Jr. , Qiong Yang , and Lendert Gelens . Mechanistic origins of temperature scaling in the early embryonic cell cycle . Nature Communications , 16 ( 8045 ): 1 – 18 , 2025 . OpenUrl CrossRef PubMed 13. ↵ ML Begasse , M Leaver , F Vazquez , SW Grill , and AA Hyman . Temperature dependence of cell division timing accounts for a shift in the thermal limits of c. elegans and c. briggsae . Cell Reports , 10 ( 5 ): 647 – 653 , 2015 . OpenUrl PubMed 14. ↵ J Crapse , N Pappireddi , M Gupta , SY Shvartsman , E Wieschaus , and M Wühr . Evaluating the arrhenius equation for developmental processes . Molecular Systems Biology , 17 ( 8 ): e9895 , 2021 . OpenUrl CrossRef PubMed 15. ↵ Louis Powsner . The effects of temperature on the durations of the developmental stages of drosophila melanogaster . Physiological Zoology , 8 ( 4 ): 474 – 520 , 1935 . OpenUrl CrossRef 16. ↵ JP DeLong , JP Gibert , TM Luhring , G Bachman , B Reed , A Neyer , and KL Montooth . The combined effects of reactant kinetics and enzyme stability explain the temperature dependence of metabolic rates . Ecology and Evolution , 7 ( 11 ): 3940 – 3950 , 2017 . OpenUrl 17. RM Daniel , MJ Danson , and R Eisenthal . The temperature optima of enzymes: A new perspective on an old phenomenon . Trends in Biochemical Sciences , 26 ( 4 ): 223 – 225 , 2001 . OpenUrl CrossRef PubMed Web of Science 18. ↵ Julian B. Voits and Ulrich S. Schwarz . The generic temperature response of large biochemical networks . arXiv, (2403.17202) , 2024 . 19. ↵ J Newport and M Kirschner . A major developmental transition in early xenopus embryos: characterization and timing of cellular changes at the midblastula stage . Cell , 30 ( 3 ): 675 – 686 , 1982 . OpenUrl CrossRef PubMed Web of Science 20. ↵ J Newport and M Kirschner . A major developmental transition in early xenopus embryos: Ii. control of the onset of transcription . Cell , 30 ( 3 ): 687 – 696 , 1982 . OpenUrl CrossRef PubMed Web of Science 21. ↵ DA Ratkowsky , J Olley , and T Ross . Unifying temperature effects on the growth rate of bacteria and the stability of globular proteins . Journal of Theoretical Biology , 233 ( 3 ): 351 – 362 , 2005 . OpenUrl CrossRef PubMed Web of Science 22. JI Arroyo , B Díez , CP Kempes , GB West , and PA Marquet . A general theory for temperature dependence in biology . Proceedings of the National Academy of Sciences , 119 ( 30 ): e2119872119 , 2022 . OpenUrl CrossRef PubMed 23. ↵ TL Wagner , RL Olson , and JL Willers . Modelling arthropod development time . 1991 . 24. ↵ Dimitris C Kontodimas , Panagiotis A Eliopoulos , George J Stathas , and Leonidas P Economou . Comparative temperature-dependent development of nephus includens (kirsch) and nephus bisignatus (boheman)(coleoptera: Coccinellidae) preying on planococcus citri (risso)(homoptera: Pseudococcidae): evaluation of a linear and various nonlinear models using specific criteria . Environmental Entomology , 33 ( 1 ): 1 – 11 , 2004 . OpenUrl CrossRef 25. ↵ Hendrik Anthony Kramers . Brownian motion in a field of force and the diffusion model of chemical reactions . physica , 7 ( 4 ): 284 – 304 , 1940 . OpenUrl CrossRef Web of Science 26. ↵ Peter Hänggi , Peter Talkner , and Michal Borkovec . Reaction-rate theory: fifty years after kramers . Reviews of modern physics , 62 ( 2 ): 251 , 1990 . OpenUrl CrossRef Web of Science 27. ↵ H Eyring . The activated complex in chemical reactions . The Journal of Chemical Physics , 3 ( 2 ): 107 – 115 , 1935 . OpenUrl CrossRef 28. ↵ Q Yang and JE Ferrell Jr . The cdk1-apc/c cell cycle oscillator circuit functions as a timedelayed, ultrasensitive switch . Nature Cell Biology , 15 ( 5 ): 519 – 525 , 2013 . OpenUrl CrossRef PubMed Web of Science 29. ↵ AW Murray and MW Kirschner . Cyclin synthesis drives the early embryonic cell cycle . Nature , 339 ( 6222 ): 275 – 280 , 1989 . OpenUrl CrossRef PubMed Web of Science 30. ↵ JR Pomerening , ED Sontag , and JE Ferrell Jr . Building a cell cycle oscillator: hysteresis and bistability in the activation of cdc2 . Nature Cell Biology , 5 ( 4 ): 346 – 351 , 2003 . OpenUrl CrossRef PubMed Web of Science 31. ↵ W Sha , J Moore , K Chen , AD Lassaletta , CS Yi , JJ Tyson , and JC Sible . Hysteresis drives cell-cycle transitions in xenopus laevis egg extracts . Proceedings of the National Academy of Sciences , 100 ( 3 ): 975 – 980 , 2003 . OpenUrl Abstract / FREE Full Text 32. ↵ Tom Evans , Eric T Rosenthal , Jim Youngblom , Dan Distel , and Tim Hunt . Cyclin: a protein specified by maternal mrna in sea urchin eggs that is destroyed at each cleavage division . Cell , 33 ( 2 ): 389 – 396 , 1983 . OpenUrl CrossRef PubMed Web of Science 33. ↵ Marie-Anne Félix , Jean-Claude Labbé , Marcel Dorée , Tim Hunt , and Eric Karsenti . Triggering of cyclin degradation in interphase extracts of amphibian eggs by cdc2 kinase . Nature , 346 ( 6282 ): 379 – 382 , 1990 . OpenUrl CrossRef PubMed Web of Science 34. ↵ B Novak and JJ Tyson . Modeling the cell division cycle: M-phase trigger, oscillations, and size control . Journal of Theoretical Biology , 165 : 101 – 134 , 1993 . OpenUrl CrossRef Web of Science 35. ↵ B Novák and JJ Tyson . Design principles of biochemical oscillators . Nature Reviews Molecular Cell Biology , 9 ( 12 ): 981 – 991 , 2008 . OpenUrl CrossRef PubMed Web of Science 36. TY Tsai , YS Choi , W Ma , JR Pomerening , C Tang , and JE Ferrell Jr . . Robust, tunable biological oscillations from interlinked positive and negative feedback loops . Science, 321 (5885):126–129 , 2008 . 37. Joseph R Pomerening , Sun Young Kim , and James E Ferrell . Systems-level dissection of the cell-cycle oscillator: bypassing positive feedback produces damped oscillations . Cell , 122 ( 4 ): 565 – 578 , 2005 . OpenUrl CrossRef PubMed Web of Science 38. J De Boeck , J Rombouts , and L Gelens . A modular approach for modeling the cell cycle based on functional response curves . PLOS Computational Biology , 17 ( 8 ): 1 – 39 , 2021 . OpenUrl CrossRef 39. P Parra-Rivas , D Ruiz-Reynés , and L Gelens . Cell cycle oscillations driven by two inter-linked bistable switches . Molecular Biology of the Cell , 34 ( 6 ): ar56 , 2023 . OpenUrl CrossRef PubMed 40. G Maryu and Q Yang . Nuclear-cytoplasmic compartmentalization of cyclin b1-cdk1 promotes robust timing of mitotic events . Cell Reports , 41 ( 13 ): 111870 , 2022 . OpenUrl PubMed 41. ↵ J Kamenz , L Gelens , and JE Ferrell Jr . Bistable, biphasic regulation of pp2a-b55 accounts for the dynamics of mitotic substrate phosphorylation . Current Biology , 31 ( 4 ): 794 – 808 , 2021 . OpenUrl CrossRef PubMed 42. ↵ Fangting Li , Tao Long , Ying Lu , Qi Ouyang , and Chao Tang . The yeast cell-cycle network is robustly designed . Proceedings of the National Academy of Sciences , 101 ( 14 ): 4781 – 4786 , 2004 . OpenUrl Abstract / FREE Full Text 43. ↵ J.R. Norris . Markov Chains . Cambridge Series on Statistical and Probabilistic Mathematics . Cambridge University Press , 1997 . ISBN 978-0-521-48181-6 . 44. ↵ Mikael Elias , Grzegorz Wieczorek , Shaked Rosenne , and Dan S Tawfik . The universality of enzymatic rate–temperature dependency . Trends in Biochemical Sciences , 39 ( 1 ): 1 – 7 , 2014 . OpenUrl CrossRef PubMed 45. ↵ Etienne Low-Décarie , Tobias G. Boatman , Noah Bennett , Will Passfield , Antonio Gavalás-Olea , Philipp Siegel , and Richard J. Geider . Predictions of response to temperature are contingent on model choice and data quality . Ecology and Evolution , 7 ( 23 ): 10467 – 10481 , 2017 . OpenUrl 46. James Gray . The mechanism of ciliary movement . Proceedings of the Royal Society of London. Series B, Containing Papers of a Biological Character , 93 ( 650 ): 104 – 121 , 1922 . OpenUrl CrossRef 47. Yoshimasa Yamamoto and Hiroyuki Nakahara . The formation and degradation of cyanobacterium aphanizomenon flos-aquae blooms: the importance of ph, water temperature, and day length . Limnology , 6 ( 1 ): 1 – 6 , 2005 . OpenUrl 48. Philip W Boyd , Tatiana A Rynearson , Evelyn A Armstrong , Feixue Fu , Kendra Hayashi , Zhangxi Hu , David A Hutchins , Raphael M Kudela , Elena Litchman , Margaret R Mulholland , et al. Marine phytoplankton temperature versus growth responses from polar to tropical waters–outcome of a scientific community-wide study . PloS one , 8 ( 5 ): e63091 , 2013 . OpenUrl CrossRef PubMed 49. C Scholtissek and R Rott . Effect of temperature on the multiplication of an influenza virus . Journal of General Virology , 5 ( 2 ): 283 – 290 , 1969 . OpenUrl PubMed 50. Yi-Hua Zheng , Jun Liu , Yong-Zhi Ma , Yi Xu , Fei Xu , and Tse-Chao Hua . Temperature effects on enzyme activity of chicken liver esterase used in calorimetric biosensor . Artificial Cells, Blood Substitutes, and Biotechnology , 40 ( 1-2 ): 125 – 131 , 2012 . OpenUrl 51. Leigh W Tait and David R Schiel . Impacts of temperature on primary productivity and respiration in naturally structured macroalgal assemblages . PLoS one , 8 ( 9 ): e74413 , 2013 . OpenUrl PubMed 52. Bryce J Brylawski and Thomas J Miller . Bioenergetic modeling of the blue crab (callinectes sapidus) using the fish bioenergetics (3.0) computer program . Bulletin of Marine Science , 72 ( 2 ): 491 – 504 , 2003 . OpenUrl 53. Ross Corkrey , Tom A McMeekin , John P Bowman , David A Ratkowsky , June Olley , and Tom Ross . Protein thermodynamics can be predicted directly from biological growth rates . PloS one , 9 ( 5 ): e96100 , 2014 . OpenUrl CrossRef PubMed 54. Kajal Gaur , Harish Kumar , Rathour PK Agarwal , KVS Baba , and SK Soonee . Analysing the electricity demand pattern . In 2016 National Power Systems Conference (NPSC) , pages 1 – 6 . IEEE , 2016 . 55. PW Terentév . Attempt at application of analysis of variation to the qualitative richness of the fauna of terrestrial vertebrates of the ussr vestnik leningradskovo universiteta 21 , 19 – 56. PW Terentév . Attempt at application of analysis of variation to the qualitative richness of the fauna of terrestrial vertebrates of the ussr vestnik leningradskovo universiteta 21, 19–26. translated by ej maly and edited by er pianka . Smithsonian Herpetological Information Services , 1968 . 57. Moriaki Yasuhara , Gene Hunt , Harry J Dowsett , Marci M Robinson , and Danielle K Stoll . Latitudinal species diversity gradient of marine zooplankton for the last three million years . Ecology letters , 15 ( 10 ): 1174 – 1179 , 2012 . OpenUrl CrossRef PubMed 58. Yang Ouyang , Jeanette M Norton , and John M Stark . Ammonium availability and temperature control contributions of ammonia oxidizing bacteria and archaea to nitrification in an agricultural soil . Soil Biology and Biochemistry , 113 : 161 – 172 , 2017 . OpenUrl 59. PD Franzmann , CM Haddad , RB Hawkes , WJ Robertson , and JJ Plumb . Effects of temperature on the rates of iron and sulfur oxidation by selected bioleaching bacteria and archaea: application of the ratkowsky equation . Minerals Engineering , 18 ( 13-14 ): 1304 – 1314 , 2005 . OpenUrl 60. GHJ Krüger and JN Eloff . The effect of temperature on specific growth rate and activation energy of microcystis and synechococcus isolates relevant to the onset of natural blooms . Journal of the Limnological Society of Southern Africa , 4 ( 1 ): 9 – 20 , 1978 . OpenUrl 61. Promita Deb , Saimon Ahmad Talukdar , Kaniz Mohsina , Palash Kumar Sarker , and SM Abu Sayem . Production and partial characterization of extracellular amylase enzyme from bacillus amyloliquefaciens p-001 . SpringerPlus , 2 ( 1 ): 154 , 2013 . OpenUrl PubMed 62. Georges Feller , Thierry Lonhienne , Christophe Deroanne , Cecille Libioulle , Jozef Van Beeumen , and Charles Gerday . Purification, characterization, and nucleotide sequence of the thermolabile alpha-amylase from the antarctic psychrotroph alteromonas haloplanctis a23 . Journal of Biological Chemistry , 267 ( 8 ): 5217 – 5221 , 1992 . OpenUrl Abstract / FREE Full Text 63. Christian Knoblauch and Bo Barker Jørgensen . Effect of temperature on sulphate reduction, growth rate and growth yield in five psychrophilic sulphate-reducing bacteria from arctic sediments . Environmental Microbiology , 1 ( 5 ): 457 – 467 , 1999 . OpenUrl CrossRef PubMed Web of Science 64. Irene H Tarpgaard , Antje Boetius , and Kai Finster . Desulfobacter psychrotolerans sp. nov., a new psychrotolerant sulfate-reducing bacterium and descriptions of its physiological response to temperature changes . Antonie Van Leeuwenhoek , 89 ( 1 ): 109 – 124 , 2006 . OpenUrl CrossRef PubMed Web of Science 65. Shinichi Sunagawa , Luis Pedro Coelho , Samuel Chaffron , Jens Roat Kultima , Karine Labadie , Guillem Salazar , Bardya Djahanschiri , Georg Zeller , Daniel R Mende , Adriana Alberti , et al. Structure and function of the global ocean microbiome . Science , 348 ( 6237 ): 1261359, 2015 . 66. Vania Battestin and Gabriela A Macedo . Effects of temperature, ph and additives on the activity of tannase produced by paecilomyces variotii . Electronic Journal of Biotechnology , 10 ( 2 ): 191 – 199 , 2007 . OpenUrl 67. Edward Mullaney , Kandan Sethumadhavan , Stephanie Boone , Xin Gen Lei , and Abul HJ Ullah . Elimination of a disulfide bridge in aspergillus niger nrrl 3135 phytase (phya) enhances heat tolerance and optimizes its temperature versus activity profile . Advances in Biological Chemistry , 2 ( 4 ): 372 – 378 , 2012 . OpenUrl 68. AK Rajasekaran and R Maheshwari . Thermophilic fungi: an assessment of their potential for growth in soil . Journal of biosciences , 18 ( 3 ): 345 – 354 , 1993 . OpenUrl CrossRef 69. RSS Yadav , KS Yadav , and HS Yadav . Luffa aegyptiaca (gourd) fruit juice as a source of peroxidase . Enzyme Research , 2011 ( 1 ): 319105 , 2011 . OpenUrl PubMed 70. ↵ HA Esechie . Interaction of salinity and temperature on the germination of sorghum . Journal of Agronomy and Crop science , 172 ( 3 ): 194 – 199 , 1994 . OpenUrl 71. SG Kuntz and MB Eisen . Drosophila embryogenesis scales uniformly across temperature in developmentally diverse species . PLOS Genetics , 10 ( 4 ): e1004293 , 2014 . OpenUrl 72. Chester I Bliss . Temperature characteristics for prepupal development in drosophila melanogaster . Journal of General Physiology , 9 ( 4 ): 467 – 495 , 1926 . OpenUrl Abstract / FREE Full Text 73. Gang Ma , Ary A Hoffmann , and Chun-Sen Ma . Daily temperature extremes play an imJacobs et al. | Beyond Arrhenius: Temperature scaling of biological rates. portant role in predicting thermal effects . The Journal of Experimental Biology , 218 ( 14 ): 2289 – 2296 , 2015 . OpenUrl Abstract / FREE Full Text 74. D Grzebyk and B Berland . Influences of temperature, salinity and irradiance on growth of prorocentrum minimum (dinophyceae) from the mediterranean sea . Journal of Plankton Research , 18 ( 10 ): 1837 – 1849 , 1996 . OpenUrl CrossRef Web of Science 75. Lina Chen , Yi Wei , Mao Shi , Zhengqun Li , and Shi-Hong Zhang . An archaeal chitinase with a secondary capacity for catalyzing cellulose and its biotechnological applications in shell and straw degradation . Frontiers in microbiology , 10 : 1253 , 2019 . OpenUrl PubMed 76. Qian Han , Tao Cai , Danilo A Tagle , Howard Robinson , and Jianyong Li . Substrate specificity and structure of human aminoadipate aminotransferase/kynurenine aminotransferase ii . Bioscience reports , 28 ( 4 ): 205 – 215 , 2008 . OpenUrl Abstract / FREE Full Text 77. Bo Barker Jørgensen , Mai F Isaksen , and Holger W Jannasch . Bacterial sulfate reduction above 100 c in deep-sea hydrothermal vent sediments . Science , 258 ( 5089 ): 1756 – 1757 , 1992 . OpenUrl Abstract / FREE Full Text 78. Anne E Taylor , Andrew T Giguere , Conor M Zoebelein , David D Myrold , and Peter J Bottomley . Modeling of soil nitrification responses to temperature reveals thermodynamic differences between ammonia-oxidizing activity of archaea and bacteria . The ISME journal , 11 ( 4 ): 896 – 908 , 2017 . OpenUrl PubMed 79. Göran Englund , Gunnar Öhlund , Catherine L Hein , and Sebastian Diehl . Temperature dependence of the functional response . Ecology letters , 14 ( 9 ): 914 – 921 , 2011 . OpenUrl CrossRef PubMed 80. Muhammad Rahim Khan and Muhammad Rafique Khan . The relationship between temperature and the functional response of coccinella septempunctata (l.)(coleoptera: Coccinellidae) . Pakistan Journal of Zoology , 42 ( 5 ), 2010 . 81. BN Dingha , AG Appel , and JT Vogt . Effects of temperature on the metabolic rates of insecticide resistant and susceptible german cockroaches, blattella germanica (l.)(dictyoptera: Blattellidae) . Midsouth Entomologist , 2 : 17 – 27 , 2009 . OpenUrl 82. Carol Hoole , Zenon J Czenze , Nigel Charles Bennett , and Andrew E McKechnie . Thermal physiology of three sympatric small mammals from southern africa . Journal of Zoology , 307 ( 1 ): 28 – 35 , 2019 . OpenUrl CrossRef 83. Maksim Zakhartsev , Xuelian Yang , Matthias Reuss , and Hans Otto Pörtner . Metabolic efficiency in yeast saccharomyces cerevisiae in relation to temperature dependent growth and biomass yield . Journal of thermal biology , 52 : 117 – 129 , 2015 . OpenUrl CrossRef PubMed 84. Pascal Hingamp , Nigel Grimsley , Silvia G Acinas , Camille Clerissi , Lucie Subirana , Julie Poulain , Isabel Ferrera , Hugo Sarmento , Emilie Villar , Gipsi Lima-Mendez , et al. Exploring nucleo-cytoplasmic large dna viruses in tara oceans microbial metagenomes . The ISME journal , 7 ( 9 ): 1678 – 1695 , 2013 . OpenUrl PubMed 85. Christy M McCain and Nathan J Sanders . Metabolic theory and elevational diversity of vertebrate ectotherms . Ecology , 91 ( 2 ): 601 – 609 , 2010 . OpenUrl CrossRef PubMed Web of Science 86. Mathias Milici , Jürgen Tomasch , Melissa L Wos-Oxley , Hui Wang , Ruy Jáuregui , Amelia Camarinha-Silva , Zhi-Luo Deng , Iris Plumeier , Helge-Ansgar Giebel , Mascha Wurst , et al. Low diversity of planktonic bacteria in the tropical ocean . Scientific reports , 6 ( 1 ): 19054 , 2016 . OpenUrl PubMed 87. Sharon M Talley , Phyllis D Coley , and Thomas A Kursar . The effects of weather on fungal abundance and richness among 25 communities in the intermountain west . BMC ecology , 2 ( 1 ): 7 , 2002 . OpenUrl PubMed 88. ↵ Erwan G Roussel , Barry A Cragg , Gordon Webster , Henrik Sass , Xiaohong Tang , Angharad S Williams , Roberta Gorra , Andrew J Weightman , and R John Parkes . Complex coupled metabolic and prokaryotic community responses to increasing temperatures in anaerobic marine sediments: critical temperatures and substrate changes . FEMS Microbiology Ecology , 91 ( 8 ): fiv084 , 2015 . OpenUrl CrossRef PubMed 89. ↵ Haizhou Li , Qunhui Yang , Jian Li , Hang Gao , Ping Li , and Huaiyang Zhou . The impact of temperature on microbial diversity and aoa activity in the tengchong geothermal field, china . Scientific Reports , 5 ( 1 ): 17056 , 2015 . OpenUrl PubMed 90. ↵ Christian Schröter , Leah Herrgen , Albert Cardona , Gary J Brouhard , Benjamin Feldman , and Andrew C Oates . Dynamics of zebrafish somitogenesis . Developmental Dynamics , 237 ( 3 ): 545 – 553 , 2008 . OpenUrl CrossRef PubMed Web of Science 91. ↵ Gabriella Martini and Hernan Garcia . Temperature variation drives coordinated scaling of temporal and dynamic features of transcription in embryonic development . bioRxiv , 2023.08.09.552666, 2023 . 92. ↵ MW Dorrity , LM Saunders , M Duran , SR Srivatsan , E Barkan , DL Jackson , SM Sattler , B Ewing , C Queitsch , J Shendure , DW Raible , D Kimelman , and C Trapnell . Proteostasis governs differential temperature sensitivity across embryonic cell types in zebrafish . Cell , 186 ( 23 ): 5015 – 5027 .e12, 2023 . OpenUrl CrossRef PubMed 93. ↵ Jess J Bourn and Michael W Dorrity . Degrees of freedom: temperature’s influence on developmental rate . Current Opinion in Genetics & Development , 85 : 102155 , 2024 . OpenUrl PubMed 94. ↵ H Falahati , W Hur , S Di Talia , and E Wieschaus . Temperature-induced uncoupling of cell cycle regulators . Developmental Biology , 470 : 147 – 153 , 2021 . OpenUrl CrossRef PubMed 95. ↵ Anatoly Rinberg , Adam L Taylor , and Eve Marder . The Effects of Temperature on the Stability of a Neuronal Oscillator . PLOS Computational Biology , 9 ( 1 ): e1002857 , jan 2013 . OpenUrl PubMed 96. ↵ Lamont S Tang , Marie L Goeritz , Jonathan S Caplan , Adam L Taylor , Mehmet Fisek , and Eve Marder . Precise Temperature Compensation of Phase in a Rhythmic Motor Pattern . PLOS Biology , 8 ( 8 ): e1000469 , aug 2010 . OpenUrl CrossRef PubMed 97. Colin S Pittendrigh . On temperature independence in the clock system controlling emergence time in drosophila . Proceedings of the National Academy of Sciences , 40 ( 10 ): 1018 – 1029 , 1954 . OpenUrl FREE Full Text 98. J Woodland Hastings and Beatrice M Sweeney . On the mechanism of temperature independence in a biological clock . Proceedings of the National Academy of Sciences , 43 ( 9 ): 804 – 811 , 1957 . OpenUrl FREE Full Text 99. ↵ Z. Huang , K. Curtin , and M. Rosbash . PER protein interactions and temperature compensation of a circadian clock in Drosophila . Science , 267 ( 5201 ): 1169 – 1172 , feb 1995 . ISSN 0036-8075 . doi: 10.1126/science.7855598 . OpenUrl Abstract / FREE Full Text 100. ↵ Philip B Kidd , Michael W Young , and Eric D Siggia . Temperature compensation and temperature sensation in the circadian clock . Proceedings of the National Academy of Sciences , 112 ( 46 ): E6284 – E6292 , 2015 . OpenUrl Abstract / FREE Full Text 101. ↵ Peter Ruoff , Ludger Rensing , Roald Kommedal , and Saadat Mohsenzadeh . Modeling temperature compensation in chemical and biological oscillators . Chronobiology international , 14 ( 5 ): 499 – 510 , 1997 . OpenUrl CrossRef PubMed Web of Science 102. ↵ Peter Ruoff , Melinda K Christensen , Jana Wolf , and Reinhart Heinrich . Temperature dependency and temperature compensation in a model of yeast glycolytic oscillations . Biophysical chemistry , 106 ( 2 ): 179 – 192 , 2003 . OpenUrl CrossRef PubMed Web of Science 103. ↵ Jeroen S van Zon , David K Lubensky , Pim RH Altena , and Pieter Rein ten Wolde . An allosteric model of circadian kaic phosphorylation . Proceedings of the National Academy of Sciences , 104 ( 18 ): 7420 – 7425 , 2007 . OpenUrl Abstract / FREE Full Text 104. ↵ Tetsuhiro S Hatakeyama and Kunihiko Kaneko . Generic temperature compensation of biological clocks by autonomous regulation of catalyst concentration . Proceedings of the National Academy of Sciences , 109 ( 21 ): 8109 – 8114 , 2012 . OpenUrl Abstract / FREE Full Text 105. Lili Wu , Qi Ouyang , and Hongli Wang . Robust network topologies for generating oscillations with temperature-independent periods . PLOS ONE , 12 ( 2 ): e0171263 , feb 2017 . ISSN 1932-6203 . doi: 10.1371/journal.pone.0171263 . OpenUrl CrossRef PubMed 106. ↵ Christian I. Hong , Emery D. Conrad , and John J. Tyson . A proposal for robust temperature compensation of circadian rhythms . Proceedings of the National Academy of Sciences of the United States of America , 104 ( 4 ): 1195 – 1200 , jan 2007 . ISSN 00278424 . doi: 10.1073/pnas.0601378104 . OpenUrl Abstract / FREE Full Text 107. ↵ Jean-Christophe Leloup and Albert Goldbeter . Temperature compensation of circadian rhythms: control of the period in a model for circadian oscillations of the per protein in drosophila . Chronobiology international , 14 ( 5 ): 511 – 520 , 1997 . OpenUrl PubMed Web of Science 108. ↵ Roy M Daniel , Mark Dines , and Helen H Petach . The denaturation and degradation of stable enzymes at high temperatures . Biochemical journal , 317 ( 1 ): 1 – 11 , 1996 . OpenUrl Abstract / FREE Full Text 109. ↵ Benjamin D Knapp and Kerwyn Casey Huang . The effects of temperature on cellular physiology . Annual Review of Biophysics , 51 ( 1 ): 499 – 526 , 2022 . OpenUrl CrossRef PubMed 110. Peter JH Sharpe and Don W DeMichele . Reaction kinetics of poikilotherm development . Journal of theoretical biology , 64 ( 4 ): 649 – 670 , 1977 . OpenUrl CrossRef PubMed Web of Science 111. ↵ Robert M Schoolfield , PJH Sharpe , and Charles E Magnuson . Non-linear regression of biological temperature-dependent rate models based on absolute reaction-rate theory . Journal of theoretical biology , 88 ( 4 ): 719 – 731 , 1981 . OpenUrl CrossRef PubMed Web of Science 112. ↵ Terence L Wagner , Hsin-I Wu , Peter JH Sharpe , Robert M Schoolfield , and Robert N Coulson . Modeling insect development rates: a literature review and application of a biophysical model . Annals of the Entomological Society of America , 77 ( 2 ): 208 – 220 , 1984 . OpenUrl CrossRef 113. ↵ Simen Jacobs , Federico Vazquez , Nikita Frolov , and Lendert Gelens . Gitlab repository . https://gitlab.kuleuven.be/gelenslab/publications/Jacobs_Beyond_Arrhenius.git , 2025 . 114. Y Guan , Z Li , S Wang , PM Barnes , X Liu , H Xu , M Jin , AP Liu , and Q Yang . A robust and tunable mitotic oscillator in artificial cells . eLife , 7 : e33549 , 2018 . OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted September 04, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Beyond Arrhenius: Nonlinear and negative temperature scaling of biological rates from multi-step mechanisms Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Beyond Arrhenius: Nonlinear and negative temperature scaling of biological rates from multi-step mechanisms Simen Jacobs , Federico Vazquez , Nikita Frolov , Lendert Gelens bioRxiv 2025.09.01.673554; doi: https://doi.org/10.1101/2025.09.01.673554 Share This Article: Copy Citation Tools Beyond Arrhenius: Nonlinear and negative temperature scaling of biological rates from multi-step mechanisms Simen Jacobs , Federico Vazquez , Nikita Frolov , Lendert Gelens bioRxiv 2025.09.01.673554; doi: https://doi.org/10.1101/2025.09.01.673554 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Systems Biology Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41913) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13372) Ecology (19889) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40365) Molecular Biology (17163) Neuroscience (88540) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15136) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9818) Zoology (2269)