Quantitative modelling of biological response dynamics reveals novel patterns in plant volatile signalling

23 Biological responses to environmental stimuli are inherently dynamic. Recent technological 24 advances enable detailed time-resolved measurements of such responses. However, a 25 unifying model for the quantitative characterisation of dynamic response curves is lacking, 26 thus limiting biological insights and comparisons. We developed an unbiased mathematical 27 modelling approach that allows for the quantitative characterisation of biological response 28 curves without a priori knowledge of underlying biochemical mechanisms. Using stress-29 induced plant volatile emissions, we show that the model performs well across a wide range 30 of datasets, including incomplete induction curves, curves with low sampling resolution and 31 complex multi-modal responses that result from overlapping stimuli. Using the model, we 32 uncover a range of previously unrecognised patterns in volatile induction curves, including i) 33 a strong light-independent impact of the time of day of wounding on the onset, duration and 34 shape of the volatile induction responses, ii) an accentuation of volatile-specific induction 35 curve shapes by herbivory-associated molecular patterns (HAMPs) and iii) independent 36 regulation of the strength and duration of volatile induction across genotypes. We conclude 37 that our model allows for the quantitative analysis of dynamic response curves and can 38 identify new biological response patterns. The presented approach will pave the way to 39 characterising dynamic organismal responses across the tree of life. 40 41 42 43 44 45 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 3

46 Inducible biological responses are transient in nature and change substantially over relatively 47 short periods of time (1–5). By consequence, real-time measurements are crucial to understand 48 physiological processes, precisely diagnose chronic and acute stress patterns and resolve 49 ecological p henomena (6–11). Although biological processes are inherently dynamic, 50 quantitatively describing and comparing temporal features of biological responses remains a 51 challenge. It is well established that static features such as mean concentration at a given 52 instance are useful for general interpretations of stimulus -response relationships, such as the 53 inverse relationship between viral load and effective immune response in humans and the 54 positive relationship between herbivore damage and volatile emissions in plants (12, 13). 55 However, static, nominal values, even if measured at different time points, ignore meaningful, 56 quantifiable features of the dynamics themselves. 57 Recent examples illustrate the potential of moving beyond measures such as concentration or 58 abundance to generate new biological insights. In plants, variation in the spatial distribution of 59 a toxic metabolite, either within a plant or between plants, can affect herbivore feeding 60 behaviour and thus fitness parameters such as growth, independently from the total 61 concentration of said toxin (14–16). In animals, dendritic signals are sensitive to temporal 62 patterns of synaptic activation; irrespective of the amount of signal, the temporal patterns of 63 signal perception shape spike outputs, which drive crucial functions such as sound localisation 64 (17, 18). Additionally, the interaction between vaccination status, viral exposure and immune 65 responses are all temporally linked, and the effects/response of each will be dependent on how 66 the others change over time (19). 67 From a mechanistic perspective, the dynamics of biological responses are determined by the 68 interplay of highly regulated signalling events (20–25) and modifications to one or more events 69 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 4 can have far-reaching consequences on response dynamics and metabolic outcomes (1, 17, 26–70 28). For diffusion -based processes such as generation of reactive oxygen species (ROS) or 71 single-cell signalling events like MAPK cascades, the underlying steps are few and well -72 characterised, and relatively simple, real-time mechanistic models have been developed (6, 29–73 32). Some models of temporal dynamics of metabolism with more complex biosynthetic 74 pathways have also been developed (33, 34). These models are valuable for elucidating 75 biochemical mechanisms and understanding the complexity of metabolic regulation. However, 76 they require substantial a priori biochemical knowledge and the ability to measure numerous 77 kinetic and regulatory parameters (steps), making their construction challenging even for well-78 characterised metabolic processes (35–37). For more complex and less -well characterised 79 pathways, such as the de novo production of specialised metabolites, complete time -resolved 80 mechanistic models are not currently feasible. Many underlying processes, including enzyme 81 kinetics, precursor and product concentrations, transport delays, degradation pathways and 82 other regulatory processes remain unknown (38–40). Because the biochemistry is not fully 83 defined, the number and nature of parameters required to describe such processes are 84 effectively arbitrary, leading to under -constrained models in which multiple parameter 85 combinations can reproduce the same behaviour and limit reliable fitting and meaningful 86 interpretation (37). Moreover, existing models are often system -specific and must be adapted 87 depending on the metabolites or pathways studied. Consequently, there is no universally 88 standardised framework or parameterisation that can be easily implemented across systems, 89 which limits comparability, broader applicability and ultimately accessibility of dynamic 90 metabolite analyses. 91 Mechanistically informed models may not always be necessary to study biological responses. 92 V olatile chemical signals emitted by organisms shape important interactions such as those 93 between insects and their hosts (41, 42). In this case, the mechanism of volatile formation 94 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 5 matters less than when this chemical information is released, how long it persists and whether 95 amounts produced cross perceptual thresholds (10, 43–45). These characteristics are also 96 important for within -organism biochemical processes. Consider protein structural dynamics 97 (i.e., moving between different conformational states in time), for which similar dynamic 98 features would determine interactions with other proteins and metabolites, and thus shape the 99 fate of biochemical processes (46, 47). Thus, extracting time-resolved features of response 100 dynamics that are intuitive, broadly comparable and important for driving biological functions 101 becomes important. Both statistical and theoretical methods for modelling dynamic biological 102 responses have been proposed (8, 13, 19). However, a universal model that enables 103 quantitatively rigorous and biologically informative comparison of dynamic response features, 104 either across stimuli or across systems, has yet to be developed. 105 106 Motivated by the lack of such a framework to quantitatively compare dynamic biological 107 responses, we set out to develop a method that allows for meaningful, robust and consistent 108 comparisons of dynamic responses without explicit knowledge of the underlying 109 physiochemical and biochemical processes. We then generated a series of plant stress 110 response datasets to test our method and evaluate whether it can generate new insights into 111 biological responses. We considered plant stress responses, including induced volatiles, to be 112 an optimal test case, as plants have evolved a highly complex array of responses to 113 environmental stimuli, that cover a broad range of physiological and ecological functions 114 (48), and show starkly distinct kinetics (42). Induced volatiles in particular can be measured 115 non-destructively, enabling diverse processes to be observed in real time (7). Finally, plant 116 metabolites serve as important chemical information for the wider environment, and we thus 117 assumed that the generated insights would be informative on a biological and ecological level 118 (10, 45). The results demonstrate that our method is useful to compare biological responses. 119 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 6 Most importantly, it allowed us to uncover differences in plant stress response dynamics that 120 would not be captured with traditional approaches, thus resulting in novel insights into the 121 specificity and variation of biological responses. As the method is not specific to either 122 response type or timeframe, it will be applicable over a broad range of dynamic processes 123 across biological systems, thus unlocking novel comparative and integrative approaches 124 across the tree of life. 125 126

127 Model development 128 To build a phenomenological model that captures dynamic responses in an unbiased manner, 129 we modelled responses as a distribution of waiting times. We asked the question: How long 130 does it take for a response unit (e.g. a volatile molecule) to appear following a stress event? In 131 doing so, we interpret each measured unit as a probabilistic occurrence with a certain delay, 132 and thus model the response curve as a probability density function. 133 This framing is advantageous for three reasons: first, the model yields a relatively small number 134 of parameters that are both identifiable ( i.e. are intuitive by eye) and interpretable (i.e. have 135 biological meaning). This minimises overfitting and, as a result, makes the model robust to 136 sample-to-sample variation. Secondly, it enables comparability across responses which involve 137 very different mechanisms and timeframes. Thirdly, the same framework can be applied to any 138 induced phenotype that can be measured quantitatively over time, such as gene expression, 139 hormone accumulation, enzymatic conversions, volatile emissions and resistance induction, 140 even if they show starkly different dynamics and biosynthetic processes, as they all typically 141 exhibit the same behaviour: induction delay following stimulation and a rapid increase in 142 response, which, after reaching a peak, gradually declines (1, 2, 28, 49–57). This is 143 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 7 mathematically intuitive as these responses are built from a sequential chain (or combinatory 144 network) of activation, synthesis, transport and eventually measurement/detection. Each step 145 adds stochastic waiting time with their sum naturally producing a skewed, unimodal 146 distribution (58). Modelling this kind of process usually involves treating each biochemical 147 step as a stochastic waiting time drawn from a defined probability distribution and then 148 combining these into a final waiting -time distribution which captures the total observed 149 response dynamics. As such, we used a reparametrised gamma distribution , which retains its 150 form when steps are combined and can thus describe both early and later phases of the response 151 in a consistent way (see materials and methods for full model details): 152 𝑅(𝑡) = 𝑅peak ⋅ ' 𝑡 − 𝑡onset 𝑡peak − 𝑡onset ) !"mean#"onset "mean#"peak #$% ⋅ exp '−-𝑡 − 𝑡peak. 𝑡mean − 𝑡peak ) 153 where: 154 • 𝑅&'() is the maximum measured response, 155 • 𝑡*+,'" is the onset delay, 156 • 𝑡&'() is the time of peak response, 157 • 𝑡-'(+ is the mean response time. 158 159 The fitting parameters can then be used to calculate additional characteristic parameters (Fig 160 1A). 161 Total response (𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙) – The integral is often difficult to measure experimentally ; due to 162 time limitations in the measurement processes, entire curves are not resolved and thus 163 incomplete integrations are used (59). When fit onto the exper imental data, our model can be 164 used to calculate the theoretical integral of the response curve even if it is incomplete: 165 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 8 Integral = 𝑅!"#$ ,𝑡%"#& − 𝑡!"#$/  ()*!"#$+*%&'"( *)"#&+*!"#$ 𝑒 *!"#$+*%&'"( *)"#&+*!"#$  ,𝑡!"#$ − 𝑡,&-"* / +*!"#$+*%&'"( *)"#&+*!"#$  Γ 31 + 𝑡!"#$ − 𝑡,&-"* 𝑡%"#& − 𝑡!"#$ 6 166 167 where Γ() is the Gamma function, a continuous generalisation of the factorial (46, 47). 168 169 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 – For response dynamics, response length is typically defined as ‘broadness’, for 170 example as full width at half max (FWHM) (6). The FWHM of the Gamma function has no 171 analytical solution and must be determined numerically (62). Additionally, FWHM, does not 172 explicitly capture the start of the response. For these reasons 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 is used as a simple 173 alternative, achieving a similar descriptive value and used as a normali sation value in the 174 ‘shape’ feature below. 175 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 = 𝑡-'(+ − 𝑡*+,'" 176 𝑆ℎ𝑎𝑝𝑒 – The ‘shape’ of a response curve is an abstract feature. Since our model is intended to 177 be used over variable time scales, we developed a time-normalised (by 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛) metric that 178 allows for comparison of curve shapes independently from response length. Further, we aimed 179 to understand the symmetry of the curve in response to different stimuli or across plant 180 genotypes/species to broadly assess temporal nuances. 𝑆ℎ𝑎𝑝𝑒 presents a new phenotypic axis 181 which may link to important biochemical or physiochemical properties that regulate induced 182 responses. 𝑆ℎ𝑎𝑝𝑒 is a normalised value , and w hen responses follow gamma distributions it 183 falls between zero and one; a 𝑆ℎ𝑎𝑝𝑒 of zero indicates a perfectly symmetrical curve (𝑡-'(+ =184 𝑡&'() ) and a shape of one represents a completely right-skewed curve (𝑡&'() = 𝑡*+,'" ): 185 𝑆ℎ𝑎𝑝𝑒 = 𝑡-'(+ − 𝑡&'() 𝑡-'(+ − 𝑡*+,'" 186 187 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 9 Taken together, the fit and derived parameters give a framework for compa ring curves, 188 particularly those that appear distinct by eye but are difficult to distinguish quantitatively (Fig 189 1B). 190 191 192 Uncovering novel patterns in plant volatile responses 193 To explore the usefulness of the model , we measured inducible volatile emission s in maize 194 (Zea mays) in real -time by PTR -ToF-MS following different stress treatments. We focused 195 principally on the homoterpene 4,8-dimethylnona-1,3,7-triene (DMNT) as a highly inducible, 196 ecologically relevant plant volatile (63–65). 197 First, we tested whether the amount of mechanical leaf damage influences DMNT induction. 198 Earlier work showed that damage intensity strongly correlates with terpene emissions (12). 199 Indeed, higher damage increased total DMNT emissions DMNT (Fig. 2A), and our model 200 Figure 1. Theoretical examination of important features of biological response curves. A) Parameters of the model include 𝑡!"#$%: onset delay, 𝑡&$'(: time of peak response, 𝑡)$'": mean response time, 𝑅&$'(: the maximum measured response . These parameters are further reparametrised into 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛: metric of the length of response and 𝑆ℎ𝑎𝑝𝑒: our symmetry feature , which for gamma -like distributions falls between zero and one. B) In silico data highlighting the utility of model parameters, namely that discrete temporal response patterns can emerge, despite identical integrals. .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 10 faithfully captured this pattern through 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙 (Fig 2B). No clear patterns were detected for 201 the other parameters, apart from differences in 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛, which was marginally longer for 202 intermediate amounts of wounding (Fig. 2D). 203 Next, we tested the influence of the circadian clo ck on wound-induced DMNT by wounding 204 plants at different times of day under continuous light. The circadian clock is known to play a 205 role in regulating how plants respond to stress, including transcriptional regulation of defence 206 genes (66, 67). Time of day had no impact on the overall amount of DMNT emission (Fig 2G). 207 However, 𝑡*+,'" , 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 and 𝑆ℎ𝑎𝑝𝑒 all varied significantly ; wounding in the evening 208 resulted in a significantly more rapid, shorter DMNT burst – a novel pattern that has not been 209 reported before (Fig. 2F-J). 210 We also tested the influence of insect oral secretions (OS), which contain both elicitors and 211 effectors which can increase and decrease responses compared to wounding alone, respectively 212 (68). We confirmed that OS enhance the total emission of DMNT (Fig 2K -L). Interestingly, 213 OS treatment also led to a more rapid and prolonged response, resulting in a significantly 214 different curve shape (Fig 2M-O). 215 Leaf size and developmental stage can influence total volatile emissions (49), but whether these 216 parameters also influence response curves in other ways is unknown . We thus reanalysed the 217 dataset from our earlier work with our model. We found that leaf 3 (largest size and 218 intermediate age) emits volatiles for the longest period (Fig 2S). Interestingly, leaf 2 (olde r 219 leaf), despite producing more volatiles than leaf 4 (youngest leaf), did not produce them for a 220 longer period (Fig 2R). 𝑆ℎ𝑎𝑝𝑒 followed a clear developmental gradient, whereby older leaves 221 had more asymmetrical DMNT emission dynamics compared to younger leaves (Fig 2T). 222 Different genotypes vary strongly in the amount and type of volatiles they emit (69), but 223 whether response curves are also different is unknown. We determined response curves in three 224 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 11 maize inbred lines that differ in their capacity to produce DMNT (Fig 2 U). As expected , 225 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙 varied for all genotypes (Fig 2 V). Interestingly, we found that high er emitting 226 genotypes also had substantially earlier 𝑡*+,'" (Fig 2 W), however higher emissions do not 227 necessarily equate to earlier onsets (Fig 2A -C). Additionally, the t wo highest emitting 228 genotypes, CML287 and NC300, had similar 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 (Fig 2X) and CML287 had the most 229 symmetrical dynamics . Thus, different genotypes show starkly different response patterns, 230 some of which are independent of total volatile quantity. Taken together, these experiments 231 illustrate the power of our approach to uncover novel, genetically determined response patterns 232 to environmental stimuli. 233 234 Unravelling differences between biochemically distinct compound classes 235 To explore the usefulness of our model to compare different types of response curves, we 236 characterised volatiles from different biosynthetic pathways, including terpenes, indole and 237 green leaf volatiles (1). In response to wounding, the different volatiles exhibited significantly 238 different response cur ves, as described before, which was visible in differences in 𝑡*+,'" , 239 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 and 𝑆ℎ𝑎𝑝𝑒 (Fig. 3A -E; Fig S1 ). Interestingly, the differences became more 240 pronounced with the application of OS (Fig. 3F-J). This pattern was most clearly visible in the 241 𝑆ℎ𝑎𝑝𝑒 parameter: while i n the absence of OS, all compounds except monoterpenes had a 242 similar 𝑠ℎ𝑎𝑝𝑒, all but DMNT and TMTT developed their own unique 𝑠ℎ𝑎𝑝𝑒 when induced by 243 OS (Fig 3E and J), illustrating that OS components act dif ferently on different volatile 244 pathways. 245 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 12 246 Figure 2. Application of the model to compare volatile emission dynamics in response to a range of stimuli. Each horizontal row of panels represents a unique experiment. A-E) responses to variable wounding intensities, F-J) responses to the same intensity of damage at different times of day (Note: 𝑡!"#$% is standardised based on time of damage), K-O) responses compared between wounded plants and plants treated with wounding and Spodoptera exigua oral secretions (OS), P-T) responses to damage in leaves of different developmental stages (Data from Waterman et al., 2025), U-Y) wounding responses compared between genotypes with highly variable volatile emission capacity. For the first column of panels, curves depict emission data, where the solid line represents mean of fitted emission across biological replicates and the translucent ribbon represents the baseline-subtracted raw emission data ± SE. For the remaining columns, solid points represent mean values across biological replicates (translucent po ints). Error bars represent SE. Within each panel, different letters indicate significant differences between groups as determined by multiple comparisons tests following significant (p < 0.05) one-way ANOV A or Kruskal-Wallis tests. n = 4-6. .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 13 247 248 249 250 251 Figure 3. Impacts of herbivore-specific stimuli on volatile emission dynamics. Emission and model parameters for wounded plants (A-E) and wounded plants treated with oral secretions (OS; F-J). For the first column of panels (A and F), curves depict emission data, where the solid line represents mean of fitted emission across biological replicates and the translucent ribbon represents the baseline-subtracted raw emission data ± SE. For the remaining columns, solid points represent mean across biological replicates (translucent points). Error bars represen t SE. Within each panel, different letters indicate significant differences between groups as determined by multiple comparisons tests following significant (p < 0.05) one-way ANOV A. n = 4-5. Abbreviations: DMNT= 4,8-dimethylnona-1,3,7-triene, MNT = monoterpenes, SQT = sesquiterpenes, TMTT = 4,8,12-trimethyltrideca-1,3,7,11-tetraene. Supplemental Figure 1. Dynamics of green leaf volatile (GLV) emissions. For A) curves depict emission data, where the solid line represents mean of fitted emission across biological replicates and the translucent ribbon represents the baseline-subtracted raw emission data ± SE. For the remaining columns, solid points represent mean across biological replicates (translucent points). Error bars represent SE. Within each panel, different letters indicate significant differences between groups as determined by multiple comparisons tests following significant (p < 0.05) one- way ANOV A. n = 3. Abbreviations: H-al = hexenal, H-ol = hexenol, HAC = hexenyl acetate. .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 14 To explore this phenomenon more deeply and test our model on a less temporally resolved 252 dataset, we modelled the gene expression dynamics of ZmCYP92C5, ZmIGL, ZmTPS2 and 253 ZmTPS10 (Fig S2), which are coding for the rate limiting enzymes of the biosynthesis pathways 254 of DMNT/TMTT, indole, monoterpenes and sesquiterpenes (49). Trends in 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 of these 255 genes matched emission of corresponding volatiles (Fig S 2C), confirming that volatile 256 emission is regulated by biosynthe tic limitations (Fig S 2; 12, 56). No differences in gene 257 expression were found for 𝑡*+,'" and 𝑠ℎ𝑎𝑝𝑒 (Fig S2B and D) , possibly due to a lack of 258 resolution of expression dynamics (Figs S3 and S4). 259 260 261 262 Exploring complex damage and response patterns 263 To test the ability of our model to characterize more complex induction patterns, we measured 264 volatile emissions following multiple wounding events (Fig 4A). The current assumption from 265 the literature is that subsequent wounding should lead to stronger defence responses beyond 266 Supplemental Figure 2. Volatile biosynthesis gene expression dynamics . For A), curves depict baseline-subtracted raw expression data, where the solid line represents mean of fitted expression across biological replicates and the translucent ribbon represents the baseline-subtracted raw expression data ± SE. For the remaining columns, solid points represent mean across biological replicates (translucent points). Error bars represent SE. Within each panel, different letters indicate significant differences between groups as determined by multiple comparisons tests following significant (p < 0.05) one-way ANOV A. n = 4-5. Abbreviations: CYP92C5 = dimethylnonatriene/trimethyltetradecatetraene synthase, IGL = indole-3-glycerol phosphate lyase, TPS2 = terpene synthase 2, TPS10 = terpene synthase 10 . .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 15 cumulative effects, as the earlier wounding will prime plants for subsequent responses. 267 However, such effects have been hard to isolate for wounding events that follow each other 268 closely in time due to overlapping response curves. Thus, we modelled the responses to three 269 wounding events as the sum of three separate curves, all of which being effectively incomplete 270 or not fully resolved in some way. By decomposing the overall emission into separate fitted 271 curves, we were able to quantify the dynamic contribution of each individual wounding event, 272 even when the peaks overlapped. This revealed a clear priming effect where the 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑙 of 273 the second and third responses were larger than that of the initial peak, demonstrating that prior 274 damage enhanced subsequent emissions beyond additive effects (Fig 4B) . At this particular 275 time interval, multiple damage events did not significantly modify other fit parameters (Fig 276 4C-E). The functionality of this model on more complex damage patterns highlights its utility 277 to explore response dynamics under complex stress regimes. 278 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 16 279 280 Herbivores do not make single wounds when they feed, but rather take many bites over time . 281 To test the application of the model under realistic conditions, we fit it to response curves to 282 feeding by a single Spodoptera exigua larva for 30 min. Natural herbivory patterns generated 283 responses that could easily be fitted for all groups of volatiles (Fig 4F). 𝑡*+,'" showed some 284 similarities to wounding response patterns, namely that sesquiterpenes and TMTT take longer 285 to begin emitting than other compounds (Fig 4G). Interestingly, 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 values more closely 286 matched those from simple wounding than wounding + OS (Fig 4H). 𝑆ℎ𝑎𝑝𝑒 values showed 287 Figure 4. Fitting responses to complex stimuli patterns. A) Fitted emission for each curve plotted over total fitted emission and baseline - subtracted raw emission B-E) Model parameters for each peak. F) Curves depict emission data and G-I depict model parameters from Spodoptera exigua-infested plants. For A and F, the translucent ribbons represent the baseline-subtracted raw emission data ± SE from the respective, colour- coded curve. For B-E and G-I, solid points represent mean across biological replicates (translucent points). Error bars represent SE. Within each panel, different letters indicate significant differences between groups as determined by multiple comparisons tests followin g significant one- way ANOV A. For A-E, n = 6 and for F -I, n = 6-8. Abbreviations: DMNT= 4,8-dimethylnona-1,3,7-triene, MNT = monoterpenes, SQT = sesquiterpenes, TMTT = 4,8,12-trimethyltrideca-1,3,7,11-tetraene. .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 17 stark differences compared to wounding treatments, suggesting that damage patterns may play 288 an important role in determining response curves (Fig 4I). The ability to effectively fit 289 parameters to responses induced by real herbivore feeding indicates that our model is functional 290 even when exact damage patterns and timing are not known, and thus it can be used widely. 291 292 Modelling incomplete curves 293 Understanding the dynamics of defence responses requires some level of continuous 294 monitoring. However, how the temporal resolution of measurements impacts the capacity to 295 understand dynamic patterns is not well known. We tested the sensitivity of our model to 296 measurement duration and sampling resolution. Firstly, in order to simulate a shorter 297 measurement window, we removed datapoints from the end of the curve (in time) and 298 compared fit parameters to those from the complete dataset. For 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑙 and 𝑑𝑢𝑟𝑎𝑡𝑖𝑜𝑛, 299 there was a consistent point around 𝑡-'(+ where removal of data dramatically reduce the 300 accuracy of model parameters (Fig S3B, G and L). 𝑡*+, '" was more sensitive to data removal 301 for indole and TMTT than it was for DMNT (Fig S3C, H and M). Similarly for 𝑠ℎ𝑎𝑝𝑒, which 302 is 𝑡*+,'" -dependent), indole and TMTT fits broke down relative to the complete dataset with 303 less data removal than DMNT (Fig S3E , J and O). Interestingly, for DMNT and indole, even 304 with removal of ca. 10-15 hr of curve the fitted parameters remained almost entirely stable, if 305 not identical, suggesting a substantial degree of flexibility and potential to understand the 306 dynamics of only partially resolved peaks. 307 Second, we tested the resolution tolerance of the model by periodically removing datapoints, 308 simulating a lower resolution measurement (Fig S4 A, F and K). With periodic removal of up 309 to 85% of data points, the emission 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑙 was able to be recovered with nearly the same 310 accuracy and precision for all compounds, suggesting that a coarse ‘outline’ of temporal 311 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 18 response patterns can be sufficient to accurately estimate total emissions over time (Fig S4B, 312 G and L). For indole, removal of 70% data in periodic intervals resulted in breakdown of 313 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑙 as the overall curve structure was nearly entirely lost (Fig S4B). Indole has the fastest 314 emission dynamics, and as such, indole has the lowest effective resolution in comparison to 315 DMNT and TMTT, which both have slower kinetics and thus more data points between 𝑡*+,'" 316 and 𝑡-'(+ . Interestingly, for other fit parameters (𝑡*+,'" , 𝑑𝑢𝑟𝑎𝑡𝑖𝑜𝑛 and 𝑠ℎ𝑎𝑝𝑒), accuracy and 317 precision began to consistently breakdown at 70% removal, suggesting there is a limit to 318 measurement frequency for certain dynamic features (Fig S4). Overall, indole appeared to be 319 the most susceptible to reduced resolution, again, likely due to indole’s faster dynamics and 320 lower effective resolution. Importantly, these results highlight that response kinetics are a 321 necessary consideration to determine the use of truly ‘continuous’ measurements; they may not 322 always be required or even provide additional dynamic information. 323 324 325 326 327 328 329 330 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 19 331 Supplemental Figure 3. Fitting incomplete response curves. Emission and fit parameters for A-E) indole, F-J) DMNT and K-O) TMTT. For the first column of figures, curves depict emission data, where the each solid line represents mean of fitted emission across biological replicates and the translucent ribbon represents the raw baseline-subtracted emission data ± SE. For the remaining columns, solid points represent mean across biological replicates (translucent points). Error bars represent SE. Colour gradient indicates the number of hours removed from the back end of the curve. n = 5-6. Abbreviations: DMNT= 4,8-dimethylnona-1,3,7-triene, TMTT = 4,8,12-trimethyltrideca-1,3,7,11-tetraene. .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 20 332 Supplemental Figure 4. Sampling resolution impacts fits. Data were removed systematically in periodic intervals across the entire curve. Emission and fit parameters for A-E) indole, F-J) DMNT and K-O) TMTT. For the first row of data, curves depict emission data, where the solid line represents mean of fitted emission across biological replicates and the translucent ribbon represents the baseline-subtracted raw emission data ± SE. For the remaining columns, solid points represent mean across biological replicates (translucent points). Error bars represent SE. n = 5-6. Abbreviations: DMNT= 4,8-dimethylnona-1,3,7-triene, TMTT = 4,8,12-trimethyltrideca-1,3,7,11-tetraene. .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 21

333 Dynamic responses to environmental stimuli are a universal property of life. Yet, unifying 334 approaches to characterize such responses are lacking. Here, we demonstrate full categorisation 335 and quantification of response features with a relatively simple model and demonstrate its 336 usefulness in revealing new biological patterns. The model improves on traditional response 337 quantification by providing true integrals and maximum response values as well as additional 338 response parameters that are not easily extractable otherwise. 339 Across biology, we have major gaps in our understanding of patterns and differences in 340 dynamic biological responses. Induced volatile emissions are an illustrative example in this 341 context. In a series of experiments, we find that our model reveals hitherto undocumented 342 differences in volatile induction along circadian, developmental and genetic axes, and in 343 response to wounding and insect -specific molecular patterns. Among the most noteworthy 344 discoveries is that the time of day of wounding has a strong impact on the onset, duration and 345 shape of the volatile induction curve independently of the total response. This provides a 346 window into the tight regulation of the speed and duration of induced volatile emissions by the 347 circadian clock (10, 71), providing new phenotypes that link clock regulation to environmental 348 interactions. Equally noteworthy is the finding that insect oral secretions (OS) enhance 349 differences in response curve shapes between different volatile classes, thus resulting in unique 350 temporal volatile fingerprints . This newly discovered pattern may explain OS -specific 351 responses that affect interactions with herbivores and herbivore natural enemies beyond 352 differences in overall volatile quantities or timing (72). Finally, our model uncovers that, 353 contrary to current expectations, the quantity and duration of induced volatile emissions can be 354 regulated independently of each other across plant genotypes. This finding paves the way to 355 identify regulators of volatile induction duration through genetic approaches, potentially 356 leading to the identification of new “late signalling components” involved in sustaining defence 357 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 22 responses. Further, these novel dynamic features might ultimately be leveraged to non -358 destructively distinguish and identify stimuli, for example in an agricultural context , where 359 many individual induced responses may be shared between multiple stressors and total 360 response is thus insufficiently informative (6, 45, 50). 361 Recent advances in continuous , real-time monitoring have greatly improved our ability to 362 measure biological response dynamics (5, 7, 10, 50), yet analyses of such data remain limited. 363 Our modelling approach addresses this gap by providing a standardised framework for 364 quantifying response dynamics in a way that enables biologically and ecologically meaningful 365 comparisons across systems. In this study, we use volatile emissions and gene expression from 366 plants as representative measures of broad biological responses. Thanks to the robustness and 367 generality of the model, it will be straightforward to apply it across biological systems, from 368 energy dynamics in microbes to electron transport in plants to neurological and immune 369 responses in animals (13, 17, 73, 74). Thereby, the framework opens opportunities for 370 ambitious endeavours, such as tree-of-life scale meta-analyses, to uncover shared and divergent 371 principles of organismal responses to environmental stimuli and to reveal new fundamental 372 biological phenomena. Furthermore, coupling model parameters with machine -learning 373 algorithms or constrained optimi sation solvers could enable the reconstruction of stimulus 374 patterns from observed dynamics, or the prediction of dynamic responses from environmental 375 inputs. We propose this framework as a universal standard for intuitive, interpretable and 376 biologically grounded quantification of dynamic responses, to advance understanding of 377 fundamental processes and to guide innovation across medicine, agriculture and beyond. 378 To ensure that the model can be used widely, we generated a set of freely accessible resources 379 that can be implemented in Python, R and Excel formats (Appendix I), and applied broadly to 380 time-resolved response data to yield all relevant readouts for further statistical analyses. 381 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 23

382 Plants and insects 383 Plant growth conditions were identical to those used in Waterman et al . (2025). In brief, V2 -384 stage (two developed leaves, one expanding leaf and one emerging leaf) maize (Zea mays) 385 seedlings were used throughout the study. At this stage maize plants are particularly susceptible 386 to agricultural pests such as Spodoptera spp. (75). Plants were grown in commercial potting 387 soil (Selmaterra, BiglerSamen, Switzerland) in 180 mL pots under greenhouse conditions and 388 supplemented with artificial lights (ca. 300 µmol m−2 s−1). The greenhouse was kept at 22 ± 389 2°C, 40%–60% relative humidity, with a 14 h: 10 h, light: dark cycle. 390 Spodoptera exigua larvae (Frontier Agricultural Sciences, USA) were reared from eggs on an 391 artificial diet (76). At least 24 hr prior to experimental treatments larvae were fed B73 maize 392 leaves. Oral secretions (OS) were collected by probing the mouths of larvae with a pipette tip 393 and stored at -20ºC until use. 394 395 Experimental Treatments 396 All wounds were inflicted using haemostat forceps, an established method of simulating 397 herbivory (49, 68). The maize inbred line, B73, was used unless otherwise stated. The standard 398 damage intensity was 120 mm2 and time of initial damage was 11:00 unless otherwise stated. 399 To understand how different damage regimes impact induced response dynamics we conducted 400 several experiments spanning a range of treatments: 401 Variable damage intensity: We damaged plants at 20, 40, 80, 120, 160 and 320 mm 2 on the 402 third-oldest leaf (leaf 3). Damage patterns always encompassed the central segment of each 403 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 24 leaf, and any damage above 40 mm2 was evenly distributed across the base, middle and tip of 404 leaf 3 (49). 405 Time of day: We damaged plants at 11:00, 11:45, 12:30, 13:15, 14:45, 15:45 and 23:00. 406 Oral secretions (OS): Damaged plants were either treated with milliQ water or 50% OS from 407 late-instar Spodoptera exigua larvae. 408 Leaf developmental stage: At the V2 stage maize has four leaves: two are fully developed (leaf 409 1 and leaf 2), one is actively expanding (leaf 3) and one is emerging (leaf 4). We used the raw 410 emission data from Waterman et al. (2025) to measure emission dynamics in each leaf. The 411 damage treatments in the previous study were 60mm2 total, and damage was inflicted in three 412 bouts of 20 mm2 damage in ca. 30 min intervals 413 Genotype: To explore the potential variation in response dynamics across maize genotypes we 414 damaged three additional maize inbred lines, with three distinct volatile emission capacities: 415 MO17 (low emission), NC300 (intermediate-high emission) and CML287 (high emission). 416 Real herbivory: 3rd instar Spodoptera exigua were placed on leaf 3 and were left to feed for 30 417 minutes. 418 419 V olatile sampling 420 V olatiles were collected as in Waterman et al. (2025). Briefly, entire seedlings were placed in 421 transparent glass chambers (Ø×H 12 × 45 cm) that were sealed other than a clean airflow inlet 422 and an outlet. Clean air was supplied at a flow‐rate of 0.8 L min −1. V olatiles were measured 423 with a high‐throughput platform comprising of a proton transfer reaction time‐of‐flight mass 424 spectrometer (PTR‐ToF‐MS; Tofwerk, Switzerland) and a custom‐made automated headspace 425 sampling system (Gonin et al. 2018). The outlet of the chamber was accessible to the 426 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 25 autosampler/PTR‐ToF‐MS system, which drew air at 0.1 L min −1. Between samples, a zero‐427 gas measurement was performed for 3 s to flush the system. At each time point (as indicated 428 by the x‐axis of the respective figure), volatiles were continuously measured for 25 –30 s and 429 averaged to a single mean per sample. Complete mass spectra (0 –500 m/z) were recorded in 430 positive mode at ca. 10 Hz. The PTR was operated at 100°C and an E/N of approximately 120 431 Td. The volatile data extraction and processing were conducted using Tofware software 432 package v3.2.2 (Tofwerk, Switzerland). Protonated compounds were identified based on their 433 molecular weight + 1. During volatile collection, LED lights (DYNA, heliospectra, Sweden) 434 were placed ca. 80 cm above the glass chambers and provided light at 300 μmol m−2 s−1. 435 Identical light:dark cycle timing as in the greenhouse for plant growth was used. 436 437 Gene expression 438 Leaf 3 tissue was harvested 0.5, 1, 1.5, 2, 3 and 5 hr after damage (see above for damage 439 treatment) and flash frozen on liquid nitrogen. Total RNA extraction and purification, genomic 440 DNA removal, cDNA synthesis and quantification of gene expression were conducted 441 identically to Waterman et al. 2024. Quantitative reverse transcription polymerase chain 442 reaction (qRT‐PCR) was performed using ORA SEE qPCR Mix (highQu GmbH, Germany) on 443 an Applied Biosystems QuantStudio 5 Real‐Time PCR system. The normalized expression 444 (NE) values were calculated as in Waterman et al. (2024) using ubiquitin (UBI1) as the 445

gene. Gene identifiers and primer sequences are listed in Table S1. 446 447 Modelling emission dynamics 448 For all experiments, volatile emissions were normalised by the biomass of leaf 3 (damaged 449 leaf). This is because we have previously shown that the size of the damaged leaf is limiting 450 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 26 for volatile emissions (49). Additionally, values were normalised to the maximum emission 451 observed in each experiment, yielding a range of positive values < 1. 452 453 Choosing the right distribution 454 To model the waiting time probability distribution, we used a Gamma distribution: 455 𝑓(𝑡; 𝛼, 𝛽) = 𝑡.#$𝑒#"/0 𝛽.Γ(𝛼) , 𝑡 ≥ 0 456 Where: 457 • 𝛼 is the number of steps, 458 • 𝛽 is the rate parameter (for each step), 459 • Γ(𝛼) is the Gamma function evaluated at 𝛼 460 461 This distribution arises from the sum of several exponential waiting times, making it suited to 462 processes where multiple sequential steps precede a single observable outcome (77). Unlike 463 the log-normal or inverse Gaussian, the Gamma distribution is closed under addition, meaning 464 that the sum of sequential Gamma is itself a Gamma distribution (78). This allows modelling 465 of both early and late phase responses, eve n if they are both part of the same chain of events. 466 In the formal definition, 𝛼 represents the effective number of steps and 𝛽 an effective average 467 delay, though we do not assume that these biological processes are truly memoryless and 468 independent as required by the strict derivation. Instead, we treat the Gamma distribution as a 469 phenomenological approximation of the process, which fits the shapes of the responses we (and 470 others) observe well (1, 2, 28, 49–57). With this in mind, we make a number of meaningful 471 reparameterisations to the standard Gamma distribution, as follows. 472 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 27 We begin with the standard Gamma distribution and introduce a shifted time variable 𝑡(12 to 473 account for a delay in the onset of the curve: 474 𝑡(12 = 𝑡 − 𝑡*+,'" . 475 The Gamma probability density function is then given by: 476 𝑓(𝑡; 𝛼, 𝛽, 𝑡*+,'" ) = 𝑡(12 .#$𝑒#"*+,/0 𝛽.Γ(𝛼) , 𝑡 ≥ 𝑡*+,'" . 477 478 Our next step is to relate the variables to observable quantities, and so we want to express 𝛼 479 and 𝛽 in terms of measurable features of the response curve. 480 For the Gamma distribution: 481 • 𝑀𝑜𝑑𝑒 = (𝛼 − 1)𝛽 = 𝑡&'() − 𝑡*+,'" 482 representing the time at which the response reaches its peak. 483 • 𝑀𝑒𝑎𝑛 = 𝛼𝛽 = 𝑡-'(+ − 𝑡*+,'" 484 representing the average time after onset that the response occurs. 485 In order to solve for 𝛼 and 𝛽 for these observable features we: 486 Subtract the first from the second: 487 𝑀𝑜𝑑𝑒 − 𝑀𝑒𝑎𝑛 = 𝛼𝛽 − (𝛼 − 1)𝛽 = (𝑡-'(+ − 𝑡*+,'" )− -𝑡&'() − 𝑡*+,'" . . 488 Simplify: 489 𝛽 = 𝑡-'(+ − 𝑡&'() . 490 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 28 491 Substitute into 𝛼𝛽 = 𝑡-'(+ − 𝑡*+,'" : 492 𝛼 = 𝑡-'(+ − 𝑡*+,'" 𝛽 = 𝑡-'(+ − 𝑡*+,'" 𝑡-'(+ − 𝑡&'() . 493 The Gamma probability density function is normalised to integrate to 1. 494 However, in experimental data the absolute amplitude (the response intensity , 𝑅&'() ) is an 495 observable quantity and so we define the response model: 496 𝑅(𝑡) = 𝑅&'() ⋅ 3(";.,0,"-./01) 38"20*3;.,0,"-./019 , 497 where 𝑅&'() is the observed response at 𝑡&'() . 498 This ensures that the curve shape follows the Gamma form, and that the maximum of 𝑅(𝑡) 499 equals 𝑅&'() . 500 The next step is to solve for 𝑓-𝑡&'() ; 𝛼, 𝛽, 𝑡*+,'" . 501 Substitute 𝑡 = 𝑡&'() into 𝑓(𝑡; 𝛼, 𝛽, 𝑡*+,'" ): 502 𝑓-𝑡peak; 𝛼, 𝛽, 𝑡onset. = -𝑡peak − 𝑡onset. .#$  exp L− 𝑡peak − 𝑡onset 𝛽 M 𝛽.Γ(𝛼) . 503 Replace 𝑡&'() − 𝑡*+,'" with (𝛼 − 1)𝛽 using the mode definition: 504 𝑓-𝑡peak; 𝛼, 𝛽, 𝑡onset. = [(𝛼 − 1)𝛽].#$ 𝑒𝑥𝑝 !#(.#$)0 0 % 𝛽. Γ(𝛼) . 505 Simplify the exponent: 506 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 29 exp '− (𝛼 − 1)𝛽 𝛽 ) = 𝑒#(.#$). 507 Final expression: 508 𝑓-𝑡&'() ; 𝛼, 𝛽, 𝑡*+,'" . = [(𝛼 − 1)𝛽].#$𝑒#(.#$) 𝛽.Γ(𝛼) . 509 We can therefore define: 510 𝑅&'() = 𝑓-𝑡&'() ; 𝛼, 𝛽, 𝑡*+,'" . = [(𝛼 − 1)𝛽].#$𝑒#(.#$) 𝛽.Γ(𝛼) . 511 Substituting back to obtain the final model: 512 Substitute 𝑓(𝑡)and 𝑓-𝑡&'() . into the ratio: 513 𝑅(𝑡) = 𝑅peak ⋅ (𝑡 − 𝑡onset).#$ exp L− 𝑡 − 𝑡onset 𝛽 M β;Γ(𝛼) -𝑡peak − 𝑡onset. .#$ exp L− 𝑡peak − 𝑡onset 𝛽 M 𝛽.Γ(𝛼) 514 515 Cancel the common terms 𝛽.Γ(𝛼): 516 𝑅(𝑡) = 𝑅peak ⋅ ' 𝑡 − 𝑡onset 𝑡peak − 𝑡onset ) .#$ ⋅ exp '− (𝑡 − 𝑡onset)− -𝑡peak − 𝑡onset. 𝛽 ) 517 Simplify the exponent: 518 𝑅(𝑡) = 𝑅&'() ⋅ ' 𝑡 − 𝑡*+,'" 𝑡&'() − 𝑡*+,'" ) .#$ ⋅ exp L− 𝑡 − 𝑡&'() 𝛽 M 519 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 30 Replace α and β with expressions in terms of -𝑡*+,'" , 𝑡&'() , 𝑡-'(+ . to give the final model: 520 𝑅(𝑡) = 𝑅peak ⋅ ' 𝑡 − 𝑡onset 𝑡peak − 𝑡onset ) !mean"!onset !mean"!peak "# ⋅ exp '− 𝑡 − 𝑡peak 𝑡mean − 𝑡peak ) 521 Where: 522 • 𝑅&'() is the maximum measured response, 523 • 𝑡*+,'" is the onset delay, 524 • 𝑡&'() is the time of peak response, 525 • 𝑡-'(+ is the mean response time. 526 527

subtraction 528 Since low levels of plant volatiles are constantly emitted e ven without wounding or other 529 measurable stimuli, a background curve was subtracted so that only induc ed emissions were 530 modelled (42). This background was estimated from undamaged plant s and either directly 531 subtracted or scaled to the pre -damage baseline to correct for instrumental drift or individual 532 plant batch variation. 533 534 Fitting 535 Each individual curve was fit using a two-stage optimisation approach. In order to avoid local 536 minima, we first used a differential evolution optimisation followed by a local refinement step 537 using SLSQP (79). Since the emission data were individually relatively noisy, we reconstructed 538 the dominant temporal trend using singular value decomposition (SVD) and incorporated this 539 as a weak prior to guide the fitting algorithm toward biologically realistic solutions (80, 81). 540 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 31 A number of constraints and bounds were enforced in order to guide the fitting of the data. 541 Temporal ordering was imposed such that 𝑡*+,'" < 𝑡&'() 0.1 and 𝑡-'(+ −543 𝑡&'() > 0.1. Bounds were also dynamically set for each curve; for instance, 𝑡&'() was 544 initialised at the time of maximum observed response and restricted to within ±20% of this 545 value, while 𝑅&'() was constrained between 80% and 120% of the observed maximum. 546 𝑡*+,'" and 𝑡-'(+ times were shifted according to experimental timing and constrained within 547 empirically informed ranges. All fitting, calculations and data manipulation were done using 548 SciPy (60), NumPy (82) and Pandas (83) in Python 3.12.4 (84). 549 550 Model testing 551 To explore the robustness of the model and fitting methods to the effects of resolution and 552 measurement time, a data set was taken and data were removed in order to simulate a lower 553 quality measurement. We chose the dataset of damage plants treated with insect OS, as a 554 representative and ecologically relevant subset. In a first test, data from the end of the 555 measurement were removed in a stepwise manner to simulate a shorter (incomplete) 556 experimental measurement period. Secondly, data were systematically down -sampled to 557 simulate a lower resolution measurement. After each step of the truncation or down sampling 558 the data was fit and the fit parameters compared. 559 560 Additional statistical analyses 561 All further statistical analyses were conducted in R version 4.4.2 (85). Differences in model 562 outputs between groups were determined using one-way ANOV A. Where necessary, data were 563 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 32 transformed to meet ANOV A assumptions. Where necessary, to obtain heteroscedasticity -564 consistent standard errors, White-adjusted ANOV As were used (86). Where data did not meet 565 the assumptions of ANOV A, even following transformation, differences were determined using 566 a Kruskal–Wallis test. Statistical test summaries are included in supplemental tables 1-4. 567 568

569 1. M. Erb, N. Veyrat, C. A. M. Robert, H. Xu, M. Frey, J. Ton, T. C. J. Turlings, Indole is 570 an essential herbivore-induced volatile priming signal in maize. Nature Communications 571 6, 6273 (2015). 572 2. Y.-Q. Gao, P. Jimenez-Sandoval, S. Tiwari, S. Stolz, J. Wang, G. Glauser, J. Santiago, E. 573 E. Farmer, Ricca’s factors as mobile proteinaceous effectors of electrical signaling. Cell 574 186, 1337-1351.e20 (2023). 575 3. H. P. Lu, L. Xun, X. S. Xie, Single-Molecule Enzymatic Dynamics. Science 282, 1877–576 1882 (1998). 577 4. M. Toyota, D. Spencer, S. Sawai-Toyota, W. Jiaqi, T. Zhang, A. J. Koo, G. A. Howe, S. 578 Gilroy, Glutamate triggers long-distance, calcium-based plant defense signaling. Science 579 361, 1112–1115 (2018). 580 5. L. Yang, V . Pathiranage, S. Zhou, X. Sun, H. Zhang, C. Lai, C. Gu, F. V . Subach, M. 581 Drobizhev, A. R. Walker, K. D. Piatkevich, Sensitive red fluorescent indicators for real-582 time visualization of potassium ion dynamics in vivo. PLOS Biology 23, e3002993 583 (2025). 584 6. M. C.-Y . Ang, J. M. Saju, T. K. Porter, S. Mohaideen, S. Sarangapani, D. T. Khong, S. 585 Wang, J. Cui, S. I. Loh, G. P. Singh, N.-H. Chua, M. S. Strano, R. Sarojam, Decoding 586 early stress signaling waves in living plants using nanosensor multiplexing. Nature 587 Communications 15, 2943 (2024). 588 7. P. Coatsworth, L. Gonzalez-Macia, A. S. P. Collins, T. Bozkurt, F. Güder, Continuous 589 monitoring of chemical signals in plants under stress. Nature Reviews Chemistry 7, 7–25 590 (2023). 591 8. R. Grote, M. Sharma, A. Ghirardo, J.-P. Schnitzler, A New Modeling Approach for 592 Estimating Abiotic and Biotic Stress-Induced de novo Emissions of Biogenic V olatile 593 Organic Compounds From Plants. Frontiers in Forests and Global Change Volume 2-594 2019 (2019). 595 9. L. Oyarte Galvez, C. Bisot, P. Bourrianne, R. Cargill, M. Klein, M. van Son, J. van 596 Krugten, V . Caldas, T. Clerc, K.-K. Lin, F. Kahane, S. van Staalduine, J. D. Stewart, V . 597 Terry, B. Turcu, S. van Otterdijk, A. Babu, M. Kamp, M. Seynen, B. Steenbeek, J. 598 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 33 Zomerdijk, E. Tutucci, M. Sheldrake, C. Godin, V . Kokkoris, H. A. Stone, E. T. Kiers, T. 599 S. Shimizu, A travelling-wave strategy for plant–fungal trade. Nature 639, 172–180 600 (2025). 601 10. J. M. Waterman, T. M. Cofer, L. Wang, G. Glauser, M. Erb, High-resolution kinetics of 602 herbivore-induced plant volatile transfer reveal clocked response patterns in neighboring 603 plants. eLife 12, RP89855 (2024). 604 11. J. Zhang, H. Yin, J. Zhang, G. Yang, J. Qin, L. He, Real-time mental stress detection 605 using multimodality expressions with a deep learning framework. Frontiers in 606 Neuroscience Volume 16-2022 (2022). 607 12. Ü. Niinemets, A. Kännaste, L. Copolovici, Quantitative patterns between plant volatile 608 emissions induced by biotic stresses and the degree of damage. Frontiers in Plant 609 Science Volume 4-2013 (2013). 610 13. J. Reyes-Silveyra, A. R. Mikler, Modeling immune response and its effect on infectious 611 disease outbreak dynamics. Theoretical Biology and Medical Modelling 13, 10 (2016). 612 14. V . S. Pan, E. Ghosh, P. J. Ode, W. C. Wetzel, K. J. Gilbert, I. S. Pearse, Large 613 Differences in Herbivore Performance Emerge From Simple Herbivore Behaviours and 614 Fine-Scale Spatial Heterogeneity in Phytochemistry. Ecology Letters 28, e70044 (2025). 615 15. I. S. Pearse, R. Paul, P. J. Ode, Variation in Plant Defense Suppresses Herbivore 616 Performance. Current Biology 28, 1981-1986.e2 (2018). 617 16. W. C. Wetzel, H. M. Kharouba, M. Robinson, M. Holyoak, R. Karban, Variability in 618 plant nutrients reduces insect herbivore performance. Nature 539, 425–427 (2016). 619 17. T. Branco, B. A. Clark, M. Häusser, Dendritic Discrimination of Temporal Input 620 Sequences in Cortical Neurons. Science 329, 1671–1675 (2010). 621 18. B. Grothe, M. Pecka, D. McAlpine, Mechanisms of Sound Localization in Mammals. 622 Physiological Reviews 90, 983–1012 (2010). 623 19. J. M. Köhler, Vaccination, Immunity and Breakthrough: Quantitative Effects in 624 Individual Immune Responses Illustrated by a Simple Kinetic Model. Applied Sciences 625 12 (2022). 626 20. P. A. Gagliardi, O. Pertz, The mitogen-activated protein kinase network, wired to 627 dynamically function at multiple scales. Current Opinion in Cell Biology 88, 102368 628 (2024). 629 21. I. Kinjyo, J. Qin, S.-Y . Tan, C. J. Wellard, P. Mrass, W. Ritchie, A. Doi, L. L. Cavanagh, 630 M. Tomura, A. Sakaue-Sawano, O. Kanagawa, A. Miyawaki, P. D. Hodgkin, W. 631 Weninger, Real-time tracking of cell cycle progression during CD8+ effector and 632 memory T-cell differentiation. Nature Communications 6, 6301 (2015). 633 22. A. Stirbet, D. Lazár, Y . Guo, G. Govindjee, Photosynthesis: basics, history and 634 modelling. Annals of Botany 126, 511–537 (2020). 635 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 34 23. J. WU, I. T. BALDWIN, Herbivory-induced signalling in plants: perception and action. 636 Plant, Cell & Environment 32, 1161–1174 (2009). 637 24. X. Yu, B. Feng, P. He, L. Shan, From Chaos to Harmony: Responses and Signaling upon 638 Microbial Pattern Recognition, Annual Review of Phytopathology. 55 (2017)pp. 109–639 137. 640 25. J. M. S. Zegers, I. Irisarri, S. de Vries, J. de Vries, Evolving circuitries in plant signaling 641 cascades. Journal of Cell Science 137, jcs261712 (2024). 642 26. A. K. Block, C. T. Hunter, C. Rering, S. A. Christensen, R. L. Meagher, Contrasting 643 insect attraction and herbivore-induced plant volatile production in maize. Planta 248, 644 105–116 (2018). 645 27. R. Moreno-Sánchez, E. Saavedra, S. Rodríguez-Enríquez, V . Olín-Sandoval, Metabolic 646 Control Analysis: A Tool for Designing Strategies to Manipulate Metabolic Pathways. 647 BioMed Research International 2008, 597913 (2008). 648 28. S. A. R. Mousavi, A. Chauvin, F. Pascaud, S. Kellenberger, E. E. Farmer, GLUTAMATE 649 RECEPTOR-LIKE genes mediate leaf-to-leaf wound signalling. Nature 500, 422–426 650 (2013). 651 29. Y . Blum, J. Mikelson, M. Dobrzyński, H. Ryu, M. Jacques, N. L. Jeon, M. Khammash, 652 O. Pertz, Temporal perturbation of ERK dynamics reveals network architecture of 653 FGF2/MAPK signaling. Molecular Systems Biology 15, e8947 (2019). 654 30. C. Dessauges, J. Mikelson, M. Dobrzyński, M. Jacques, A. Frismantiene, P. A. 655 Gagliardi, M. Khammash, O. Pertz, Optogenetic actuator – ERK biosensor circuits 656 identify MAPK network nodes that shape ERK dynamics. Molecular Systems Biology 657 18, e10670 (2022). 658 31. T. K. Porter, M. N. Heinz, D. J. Lundberg, A. M. Brooks, T. T. S. Lew, K. S. Silmore, V . 659 B. Koman, M. C.-Y . Ang, D. T. Khong, G. P. Singh, J. W. Swan, R. Sarojam, N.-H. 660 Chua, M. S. Strano, A theory of mechanical stress-induced H2O2 signaling waveforms 661 in Planta. Journal of Mathematical Biology 86, 11 (2022). 662 32. T. K. Porter, G. Sánchez-Velázquez, M. S. Strano, The Role of Basal H2O2 663 Concentration in ROS Stress Signaling Waveforms In Planta. ACS Agric. Sci. Technol. 664 5, 1434–1441 (2025). 665 33. P. A. Saa, L. K. Nielsen, Construction of feasible and accurate kinetic models of 666 metabolism: A Bayesian approach. Scientific Reports 6, 29635 (2016). 667 34. L. M. Tran, M. L. Rizk, J. C. Liao, Ensemble Modeling of Metabolic Networks. 668 Biophysical Journal 95, 5606–5617 (2008). 669 35. A. Khodayari, C. D. Maranas, A genome-scale Escherichia coli kinetic metabolic model 670 k-ecoli457 satisfying flux data for multiple mutant strains. Nature Communications 7, 671 13806 (2016). 672 36. H. Link, D. Christodoulou, U. Sauer, Advancing metabolic models with kinetic 673 information. Current Opinion in Biotechnology 29, 8–14 (2014). 674 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 35 37. R. N. Gutenkunst, J. J. Waterfall, F. P. Casey, K. S. Brown, C. R. Myers, J. P. Sethna, 675 Universally Sloppy Parameter Sensitivities in Systems Biology Models. PLOS 676 Computational Biology 3, e189 (2007). 677 38. M. E. Bergman, R. W. J. Kortbeek, M. Gutensohn, N. Dudareva, Plant terpenoid 678 biosynthetic network and its multiple layers of regulation. Progress in Lipid Research 679 95, 101287 (2024). 680 39. J.-K. Weng, The evolutionary paths towards complexity: a metabolic perspective. New 681 Phytologist 201, 1141–1149 (2014). 682 40. J. P. Torres, Z. Lin, J. M. Winter, P. J. Krug, E. W. Schmidt, Animal biosynthesis of 683 complex polyketides in a photosynthetic partnership. Nature Communications 11, 2882 684 (2020). 685 41. M. E. De Obaldia, T. Morita, L. C. Dedmon, D. J. Boehmler, C. S. Jiang, E. V . Zeledon, 686 J. R. Cross, L. B. V osshall, Differential mosquito attraction to humans is associated with 687 skin-derived carboxylic acid levels. Cell 185, 4099-4116.e13 (2022). 688 42. R. Escobar-Bravo, P.-A. Lin, J. M. Waterman, M. Erb, Dynamic environmental 689 interactions shaped by vegetative plant volatiles. Nat. Prod. Rep. 40, 840–865 (2023). 690 43. A. M. Angioy, A. Desogus, I. T. Barbarossa, P. Anderson, B. S. Hansson, Extreme 691 Sensitivity in an Olfactory System. Chemical Senses 28, 279–284 (2003). 692 44. H. Guo, D. P. Smith, Time-Dependent Odorant Sensitivity Modulation in Insects. 693 Insects 13 (2022). 694 45. P. Zu, K. Boege, E. del-Val, M. C. Schuman, P. C. Stevenson, A. Zaldivar-Riverón, S. 695 Saavedra, Information arms race explains plant-herbivore chemical communication in 696 ecological communities. Science 368, 1377–1381 (2020). 697 46. Y . Qin, L. Wang, D. Zhong, Dynamics and mechanism of ultrafast water–protein 698 interactions. Proceedings of the National Academy of Sciences 113, 8424–8429 (2016). 699 47. M. Shibata, H. Nishimasu, N. Kodera, S. Hirano, T. Ando, T. Uchihashi, O. Nureki, 700 Real-space and real-time dynamics of CRISPR-Cas9 visualized by high-speed atomic 701 force microscopy. Nature Communications 8, 1430 (2017). 702 48. E. T. Wurtzel, T. M. Kutchan, Plant metabolism, the diverse chemistry set of the future. 703 Science 353, 1232–1236 (2016). 704 49. J. M. Waterman, T. M. Cofer, O. M. V on Laue, P. Mateo, L. Wang, M. Erb, Leaf Size 705 Determines Damage- and Herbivore-Induced V olatile Emissions in Maize. Plant, Cell & 706 Environment 10.1111/pce.15355 (2025). 707 50. T. T. S. Lew, V . B. Koman, K. S. Silmore, J. S. Seo, P. Gordiichuk, S.-Y . Kwak, M. Park, 708 M. C.-Y . Ang, D. T. Khong, M. A. Lee, M. B. Chan-Park, N.-H. Chua, M. S. Strano, 709 Real-time detection of wound-induced H2O2 signalling waves in plants with optical 710 nanosensors. Nature Plants 6, 404–415 (2020). 711 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 36 51. A. D. Steinbrenner, M. Muñoz-Amatriaín, A. F. Chaparro, J. M. Aguilar-Venegas, S. Lo, 712 S. Okuda, G. Glauser, J. Dongiovanni, D. Shi, M. Hall, D. Crubaugh, N. Holton, C. 713 Zipfel, R. Abagyan, T. C. J. Turlings, T. J. Close, A. Huffaker, E. A. Schmelz, A 714 receptor-like protein mediates plant immune responses to herbivore-associated 715 molecular patterns. Proceedings of the National Academy of Sciences 117, 31510–31518 716 (2020). 717 52. L. Wang, E. Einig, M. Almeida-Trapp, M. Albert, J. Fliegmann, A. Mithöfer, H. 718 Kalbacher, G. Felix, The systemin receptor SYR1 enhances resistance of tomato against 719 herbivorous insects. Nature Plants 4, 152–156 (2018). 720 53. J. M. Waterman, C. R. Hall, M. Mikhael, C. I. Cazzonelli, S. E. Hartley, S. N. Johnson, 721 Short-term resistance that persists: Rapidly induced silicon anti-herbivore defence 722 affects carbon-based plant defences. Functional Ecology 35, 82–92 (2021). 723 54. K. C. Sridhar, N. Hersch, G. Dreissen, R. Merkel, B. Hoffmann, Calcium mediated 724 functional interplay between myocardial cells upon laser-induced single-cell injury: an 725 in vitro study of cardiac cell death signaling mechanisms. Cell Communication and 726 Signaling 18, 191 (2020). 727 55. X. Wang, H. Fang, Z. Huang, W. Shang, T. Hou, A. Cheng, H. Cheng, Imaging ROS 728 signaling in cells and animals. Journal of Molecular Medicine 91, 917–927 (2013). 729 56. T.-O. Kiess, J. Kockskämper, SERCA Activity Controls the Systolic Calcium Increase in 730 the Nucleus of Cardiac Myocytes. Frontiers in Physiology Volume 10-2019 (2019). 731 57. R. V osoughi, Z. Sadeghi Goghari, A. H. Jafari, Modelling system of two insulin-glucose 732 delays to achieve the dynamics of glucose changes. Journal of Biomedical Physics and 733 Engineering 12, 189–204 (2022). 734 58. D. Morris, J. Maclean, A. J. Black, Computation of random time-shift distributions for 735 stochastic population models. Journal of Mathematical Biology 89, 33 (2024). 736 59. Y . Jiang, J. Ye, B. Rasulov, Ü. Niinemets, Role of Stomatal Conductance in Modifying 737 the Dose Response of Stress-V olatile Emissions in Methyl Jasmonate Treated Leaves of 738 Cucumber (Cucumis Sativa). International Journal of Molecular Sciences 21 (2020). 739 60. P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. 740 Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, 741 K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. 742 Polat, Y . Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. 743 Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. 744 van Mulbregt, A. Vijaykumar, A. P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. 745 Scopatz, A. Lee, A. Rokem, C. N. Woods, C. Fulton, C. Masson, C. Häggström, C. 746 Fitzgerald, D. A. Nicholson, D. R. Hagen, D. V . Pasechnik, E. Olivetti, E. Martin, E. 747 Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G. A. Price, G.-L. Ingold, G. E. 748 Allen, G. R. Lee, H. Audren, I. Probst, J. P. Dietrich, J. Silterra, J. T. Webber, J. Slavič, 749 J. Nothman, J. Buchner, J. Kulick, J. L. Schönberger, J. V . de Miranda Cardoso, J. 750 Reimer, J. Harrington, J. L. C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. 751 Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N. J. Smith, 752 N. Nowaczyk, N. Shebanov, O. Pavlyk, P. A. Brodtkorb, P. Lee, R. T. McGibbon, R. 753 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 37 Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. 754 Oshima, T. J. Pingel, T. P. Robitaille, T. Spura, T. R. Jones, T. Cera, T. Leslie, T. Zito, T. 755 Krauss, U. Upadhyay, Y . O. Halchenko, Y . Vázquez-Baeza, SciPy 1.0 Contributors, 756 SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods 757 17, 261–272 (2020). 758 61. M. Abramowitz, Stegun, I.A., Handbook of Mathematical Functions with Formulas, 759 Graphs and Mathematical Tables (Dover, New York, ed. 9). 760 62. A. LoPrete, J. Burge, Bandwidth of Gamma-Distribution-Shaped Functions via Lambert 761 W Function (2025). https://arxiv.org/abs/2509.19307. 762 63. C. Chen, H. Chen, S. Huang, T. Jiang, C. Wang, Z. Tao, C. He, Q. Tang, P. Li, V olatile 763 DMNT directly protects plants against Plutella xylostella by disrupting the peritrophic 764 matrix barrier in insect midgut. eLife 10, e63938 (2021). 765 64. G. M. Gurr, J. Liu, J. A. Pickett, P. C. Stevenson, Review of the chemical ecology of 766 homoterpenes in arthropod–plant interactions. Austral Entomology 62, 3–14 (2023). 767 65. E. Hatano, A. M. Saveer, F. Borrero-Echeverry, M. Strauch, A. Zakir, M. Bengtsson, R. 768 Ignell, P. Anderson, P. G. Becher, P. Witzgall, T. Dekker, A herbivore-induced plant 769 volatile interferes with host plant and mate location in moths through suppression of 770 olfactory signalling pathways. BMC Biology 13, 75 (2015). 771 66. J. Hua, Modulation of plant immunity by light, circadian rhythm, and temperature. 772 Current Opinion in Plant Biology 16, 406–413 (2013). 773 67. N. Guayazán Palacios, T. Imaizumi, A. D. Steinbrenner, The Circadian Clock Regulates 774 Receptor-Mediated Immune Responses to a Herbivore-Associated Molecular Pattern. 775 Plant, Cell & Environment n/a (2025). 776 68. J. M. Waterman, C. I. Cazzonelli, S. E. Hartley, S. N. Johnson, Simulated Herbivory: 777 The Key to Disentangling Plant Defence Responses. Trends in Ecology & Evolution 34, 778 447–458 (2019). 779 69. T. Degen, C. Dillmann, F. Marion-Poll, T. C. J. Turlings, High Genetic Variability of 780 Herbivore-Induced V olatile Emission within a Broad Range of Maize Inbred Lines. 781 Plant Physiology 135, 1928–1938 (2004). 782 70. P.-A. Lin, Y . Chen, D. Chaverra-Rodriguez, C. C. Heu, N. B. Zainuddin, J. S. Sidhu, M. 783 Peiffer, C.-W. Tan, A. Helms, D. Kim, J. Ali, J. L. Rasgon, J. Lynch, C. T. Anderson, G. 784 W. Felton, Silencing the alarm: an insect salivary enzyme closes plant stomata and 785 inhibits volatile release. New Phytologist 230, 793–803 (2021). 786 71. Y . Joo, J. K. Goldberg, L. T. S. Chrétien, S.-G. Kim, I. T. Baldwin, M. C. Schuman, The 787 circadian clock contributes to diurnal patterns of plant indirect defense in nature. 788 Journal of Integrative Plant Biology 61, 924–928 (2019). 789 72. N. G. Palacios, P. Grof-Tisza, B. Behnken, C. M. Arce, D. Wu, A. F. Chaparro, E. A. 790 Schmelz, T. C. J. Turlings, B. Benrey, A. D. Steinbrenner, A plant immune receptor 791 mediates tritrophic interactions by linking caterpillar detection to predator recruitment. 792 bioRxiv, 2025.07.29.667524 (2025). 793 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 38 73. X. Mu, T. D. Evans, F. Zhang, ATP biosensor reveals microbial energetic dynamics and 794 facilitates bioproduction. Nature Communications 15, 5299 (2024). 795 74. C. Chondrogiannis, G. Grammatikopoulos, Transition from juvenility to maturity 796 strengthens photosynthesis in sclerophyllous and deciduous but not in semi-deciduous 797 Mediterranean shrubs. Environmental and Experimental Botany 181, 104265 (2021). 798 75. J. van den Berg, C. Britz, H. du Plessis, Maize Yield Response to Chemical Control of 799 Spodoptera frugiperda at Different Plant Growth Stages in South Africa. Agriculture 11 800 (2021). 801 76. D. Maag, C. Dalvit, D. Thevenet, A. Köhler, F. C. Wouters, D. G. Vassão, J. 802 Gershenzon, J.-L. Wolfender, T. C. J. Turlings, M. Erb, G. Glauser, 3-β-d-803 Glucopyranosyl-6-methoxy-2-benzoxazolinone (MBOA-N-Glc) is an insect 804 detoxification product of maize 1,4-benzoxazin-3-ones. Phytochemistry 102, 97–105 805 (2014). 806 77. W. Feller, An Introduction to Probability Theory and Its Applications (Wiley, ed. 2, 807 1971)vol. 2. 808 78. N. L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distributions (Wiley, 809 ed. 2, 1994)vol. 1. 810 79. D. Kraft, A Software Package for Sequential Quadratic Programming (Wiss. 811 Berichtswesen d. DFVLR, 1988; 812 https://books.google.ie/books?id=4rKaGwAACAAJ)Deutsche Forschungs- und 813 Versuchsanstalt für Luft- und Raumfahrt Köln: Forschungsbericht. 814 80. A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, D. B. Rubin, Bayesian 815 Data Analysis (CRC Press, 2013; 816 https://books.google.ie/books?id=eSHSBQAAQBAJ)Chapman & Hall/CRC Texts in 817 Statistical Science. 818 81. G. H. Golub, C. Reinsch, Singular value decomposition and least squares solutions. 819 Numerische Mathematik 14, 403–420 (1970). 820 82. C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. 821 Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. 822 H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-823 Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, T. E. Oliphant, 824 Array programming with NumPy. Nature 585, 357–362 (2020). 825 83. The Pandas Development Team, pandas-dev/pandas: Pandas, version v2.2.2, Zenodo 826 (2024); https://doi.org/10.5281/zenodo.10957263. 827 84. G. Van Rossum, F. L. Drake, Python 3 Reference Manual (CreateSpace, Scotts Valley, 828 CA, 2009). 829 85. R Core Team, R: A Language and Environment for Statistical Computing., version 4.4.2, 830 R Foundation for Statistical Computing (2024); https://www.R-project.org/. 831 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 39 86. H. White, A Heteroskedasticity-Consistent Covariance Matrix Estimator and a Direct 832 Test for Heteroskedasticity. Econometrica 48, 817–838 (1980). 833 834 Funding 835 This work was supported by the Swiss National Science Foundation (Grant Nr. 201651), the 836 State Secretariat for Education, Research and Innovation (CANWAS), Trinity College Dublin 837 and the University of Bern. 838 839 Author contributions 840 Conceptualization: JMW, GM, LAC, ME 841 Methodology: JMW, GM, LAC, ME 842 Investigation: JMW, GM, LAC, SH, ME 843 Visualization: JMW, GM 844 Funding acquisition: JMW, ME 845 Project administration: JMW, ME 846 Writing (all drafts): JMW, GM, LAC, ME 847 848 Competing Interests 849 The authors have no competing interests to declare. 850 851 Data and Materials Availability 852 All data and code are available upon request from the authors. Additionally, we provide 853 examples of how to adopt the presented approach in Python, R and Excel formats (see 854 Appendix I). 855 856 857 858 859 860 861 862 863 864 865 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 40 Supplemental Tables 866 Supplemental Table 1. Summary of statistical analyses presented in Figure 2 (main text). Bold values: p < 0.05. 867 a = analysed using one-way ANOV A, b = analysed using a Kruskal-Wallis test, * = analysed on log-transformed 868 data, † = analysed using white-adjusted ANOV As. 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 Response F/χ2 p-value df Damage Intensity 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙*a 13.10 < 0.001 5 𝑡*+,'" a 0.37 0.863 5 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 a 4.69 0.004 5 𝑆ℎ𝑎𝑝𝑒 a 2.00 0.115 5 Circadian clock 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙*a 0.502 0.772 5 𝑡*+,'" a 14.77 < 0.001 5 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 b 24.42 < 0.001 5 𝑆ℎ𝑎𝑝𝑒 *a 11.70 < 0.001 5 Oral secretions 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙a 22.84 0.001 1 𝑡*+,'" a 7.49 0.026 1 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 a 94.53 < 0.001 1 𝑆ℎ𝑎𝑝𝑒 a 19.53 0.002 1 Leaf number 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙a 9.18 0.004 2 𝑡*+,'" †a 0.47 0.636 2 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 a 235.36 < 0.001 2 𝑆ℎ𝑎𝑝𝑒 a 126.62 < 0.001 2 Genotype 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙†a 46.98 < 0.001 2 𝑡*+,'" a 115.36 < 0.001 2 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 a 94.85 < 0.001 2 𝑆ℎ𝑎𝑝𝑒 *a 58.92 < 0.001 2 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 41 Supplemental Figure 2. Summary of statistical analyses presented in Figure 3 (main text). Bold values: p < 0.05. a = analysed using one -way ANOV A, b = analysed using a Kruskal-Wallis test, * = analysed on log-transformed data, † = analysed using white -adjusted ANOV As, ‡ = analysed on square -root- transformed data. Supplemental Figure 3. Summary of statistical analyses presented in Figure 4 (main text). Bold values: p < 0.05 and underlined values: p < 0.1. a = analysed using one-way ANOV A, b = analysed using a Kruskal- Wallis test, * = analysed on log -transformed data, † = analysed using white -adjusted ANOV As, ‡ = analysed on square-root-transformed data. Response F/χ2 p-value df Complex damage 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙†a 42.39 < 0.001 2 𝑡*+,'" a 0.99 0.396 2 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛a 1.58 0.239 2 𝑆ℎ𝑎𝑝𝑒a 2.58 0.109 2 Real herbivory 𝑡*+,'" b 18.07 0.001 4 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 †‡a 75.61 < 0.001 4 𝑆ℎ𝑎𝑝𝑒 a 20.95 < 0.001 4 Response F/χ2 p-value df Without oral secretions 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙†a 1.69 0.19 4 𝑡*+,'" ‡a 39.78 < 0.001 4 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 a 95.42 < 0.001 4 𝑆ℎ𝑎𝑝𝑒 *a 34.31 < 0.001 4 With oral secretions 𝐼𝑛𝑡𝑒𝑔𝑟𝑎𝑙a 6.12 0.002 4 𝑡*+,'" *a 1416.00 < 0.001 4 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 †a 274.58 < 0.001 4 𝑆ℎ𝑎𝑝𝑒 †a 291.6 < 0.001 4 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint 42 Supplemental Table 4. Summary of statistical analyses presented in Supplemental figures 1 and 2. Bold values: p < 0.05 and underlined values: p < 0.1. a = analysed using one-way ANOV A, b = analysed using a Kruskal-Wallis test, * = analysed on log-transformed data, † = analysed using white- adjusted ANOV As, ‡ = analysed on square-root-transformed data. Response F/χ2 p-value df Green leaf volatiles 𝑡*+,'" ‡a 0.59 0.584 2 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 a 34.02 < 0.001 2 𝑆ℎ𝑎𝑝𝑒 a 6.78 0.029 2 Gene expression 𝑡*+,'" ‡a 2.18 0.133 3 𝐷𝑢𝑟𝑎𝑡𝑖𝑜𝑛 †a 29.92 < 0.001 3 𝑆ℎ𝑎𝑝𝑒 a 2.95 0.067 3 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Appendix 1. Title: Quantitative modelling of biological response dynamics reveals novel patterns in plant volatile signalling Authors: Jamie M. Waterman*†1,2, Gareth J. Moore†3, Loren K. Amdahl-Culleton3, Sara Hoefer2, Matthias Erb2 Affiliations: 1Discipline of Botany, School of Natural Sciences, Trinity College Dublin, Dublin, Ireland 2Institute of Plant Sciences, University of Bern, Bern, Switzerland 3Independent *Corresponding author. Email: [email protected] †Authors contributed equally to this work .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Practical Guide to Fitting the Response Model in Python, R, and Excel This provides a straightforward workflow to apply the response model to experimental time-series data. Each section (Python, R, Excel) can be read independently. The goal is simple: given a time vector and a response vector, fit the model and extract the derived quantities (shape, duration, total integral). To apply the model, you only need two data vectors: time and response. Provide rough starting values or bounds for the four parameters, then run the fitting procedure in your chosen environment (Python, R, or Excel). The steps are the same everywhere: 1. Load or enter your time and response data 2. Set initial parameter guesses or bounds 3. Run the optimiser (Python: DE→SLSQP; R: DEoptim→nloptr; Excel: Solver) 4. Check the fitted curve against your raw data 5. Use the fitted parameters to compute the derived quantities (shape, duration, and total integral) .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Python Define the Model and Fitting Function This section defines the components needed to run the model: The model is a direct implementation of the form described in the manuscript. A gamma-shaped rise and decay captures the response dynamics, while a logistic onset term smooths the beginning of the curve to improve numerical stability. The helper functions compute the analytical quantities defined in the paper. The main fitting routine fit_response_curve takes three inputs: time, response, and a set of parameter bounds. These bounds can be estimated visually: The fitting proceeds in two steps. A Differential Evolution search explores the full parameter space without relying on good initial guesses. A constrained SLSQP refinement then improves accuracy while ensuring valid parameter ordering (t_onset < t_peak < t_mean). The result is a dictionary containing the four optimised parameters. 1. The model, which evaluates the response for any set of parameters. 2. Helper functions that compute the derived quantities (shape, duration, total integral). 3. The fitting function, which estimates the four parameters from data.t_onset : where the rise begins t_peak : where the maximum occurs t_mean: a point on the decay tail R_peak : approximate peak height import numpy as np from scipy.optimize import differential_evolution, minimize, Bounds from scipy.special import gamma # ----------------------------------------------------- # Model and helper functions # ----------------------------------------------------- defmodel(t, R_peak, t_peak, t_onset, t_mean, epsilon=0.1): alpha =(t_mean - t_onset)/(t_mean - t_peak) beta = t_mean - t_peak t_adj = t - t_onset # np.clip prevents negative or zero values inside the power-law term. # Without this, early-time evaluations can cause (-ve)**fraction -> NaNs. core =((np.clip(t_adj,1e-10,None)/(t_peak - t_onset))**(alpha -1))* \ np.exp(-(t_adj -(t_peak - t_onset))/ beta) # Logistic "ramp" smooths the onset (avoids a discontinuous corner at t_onset) # making the function differentiable for gradient-based optimizers. ramp =1/(1+ np.exp(-(t - t_onset)/ epsilon)) model_pred = R_peak * ramp * core # Replace NaN/inf values that appear with bad parameters. # Forces the model to stay finite, stabilizing optimization. model_pred = np.nan_to_num(model_pred, nan=0.0, posinf=0.0, neginf=0.0) return model_pred .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint defshape(t_onset, t_peak, t_mean): denominator = t_mean - t_onset if denominator ==0: raise ValueError("t_mean and t_onset must not be equal to avoid division by zero.") return(t_mean - t_peak)/ denominator defduration(t_onset, t_peak, t_mean): return t_mean - t_onset deftotal_integral_model(R_peak, t_peak, t_onset, t_mean): alpha =(t_mean - t_onset)/(t_mean - t_peak) beta = t_mean - t_peak t0 = t_peak - t_onset # gamma(alpha) is the Gamma function, generalizes factorial. A = R_peak *(t0 **(1- alpha))* np.exp(t0 / beta) return A *(beta ** alpha)* gamma(alpha) deffit_response_curve(time, response, param_config): param_names =['R_peak','t_peak','t_onset','t_mean'] bounds_list =[(param_config[n]['min'], param_config[n]['max'])for n in param_names] defobjective(x): params =dict(zip(param_names, x)) return np.sum((model(time,**params)- response)**2) # Constraints enforce temporal ordering: # t_onset < t_peak = 10% of mean - onset lambda x: x[1]- x[2]-0.1, # peak must follow onset lambda x: x[3]- x[1]-0.1, # mean must follow peak ] constraints =[{'type':'ineq','fun': c}for c in constraints] try: # Step 1: Differential Evolution (global) de_result = differential_evolution( objective, bounds=bounds_list, strategy='best1bin', maxiter=100, polish=False ) # Step 2: SLSQP (local, constrained) local_result = minimize( objective, de_result.x, method='SLSQP', bounds=Bounds(*zip(*bounds_list)), constraints=constraints, options={'maxiter':2000} .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Fit the model To run the workflow, provide arrays for time and response, set reasonable bounds, and call fit_response_curve. The fitted parameters can then be used to generate a model prediction and to compute shape, duration, and total integral. A plot of raw versus fitted data helps assess quality. A good fit should capture the onset, peak position, and overall decay shape. Large systematic deviations usually indicate that bounds were too restrictive or the data deviate from model assumptions. ) result_x = local_result.x if local_result.success else de_result.x except Exception as e: print(f"Optimization error: {e}") returnNone returndict(zip(param_names, result_x)) import numpy as np import matplotlib.pyplot as plt # ----------------------------------------------------- # Example data # ----------------------------------------------------- time = np.array([ -1.0,-0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5, 4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5, 9.0,9.5,10.0,10.5,11.0,11.5,12.0,12.5,13.0,13.5, 14.0,14.5,15.0,15.5,16.0,16.5,17.0,17.5,18.0,18.5, 19.0,19.5,20.0,20.5,21.0,21.5,22.0 ]) response = np.array([ -0.54,1.062,0.205,0.158,-0.82,3.233,15.25,33.215,47.672,53.405, 57.718,55.81,52.93,53.346,44.75,42.752,36.476,36.371,32.15,29.553, 27.149,23.344,18.599,15.548,14.316,11.673,9.89,8.171,8.574,7.866, 7.178,5.328,5.651,4.508,2.82,2.063,3.402,2.15,-0.062,0.83, 0.435,2.143,0.198,1.669,0.046,-0.256,0.709 ]) # ----------------------------------------------------- # Parameter configuration # ----------------------------------------------------- param_config ={ 'R_peak': {'min':40,'max':80}, 't_peak': {'min':2, 'max':7}, 't_onset': {'min':0, 'max':3}, 't_mean': {'min':5, 'max':12}, } # ----------------------------------------------------- # Fit the curve .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint # ----------------------------------------------------- fit = fit_response_curve(time, response, param_config) R_peak, t_peak, t_onset, t_mean = fit.values() # Compute metrics shape_val = shape(t_onset, t_peak, t_mean) duration_val = duration(t_onset, t_peak, t_mean) integral_val = total_integral_model(R_peak, t_peak, t_onset, t_mean) # Create fitted curve fitted = model(time, R_peak, t_peak, t_onset, t_mean) # ----------------------------------------------------- # Print results # ----------------------------------------------------- print("\n=== Fitted Parameters ===") for k, v in fit.items(): print(f"{k:8s}: {v:8.4f}") print("\n=== Calculated Parameters ===") print(f"Shape: {shape_val:.4f}") print(f"Duration: {duration_val:.4f}") print(f"Integral: {integral_val:.4f}") # ----------------------------------------------------- # Plot results # ----------------------------------------------------- plt.figure(figsize=(10,5)) plt.plot(time, response,'o', label="Data", ms=5) plt.plot(time, fitted,'-', label="Fitted Model", lw=2) plt.xlabel("Time") plt.ylabel("Response") plt.legend() plt.grid(alpha=0.3) plt.title("Model Fit") plt.show() .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Output: Python The printed results show the four fitted parameters and the derived quantities: The accompanying plot confirms whether the fitted curve follows the raw data. R_peak : fitted peak height t_onset, t_peak, t_mean: timing of the main phases Shape, Duration, Integral: derived measures of the response === Fitted Parameters === R_peak :56.0646 t_peak :4.2280 t_onset :1.7252 t_mean :6.8144 === Calculated Parameters === Shape: 0.5082 Duration:5.0892 Integral:388.7458 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint R: Model and Fitting Functions This section provides a complete workflow for fitting the response model in R. It includes three components: The model is a direct implementation of the analytical form described in the manuscript. It combines a gamma-shaped rise and decay with a logistic onset term, which smooths the beginning of the curve and stabilises numerical optimisation. The helper functions use the fitted parameters to compute the analytical descriptors of interest. The fitting routine fit_response_curve takes three inputs: These bounds can be chosen by eye: The optimisation proceeds in two stages. A DEoptim global search explores the parameter space without relying on good initial guesses. A local refinement using optim (L-BFGS-B) then improves accuracy while keeping parameters within their bounds. Simple inequality checks enforce the required temporal order (t_onset < t_peak < t_mean). The function returns all four fitted parameters in a named list. 1. The model function, which computes the predicted response at each time point. 2. Helper functions that calculate the derived quantities (shape, duration, total integral). 3. A fitting function that estimates the four parameters from experimental data. time: numeric vector of time values response: numeric vector of observed responses param_config: named list giving lower and upper bounds for each parameter t_onset: near the initial rise t_peak: near the maximum t_mean: in the decay region R_peak: near the observed peak height library(DEoptim) library(stats) library(pracma) library(ggplot2) # ===================================================== # Model and helper functions # ===================================================== model <-function(t, R_peak, t_peak, t_onset, t_mean, epsilon =0.1) { alpha <-(t_mean - t_onset)/(t_mean - t_peak) beta <- t_mean - t_peak t_adj <- t - t_onset core <-((pmax(t_adj,1e-10)/(t_peak - t_onset))^(alpha -1))* exp(-(t_adj -(t_peak - t_onset))/ beta) ramp <-1/(1+ exp(-(t - t_onset)/ epsilon)) pred <- R_peak * ramp * core pred[!is.finite(pred)]<-0 return(pred) } shape_fun <-function(t_onset, t_peak, t_mean) { .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint denom <- t_mean - t_onset if(denom ==0) stop("t_mean == t_onset") return((t_mean - t_peak)/ denom) } duration_fun <-function(t_onset, t_peak, t_mean) { return(t_mean - t_onset) } total_integral_model <-function(R_peak, t_peak, t_onset, t_mean) { alpha <-(t_mean - t_onset)/(t_mean - t_peak) beta <- t_mean - t_peak t0 <- t_peak - t_onset A <- R_peak *(t0^(1- alpha))* exp(t0 / beta) return(A *(beta^alpha)* gamma(alpha)) } # ===================================================== # Fitting function # ===================================================== fit_response_curve <-function(time, response, param_config) { param_names <- c("R_peak","t_peak","t_onset","t_mean") lower <- sapply(param_config,function(x) x$min) upper <- sapply(param_config,function(x) x$max)

<-function(x) { names(x)<- param_names pred <- model(time, x["R_peak"], x["t_peak"], x["t_onset"], x["t_mean"]) sum((pred - response)^2) } constraint_valid <-function(x) { t_onset <- x[3] t_peak <- x[2] t_mean <- x[4] c1 <-(t_mean - t_peak)-0.1*(t_mean - t_onset) c2 <- t_peak - t_onset -0.1 c3 =0)) } objective_constrained <-function(x) { if(!constraint_valid(x)) return(1e12) objective(x) } de_out <- DEoptim( fn = objective_constrained, lower = lower, upper = upper, DEoptim.control( itermax =100, trace =FALSE ) ) x_de <- de_out$optim$bestmem .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Usage Example: R The following example shows the complete workflow: providing time and response vectors, defining parameter bounds, running the fitting routine, and computing the derived quantities. After fitting, the parameters are used both to generate the model prediction and to calculate shape, duration, and integral. A plot comparing the raw data with the fitted curve provides an immediate check on fit quality. Good fits typically recover the onset, location of the peak, and overall decay structure, though small deviations in noisy regions are expected. penalty_obj <-function(x) { penalty <-if(!constraint_valid(x))1e10else0 objective(x)+ penalty } local_out <- optim( par = x_de, fn = penalty_obj,

="L-BFGS-B", lower = lower, upper = upper, control = list(maxit =2000) ) x_final <- local_out$par names(x_final)<- param_names return(as.list(x_final)) }# ===================================================== # Example data # ===================================================== time <- c( -1.0,-0.5,0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5, 4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5, 9.0,9.5,10.0,10.5,11.0,11.5,12.0,12.5,13.0,13.5, 14.0,14.5,15.0,15.5,16.0,16.5,17.0,17.5,18.0,18.5, 19.0,19.5,20.0,20.5,21.0,21.5,22.0 ) response <- c( -0.54,1.062,0.205,0.158,-0.82,3.233,15.25,33.215,47.672,53.405, 57.718,55.81,52.93,53.346,44.75,42.752,36.476,36.371,32.15,29.553, 27.149,23.344,18.599,15.548,14.316,11.673,9.89,8.171,8.574,7.866, 7.178,5.328,5.651,4.508,2.82,2.063,3.402,2.15,-0.062,0.83, 0.435,2.143,0.198,1.669,0.046,-0.256,0.709 ) # ===================================================== # Parameter config # ===================================================== param_config <- list( R_peak = list(min =40, max =80), .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint t_peak = list(min =2, max =7), t_onset = list(min =0, max =3), t_mean = list(min =5, max =12) ) # ===================================================== # Fit curve # ===================================================== fit <- fit_response_curve(time, response, param_config) R_peak <- fit$R_peak t_peak <- fit$t_peak t_onset <- fit$t_onset t_mean <- fit$t_mean shape_val <- shape_fun(t_onset, t_peak, t_mean) duration_val <- duration_fun(t_onset, t_peak, t_mean) integral_val <- total_integral_model(R_peak, t_peak, t_onset, t_mean) fitted <- model(time, R_peak, t_peak, t_onset, t_mean) # ===================================================== # Print results # ===================================================== cat("\n=== Fitted Parameters ===\n") print(fit) cat("\n=== Calculated Parameters ===\n") cat("Shape: ", shape_val,"\n") cat("Duration: ", duration_val,"\n") cat("Integral: ", integral_val,"\n") # ===================================================== # Plot # ===================================================== df <- data.frame(time = time, response = response, fitted = fitted) ggplot(df, aes(time))+ geom_point(aes(y = response), color ="blue")+ geom_line(aes(y = fitted), color ="red")+ theme_minimal()+ labs( title ="Model Fit", x ="Time", y ="Response" ) .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Output: R The fitted parameters and derived quantities are printed in a straightforward format: The plot confirms visually whether the fitted curve captures the key features of the observed data. R_peak: fitted maximum response t_onset, t_peak, t_mean: timing of the response phases Shape, Duration, Integral: analytical quantities derived from the fitted model === Fitted Parameters === $R_peak [1]56.06604 $t_peak [1]4.228 $t_onset [1]1.725197 $t_mean [1]6.814355 === Calculated Parameters === Shape: 0.5082089 Duration:5.089157 Integral:388.7487 .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Excel: Model and Fitting Functions Excel provides a simple, visual way to use the model. The full workflow is: Step 1: Place your measurements into two columns: Step 2: Reserve cells for the four parameters and their bounds, entering rough starting velues: Step 3: In C2, enter the model formula and drag down: 1. Enter your time and response data 2. Enter initial parameter guesses and provide min/max bounds for each parameter 3. Compute the model prediction 4. Compute squared errors and SSE 5. Use Solver to optimise the parameters 6. Compute the derived quantities (shape, duration, integral) "Our raw data time and response values in the following columns." time = A response = B "Our parameter values in the following cells." R_peak = G2, R_peak_min = H2, R_peak_max = I2 t_peak = G3, t_peak_min = H3, t_peak_max = I3 t_onset = G4, t_onset_min = H4, t_onset_max = I4 t_mean = G5, t_mean_min = H5, t_mean_max = I5 "Model prediction for each time value A2:A…." model = $G$2* ( MAX( A - $G$4,1E-10)/($G$3- $G$4) )^ ( ( $G$5- $G$4)/( $G$5- $G$3)-1)* EXP(-( (A2 - $G$4)-($G$3- $G$4) )/($G$5- $G$3) ) "Model can be used to fill a column with predicted values." model = C .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Column C now contains the predicted values. Step 4: In order to fit the model, we first need a way to measure how far the current parameter guesses are from the actual data. We do this by calculating the squared error between the observed response and the model prediction at each time point. Drag this down to fill the column. To combine all point-wise errors into a single value that Solver can minimise, compute the sum of squared errors (SSE) in G10: The SSE reflects how well the current parameters match the data: lower values indicate a better fit. Solver will adjust the parameters to minimise this value. Step 5: "The square error is calculated from each raw and model response value" sq_err = D =(B - C)^2 "Sum the square errors." SSE = $G$10= SUM(D2:D...) .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Open Data → Solver and configure the following:

Constraints Add both parameter bounds and temporal ordering: Subject to the Constraints: Temporal ordering Set Objective: G10 To: Min Variable cells By Changing Variable Cells: G2:G5 G2>=H2 and G2 =H3 and G3 =H4 and G4 =H5 and G5 <= I5 G4< G3 (onset before peak) G3< G5 (peak before mean) .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint

Click Solve. If successful, Solver updates cells G2–G5 with the best-fitting parameters. Step 6: Once Solver converges, compute the analytical quantities that characterise the fitted response. In any convenient column or cells: These update automatically if the parameter cells change. Select a Solving Method: GRG Nonlinear "Calculated parameter cells filled based on the parameter values." shape =($G$5- $G$3)/($G$5- $G$4) duration = $G$5- $G$4 integral = $G$2* ( ($G$3- $G$4)^(1-( ($G$5- $G$4)/($G$5- $G$3) ) ) )* EXP( ( $G$3- $G$4)/( $G$5- $G$3) )* ( ($G$5- $G$3)^( ($G$5- $G$4)/($G$5- $G$3) ) )* GAMMA( ( $G$5- $G$4)/( $G$5- $G$3) ) .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint Output: Excel After running Solver, the fitted parameter values appear in the parameter cells, and the SSE decreases accordingly. The model prediction column updates automatically, as do the derived quantities (shape, duration, integral). These values provide a direct summary of the fitted response. A plot comparing the observed data with the model prediction gives a visual check of fit quality. A successful fit will show the model curve following the rise, peak, and decay of the measured response, as illustrated in the example screenshot. .CC-BY-NC 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted November 28, 2025. ; https://doi.org/10.1101/2025.11.26.690448doi: bioRxiv preprint