Introduction

Inducible biological responses are transient in nature, changing substantially over relatively short periods of time (1–5). By consequence, real-time measurements are crucial to understand physiological processes, precisely diagnose chronic and acute stress patterns and resolve the mechanisms behind ecological phenomena (6–11). Although biological processes are inherently dynamic, quantitatively describing and comparing temporal features of biological responses remains a challenge. It is well established that static features such as mean concentration at a given instance are useful for general interpretations of stimulus-response relationships, such as the inverse relationship between viral load and effective immune response in humans and the positive relationship between herbivore damage and volatile emissions in plants (12, 13). However, static, nominal values, even if measured at different time points, ignore meaningful, quantifiable features of the dynamics themselves.

Recent examples illustrate the potential of moving beyond measures such as concentration or abundance to generate new biological insights. In plants, variation in the spatial distribution of a toxic metabolite, either within a plant or between plants, can affect herbivore feeding behaviour and thus fitness parameters such as growth, independently from the total concentration of said toxin (14–16). In animals, dendritic signals are sensitive to temporal patterns of synaptic activation; irrespective of the amount of signal, the temporal patterns of signal perception shape spike outputs, which drive crucial functions such as sound localisation (17, 18). Additionally, the interaction between vaccination status, viral exposure and immune responses are all temporally linked, and the effects/response of each will be dependent on how the others change over time (19).

From a mechanistic perspective, the dynamics of biological responses are determined by the interplay of highly regulated signalling events (20–25) and modifications to one or more events can have far-reaching consequences on response dynamics and metabolic outcomes (1, 17, 26– 28). For diffusion-based processes such as generation of reactive oxygen species (ROS) or single-cell signalling events like MAPK cascades, the underlying steps are few and well-characterised, and relatively simple, real-time mechanistic models have been developed (6, 29– 32). Some models of temporal dynamics of metabolism with more complex biosynthetic pathways have also been developed (33, 34). These models are valuable for elucidating biochemical mechanisms and understanding the complexity of metabolic regulation. However, they require substantial a priori biochemical knowledge and the ability to measure numerous kinetic and regulatory parameters (steps), making their construction challenging even for well-characterised metabolic processes (35–37). For more complex and less-well characterised pathways, such as the de novo production of specialised metabolites, complete time-resolved mechanistic models are not currently feasible. Many underlying processes, including enzyme kinetics, precursor and product concentrations, transport delays, degradation pathways and other regulatory processes remain unknown (38–40). Because the biochemistry is not fully defined, the number and nature of parameters required to describe such processes are effectively arbitrary, leading to under-constrained models in which multiple parameter combinations can reproduce the same behaviour and limit reliable fitting and meaningful interpretation (37). Moreover, existing models are often system-specific and must be adapted depending on the metabolites or pathways studied. Consequently, there is no universally standardised framework or parameterisation that can be easily implemented across systems, which limits comparability, broader applicability and ultimately accessibility of dynamics quantification.

Mechanistically informed models may not always be necessary to study biological responses. Volatile chemical signals emitted by organisms shape important interactions such as those between insects and their hosts (41, 42). In this case, the mechanism of volatile formation matters less than when this chemical information is released, how long it persists and whether amounts produced cross perceptual thresholds (10, 43–45). These characteristics are also important for within-organism biochemical processes. Consider protein structural dynamics (i.e., moving between different conformational states in time), for which similar dynamic features would determine interactions with other proteins and metabolites, and thus shape the fate of biochemical processes (46, 47). Thus, extracting time-resolved features of response dynamics that are intuitive, broadly comparable and important for driving biological functions becomes important. Both statistical and theoretical methods for modelling dynamic biological responses have been proposed (8, 13, 19). However, a universal model that enables quantitatively rigorous and biologically informative comparison of dynamic response features, either across stimuli or across systems, has yet to be developed.

Motivated by the need to extract biologically meaningful information from dynamic responses, and the lack of a robust framework to do so, we set out to develop a method that allows for biologically meaningful and consistent comparisons of dynamic responses without explicit knowledge of the underlying physiochemical and biochemical processes. We used our model quantify the dynamics of diverse plant volatile responses to different stimuli to identify and evaluate whether meaningful differences exist in the dynamics themselves, generating novel biological insights. We considered plant responses to stress to be an optimal test case, as plants have evolved a highly complex array of responses to environmental stimuli, that cover a broad range of physiological and ecological functions (48), and show starkly distinct kinetics (42). Induced volatiles in particular can be measured non-destructively, enabling diverse processes to be observed in real time (7). Finally, these volatiles provide important chemical information to the wider environment, and we thus assume that the generated insights will be ecologically informative (10, 45).

Through this approach, we uncovered that all measured responses conform to a common model structure, while the parameter values varied considerably between stimuli and volatiles. Not only would these phenomena go undetected without a robust quantitative framework, but they reveal that biologically meaningful information is encoded in induced response dynamics and that temporal response structure carries functional significance. Importantly, as the model structure is not specific to either response type or timeframe, it will be applicable over a broad range of dynamic processes and biological systems, thus unlocking novel comparative and integrative approaches across the tree of life.

Results

Model development

To build a phenomenological model that captures dynamic responses in an unbiased manner, we modelled responses as a distribution of waiting times. We asked the question: How long does it take for a response unit (e.g. a volatile molecule) to appear following a stress event? In doing so, we interpret each measured unit as a probabilistic occurrence with a certain delay, and thus model the response curve as a probability density function.

This framing is advantageous for three reasons: first, the model yields a relatively small number of parameters that are both identifiable (i.e. are intuitive by eye) and interpretable (i.e. have biological meaning). This minimises overfitting and, as a result, makes the model robust to sample-to-sample variation. Secondly, it enables comparability across responses which involve very different mechanisms and timeframes. Thirdly, the same framework can be applied to any induced phenotype that can be measured quantitatively over time, such as gene expression, hormone accumulation, enzymatic conversions, volatile emissions and resistance induction, even if they show starkly different dynamics and biosynthetic processes, as they all typically exhibit the same behaviour: induction delay following stimulation and a rapid increase in response, which, after reaching a peak, gradually declines (1, 2, 28, 49–57). This is mathematically intuitive as these responses are built from a sequential chain (or combinatory network) of activation, synthesis, transport and eventually measurement/detection. Each step adds stochastic waiting time with their sum naturally producing a skewed, unimodal distribution (58). Modelling this kind of process usually involves treating each biochemical step as a stochastic waiting time drawn from a defined probability distribution and then combining these into a final waiting-time distribution which captures the total observed response dynamics. As such, we used a reparametrised gamma distribution, which retains its form when steps are combined and can thus describe both early and later phases of the response in a consistent way (see materials and methods for full model details):

where:

• R_peak is the maximum measured response,

• t_onset is the onset delay,

• t_peak is the time of peak response,

• t_mean is the mean response time.

The fitting parameters can then be used to calculate additional characteristic parameters (Fig 1A).

Figure 1.

Total response (Integral) – The integral is often difficult to measure experimentally; due to time limitations in the measurement processes, entire curves are not resolved and thus incomplete integrations are used (59). When fit onto the experimental data, our model can be used to calculate the theoretical integral of the response curve even if it is incomplete:

where Γ() is the Gamma function, a continuous generalisation of the factorial (46, 47).

Duration

For response dynamics, response length is typically defined as ‘broadness’, for example as full width at half max (FWHM) (6). The FWHM of the Gamma function has no analytical solution and must be determined numerically (62). Additionally, FWHM, does not explicitly capture the start of the response. For these reasons Duration is used as a simple alternative, achieving a similar descriptive value and used as a normalisation value in the ‘shape’ feature below.

Shape

The ‘shape’ of a response curve is an abstract feature. Since our model is intended to be used over variable time scales, we developed a time-normalised (by Duration) metric that allows for comparison of curve shapes independently from response length. Further, we aimed to understand the symmetry of the curve in response to different stimuli or across plant genotypes/species to broadly assess temporal nuances. Shape presents a new phenotypic axis which may link to important biochemical or physiochemical properties that regulate induced responses. Shape is a normalised value, and when responses follow gamma distributions it falls between zero and one; a Shape of zero indicates a perfectly symmetrical curve (t_mean = t_peak) and a shape of one represents a completely right-skewed curve (t_peak= t_onset):

Taken together, the fit and derived parameters give a framework for comparing curves, particularly those that appear distinct by eye but are difficult to distinguish quantitatively (Fig 1B).

Uncovering novel patterns in plant volatile responses

We applied the model to inducible volatile emissions in maize (Zea mays) in real time by PTR-ToF-MS following different stress treatments. First, we focused on the homoterpene 4,8-dimethylnona-1,3,7-triene (DMNT) as a highly inducible, ecologically relevant plant volatile (63–65). We tested whether the amount of mechanical leaf damage influences DMNT induction. Earlier work showed that damage intensity strongly correlates with terpene emissions (12). Indeed, higher damage increased total DMNT emissions DMNT (Fig. 2A), and our model faithfully captured this pattern through Integral (Fig 2B). No clear patterns were detected for the other parameters, apart from differences in Duration, which was marginally longer for intermediate amounts of wounding (Fig. 2D).

Figure 2.

Next, we tested the influence of the circadian clock on wound-induced DMNT by wounding plants at different times of day under continuous light. The circadian clock is known to play a role in regulating how plants respond to stress, including transcriptional regulation of defence genes (66, 67). Time of day had no impact on the overall amount of DMNT emission (Fig 2G). However, t_onset, Duration and Shape all varied significantly; wounding in the evening resulted in a significantly more rapid, shorter DMNT burst – a novel pattern that has not been reported before (Fig. 2F-J).

We also tested the influence of insect oral secretions (OS), which contain both elicitors and effectors which can increase and decrease responses compared to wounding alone, respectively (68). We confirmed that OS enhance the total emission of DMNT (Fig 2K-L). Interestingly, OS treatment also led to a more rapid and prolonged response, resulting in a significantly different curve shape (Fig 2M-O).

Leaf size and developmental stage can influence total volatile emissions (49), but whether these parameters also influence response curves in other ways is unknown. We thus reanalysed the dataset from our earlier work with our model. We found that leaf 3 (largest size and intermediate age) emits volatiles for the longest period (Fig 2S). Interestingly, leaf 3 did not produce more volatiles than leaf 2 (older leaf) but produced them for a longer period of time (Fig 2Q and S). Both t_onset and Shape followed a clear developmental gradient, whereby older leaves exhibited slower and more symmetrical emission dynamics compared to younger leaves (Fig 2R and 2T).

Different genotypes vary strongly in the amount and type of volatiles they emit (69), but whether response curves are also different is unknown. We determined response curves in three maize inbred lines that differ in their capacity to produce DMNT (Fig 2U). As expected, Integral varied for all genotypes (Fig 2V). Interestingly, we found that higher emitting genotypes also had substantially earlier t_onset (Fig 2W), however higher emissions do not necessarily equate to earlier onsets (Fig 2A-C). Additionally, the two highest emitting genotypes, CML287 and NC300, had similar Duration (Fig 2X) and CML287 had the most symmetrical dynamics. Thus, different genotypes show starkly different response patterns, some of which are independent of total volatile quantity. Taken together, these experiments illustrate the power of our approach to uncover novel, genetically determined response patterns to environmental stimuli.

Unravelling differences between biochemically distinct compound classes

To explore how different volatiles response to the same stimuli, we characterised volatiles from different biosynthetic pathways, including terpenes, indole and green leaf volatiles (1). In response to wounding, the different volatiles exhibited significantly different response curves, as described before, which was visible in differences in t_onset, Duration and Shape (Fig. 3A-E; Fig S1). Interestingly, the differences became more pronounced with the application of OS (Fig. 3F-J). This pattern was evident across all parameters but was strongest for Duration, whereby in the presence of OS each compound was produced for a unique amount of time, which was not the case in the absence of OS (Fig 3E and 3J). This illustrates how OS components differentially modulate volatile pathways compared to mechanical wounding across volatile groups.

Figure 3.

To explore this phenomenon more deeply and test our model on a less temporally resolved dataset, we modelled the wound-responsive gene expression dynamics of ZmCYP92C5, ZmIGL, ZmTPS2 and ZmTPS10 (Fig S2), which are coding for the rate limiting enzymes of the biosynthesis pathways of DMNT/TMTT, indole, monoterpenes and sesquiterpenes (49). Trends in Duration of the expression of these genes matched emission of corresponding volatiles (Fig S2C and Fig 3E), confirming that volatile emission is, at least in part, regulated by biosynthetic limitations (Fig S2). There were clear differences in t_onset detected (Fig S2B), however the lack of temporal resolution in measurements of gene expression dynamics makes clear interpretations of t_onset difficult (Figs S3 and S4). Interestingly, although clear, differences in Shapedid not match corresponding emission values (Fig S2D), suggesting unique parameterisation across levels of organisation. However this may, again, be partially explained by low resolution.

Exploring complex damage and response patterns

To characterize more complex induction patterns, we quantified volatile emissions following multiple wounding events (Fig 4A). The current assumption from the literature is that subsequent wounding should lead to stronger defence responses beyond cumulative effects, as the earlier wounding will prime plants for subsequent responses. However, such effects have been hard to isolate for wounding events that follow each other closely in time due to overlapping response curves. Thus, we modelled the responses to three wounding events as the sum of three separate curves, all of which being effectively incomplete or not fully resolved in some way. By decomposing the overall emission into separate fitted curves, we were able to quantify the dynamic contribution of each individual wounding event, even when the peaks overlapped. This revealed a clear priming effect where the integral of the second and third responses were larger than that of the initial peak, demonstrating that prior damage enhanced subsequent emissions beyond additive effects (Fig 4B). At this particular time interval, multiple damage events did not significantly modify other fit parameters (Fig 4C-E). The functionality of this model on more complex damage patterns highlights its utility to explore response dynamics under complex stress regimes.

Figure 4.

Herbivores do not make single wounds when they feed, but rather take many bites over time. As such, we quantified emission responses to feeding by a single Spodoptera exigua larva for 30 min. Natural herbivory patterns generated responses that could be easily fit for all groups of volatiles (Fig 4F). t_onset showed some similarities to wounding response patterns, namely that sesquiterpenes and TMTT take longer to begin emitting than other compounds (Fig 4G). Interestingly, Duration values more closely matched those from mechanical wounding than wounding + OS, with fewer unique durations forming (Fig 4H). Shape values showed stark differences compared to wounding treatments, suggesting that damage patterns may play an important role in determining response curves (Fig 4I). The ability to effectively fit parameters to responses induced by real herbivore feeding indicates that our model can be used to detect differences in response dynamics even when exact damage patterns and timing are not known, and thus it can be used widely.

Quantifying incomplete curves

Understanding the dynamics of defence responses requires some level of continuous monitoring. However, how the temporal resolution of measurements impacts the capacity to understand dynamic patterns is not well known. We tested how measurement duration and sampling resolution affected quantification. Firstly, in order to simulate a shorter measurement window, we removed datapoints from the end of the curve (in time) and compared fit parameters to those from the complete dataset. For integral and Duration, there was a consistent point around t_mean where removal of data dramatically reduce the accuracy of model parameters (Fig S3B, S3D, S3G, S3I, S3L and S3N). t_onset was more sensitive to data removal for indole and TMTT than it was for DMNT (Fig S3C, H and M). Similarly for shape, which is t_onset-dependent), indole and TMTT fits broke down relative to the complete dataset with less data removal than DMNT (Fig S3E, J and O). Interestingly, for DMNT and indole, even with removal of ca. 10-15 hr of curve the fitted parameters remained almost entirely stable, if not identical, suggesting a substantial degree of flexibility and potential to understand the dynamics of only partially resolved peaks.

Second, we tested the resolution tolerance of the model by periodically removing datapoints, simulating a lower resolution measurement (Fig S4A, F and K). With periodic removal of up to 85% of data points, the emission integral was able to be recovered with accuracy and precision for all compounds, suggesting that a coarse ‘outline’ of temporal response patterns can be sufficient to accurately estimate total emissions over time (Fig S4B, G and L). For indole, removal of 70% data in periodic intervals resulted in breakdown of integral as the overall curve structure was nearly entirely lost (Fig S4B). Indole has the fastest emission dynamics, and as such, indole has the lowest effective resolution in comparison to DMNT and TMTT, which both have slower kinetics and thus more data points between t_onset and t_mean. Interestingly, for other fit parameters (t_onset, Duration and Shape), accuracy and precision began to consistently breakdown at 70% removal, suggesting there is a limit to measurement frequency for certain dynamic features (Fig S4). Importantly, these results highlight that response kinetics are a necessary consideration to determine the use of truly ‘continuous’ measurements; they may not always be required or even provide additional dynamic information.

Discussion

Dynamic responses to environmental stimuli are a universal property of life. Yet, unifying approaches to characterize such responses are lacking. Across biology, we have major gaps in our understanding of patterns and differences in dynamic biological responses. Induced volatile emissions are an illustrative example in this context. In a series of experiments, we find that our model reveals novel differences in volatile induction along circadian, developmental and genetic axes, and in response to wounding and insect-specific molecular patterns. Additionally, we uncouple convoluted response patterns to identify rapid priming events that can inform how complex stress regimes might change defences trajectories compared to simple ones. Among the most noteworthy discoveries is that, independent of light, the time of day of wounding has a strong impact on the onset, duration and shape of the volatile induction curve but not the total response intensity. This provides a window into the tight regulation of the speed and duration of induced volatile emissions by the circadian clock (10, 70), providing new phenotypes that link clock regulation to environmental interactions. Equally noteworthy is the finding that insect oral secretions (OS) enhance differences in response curves between different volatile classes, thus resulting in unique temporal volatile fingerprints. This newly discovered pattern may explain OS-specific responses that affect interactions with herbivores and herbivore natural enemies beyond differences in overall volatile quantities or timing (71). Finally, we uncover that, contrary to current expectations, the quantity and duration of induced volatile emissions can be regulated independently of each other across plant genotypes. This finding paves the way to identify regulators of volatile induction duration through genetic approaches, potentially leading to the identification of new “late signalling components” involved in sustaining defence responses. Further, these novel dynamic features might ultimately be leveraged to non-destructively distinguish and identify stimuli, for example in an agricultural context, where many individual induced responses may be shared between multiple stressors, and total response is thus insufficiently informative (6, 45, 50).

Recent advances in continuous, real-time monitoring have greatly improved our ability to measure biological response dynamics (5, 7, 10, 50), yet analyses of such data remain limited. Our modelling approach addresses this gap by providing a standardised framework for quantifying response dynamics in a way that enables biologically and ecologically meaningful comparisons across systems. The model improves on traditional response quantification by providing true integrals and maximum response values as well as additional response parameters that are not easily extractable otherwise. In this study, we use volatile emissions and gene expression from plants as representative measures of broad biological responses. Thanks to the robustness and generality of the model, it will be straightforward to apply it across biological systems, from energy dynamics in microbes to electron transport in plants to neurological and immune responses in animals (13, 17, 72, 73). Thereby, the framework opens opportunities for ambitious endeavours, such as tree-of-life scale meta-analyses, to uncover shared and divergent principles of organismal responses to environmental stimuli and to reveal new fundamental biological phenomena. Furthermore, coupling model parameters with machine-learning algorithms or constrained optimisation solvers could enable the reconstruction of stimulus patterns from observed dynamics, or the prediction of dynamic responses from environmental inputs. We propose this framework as a universal standard for intuitive, interpretable and biologically grounded quantification of dynamic responses, to advance understanding of fundamental processes and to guide innovation across medicine, agriculture and beyond.

To ensure that the model can be used widely, we generated a set of freely accessible resources that can be implemented in Python, R and Excel formats (Appendix I), and applied broadly to time-resolved response data to yield all relevant readouts for further statistical analyses.

Materials and Methods

Plants and insects

Plant growth conditions were identical to those used in Waterman et al. (2025). In brief, V2-stage (two developed leaves, one expanding leaf and one emerging leaf) maize (Zea mays) seedlings were used throughout the study. At this stage maize plants are particularly susceptible to agricultural pests such as Spodoptera spp. (74). Plants were grown in commercial potting soil (Selmaterra, BiglerSamen, Switzerland) in 180 mL pots under greenhouse conditions and supplemented with artificial lights (ca. 300 µmol m−2 s−1). The greenhouse was kept at 22 ± 2°C, 40%–60% relative humidity, with a 14 h: 10 h, light: dark cycle.

Spodoptera exigua larvae (Frontier Agricultural Sciences, USA) were reared from eggs on an artificial diet (75). At least 24 hr prior to experimental treatments larvae were fed B73 maize leaves. Oral secretions (OS) were collected by probing the mouths of larvae with a pipette tip and stored at -20ºC until use.

Experimental Treatments

All wounds were inflicted using haemostat forceps, an established method of simulating herbivory (49, 68). The maize inbred line, B73, was used unless otherwise stated. The standard damage intensity was 120 mm² and time of initial damage was 11:00 unless otherwise stated. To understand how different damage regimes impact induced response dynamics we conducted several experiments spanning a range of treatments:

Variable damage intensity

We damaged plants at 20, 40, 80, 120, 160 and 320 mm² on the third-oldest leaf (leaf 3). Damage patterns always encompassed the central segment of each leaf, and any damage above 40 mm² was evenly distributed across the base, middle and tip of leaf 3 (49).

Time of day

We damaged plants at 11:00, 11:45, 12:30, 13:15, 14:45, 15:45 and 23:00.

Oral secretions (OS)

Damaged plants were either treated with milliQ water or 50% OS from late-instar Spodoptera exigua larvae.

Leaf developmental stage

At the V2 stage maize has four leaves: two are fully developed (leaf 1 and leaf 2), one is actively expanding (leaf 3) and one is emerging (leaf 4). We used the raw emission data from Waterman et al. (2025) to measure emission dynamics in each leaf. The damage treatments in the previous study were 60mm² total, and damage was inflicted in three bouts of 20 mm² damage in ca. 30 min intervals

Genotype

To explore the potential variation in response dynamics across maize genotypes we damaged three additional maize inbred lines, with three distinct volatile emission capacities: MO17 (low emission), NC300 (intermediate-high emission) and CML287 (high emission).

Real herbivory

3^rd instar Spodoptera exigua were placed on leaf 3 and were left to feed for 30 minutes.

Volatile sampling

Volatiles were collected as in Waterman et al. (2025). Briefly, entire seedlings were placed in transparent glass chambers (Ø×H 12 × 45 cm) that were sealed other than a clean airflow inlet and an outlet. Clean air was supplied at a flow-rate of 0.8 L min⁻¹. Volatiles were measured with a high-throughput platform comprising of a proton transfer reaction time-of-flight mass spectrometer (PTR-ToF-MS; Tofwerk, Switzerland) and a custom-made automated headspace sampling system. The outlet of the chamber was accessible to the autosampler/PTR-ToF-MS system, which drew air at 0.1 L min⁻¹. Between samples, a zero-gas measurement was performed for 3 s to flush the system. At each time point (as indicated by the x-axis of the respective figure), volatiles were continuously measured for 25–30 s and averaged to a single mean per sample. Complete mass spectra (0–500 m/z) were recorded in positive mode at ca. 10 Hz. The PTR was operated at 100°C and an E/N of approximately 120 Td. The volatile data extraction and processing were conducted using Tofware software package v3.2.2 (Tofwerk, Switzerland). Protonated compounds were identified based on their molecular weight + 1. During volatile collection, LED lights (DYNA, heliospectra, Sweden) were placed ca. 80 cm above the glass chambers and provided light at 300 μmol m−2 s−1. Identical light:dark cycle timing as in the greenhouse for plant growth was used.

Gene expression

Leaf 3 tissue was harvested 0.5, 1, 1.5, 2, 3 and 5 hr after damage (see above for damage treatment) and flash frozen on liquid nitrogen. Total RNA extraction and purification, genomic DNA removal, cDNA synthesis and quantification of gene expression were conducted identically to Waterman et al. 2024. Quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed using ORA SEE qPCR Mix (highQu GmbH, Germany) on an Applied Biosystems QuantStudio 5 Real-Time PCR system. The normalised expression (NE) values were calculated as in Waterman et al. (2024) using ubiquitin (UBI1) as the reference gene. Gene identifiers and primer sequences are listed in Table S1.

Modelling emission dynamics

For all experiments, volatile emissions were normalised by the biomass of leaf 3 (damaged leaf). This is because we have previously shown that the size of the damaged leaf is limiting for volatile emissions (49). Additionally, values were normalised to the maximum response observed in each experiment, yielding a range of positive values < 1.

Choosing the right distribution

To model the waiting time probability distribution, we used a Gamma distribution:

Where:

• α is the number of steps,

• β is the rate parameter (for each step),

• Γ(α) is the Gamma function evaluated at α

This distribution arises from the sum of several exponential waiting times, making it suited to processes where multiple sequential steps precede a single observable outcome (76). Unlike the log-normal or inverse Gaussian, the Gamma distribution is closed under addition, meaning that the sum of sequential Gamma is itself a Gamma distribution (77). This allows modelling of both early and late phase responses, even if they are both part of the same chain of events. In the formal definition, α represents the effective number of steps and β an effective average delay, though we do not assume that these biological processes are truly memoryless and independent as required by the strict derivation. Instead, we treat the Gamma distribution as a phenomenological approximation of the process, which fits the shapes of the responses we (and others) observe well (1, 2, 28, 49–57). With this in mind, we make a number of meaningful reparameterisations to the standard Gamma distribution, as follows.

We begin with the standard Gamma distribution and introduce a shifted time variable t_adj to account for a delay in the onset of the curve:

The Gamma probability density function is then given by:

Our next step is to relate the variables to observable quantities, and so we want to express α and β in terms of measurable features of the response curve.

For the Gamma distribution:

• Mode = (α − 1)β = t_peak − t_onset

representing the time at which the response reaches its peak.

• Mean = αβ = t_mean − t_onset

representing the average time after onset that the response occurs.

In order to solve for α and β for these observable features we: Subtract the first from the second:

Simplify:

Substitute into αβ = t_mean − t_onset:

The Gamma probability density function is normalised to integrate to 1.

However, in experimental data the absolute amplitude (the response intensity, R_peak) is an observable quantity and so we define the response model:

where R_peak is the observed response at t_peak.

This ensures that the curve shape follows the Gamma form, and that the maximum of R(t) equals R_peak.

The next step is to solve for f-t_peak; α, β, t_onset.

Substitute t = t_peak into f(t; α, β, t_onset):

Replace t_peak − t_onset with (α − 1)β using the mode definition:

Simplify the exponent:

Final expression:

We can therefore define:

Substitute f(t)and (f-t_peak). into the ratio:

Cancel the common terms β.Γ(α):

Simplify the exponent:

Replace α and β with expressions in terms of (t_onset, t_peak, t_mean) to give the final model:

Where:

• R_peak is the maximum measured response,

• t_onset is the onset delay,

• t_peak is the time of peak response,

• t_mean is the mean response time.

Background subtraction

Since low levels of plant volatiles are constantly emitted even without wounding or other measurable stimuli, a background curve was subtracted so that only induced emissions were modelled (42). This background was estimated from undamaged plants and either directly subtracted or scaled to the pre-damage baseline to correct for instrumental drift or individual plant batch variation.

Fitting

Each individual curve was fit using a two-stage optimisation approach. In order to avoid local minima, we first used a differential evolution optimisation followed by a local refinement step using SLSQP (78). Since the emission data were individually relatively noisy, we reconstructed the dominant temporal trend using singular value decomposition (SVD) and incorporated this as a weak prior to guide the fitting algorithm toward biologically realistic solutions (79, 80). A number of constraints and bounds were enforced in order to guide the fitting of the data. Temporal ordering was imposed such that t_onset < t_peak < t_mean. To prevent degenerate solutions, minimum separations were required, with t_peak − t_onset > 0.1 and t_mean − t_peak > 0.1. Bounds were dynamically set for each curve; for instance, t_peak was initialised at the time of maximum observed response and restricted to within ±20% of this value, while R_peak was constrained between 80% and 120% of the observed maximum. t_onset and t_mean times were loosely defined based on empirically informed ranges so fitted values were virtually unconstrained. All fitting, calculations and data manipulation were done using SciPy (60), NumPy (81) and Pandas (82) in Python 3.12.4 (83).

Model testing

To explore the robustness of the model and fitting methods to the effects of resolution and measurement time, a data set was taken and data were removed in order to simulate a lower quality measurement. We chose the dataset of damage plants treated with insect OS, as a representative and ecologically relevant subset. In a first test, data from the end of the measurement were removed in a stepwise manner to simulate a shorter (incomplete) experimental measurement period. Secondly, data were systematically down-sampled to simulate a lower resolution measurement. After each step of the truncation or down sampling the data was fit and the fit parameters compared.

Additional statistical analyses

All further statistical analyses were conducted in R version 4.4.2 (84). Differences in model outputs between groups were determined using one-way ANOVA. Where necessary, data were transformed to meet ANOVA assumptions. Where necessary, to obtain heteroscedasticity-consistent standard errors, White-adjusted ANOVAs were used (85).Where data did not meet normality assumptions, Kruskal-Wallis tests were used. Where data did not meet the homogeneity of variance and normality assumptions, even following transformations, differences were determined using Welch’s ANOVA. Statistical test summaries are included in supplemental tables 1-4.

Supplemental material

Supplemental Figure 1.

Supplemental Figure 2.

Supplemental Figure 3.

Supplemental Figure 4.

Supplemental Table 1.

Supplemental Table 2.

Supplemental Table 3.

Supplemental Table 4.

Data availability

The Python scripts used to model response dynamics and the datasets are available at the GitHub repository: https://github.com/watermaj/Quant-dynamics. Additionally, we provide examples of how to adopt the presented approach in Python, R and Excel formats (see Appendix I).

Acknowledgements

We would like to thank Sarah Holder for assistance with plant growth and Dr. Tristan Cofer for experimental assistance. This work was supported by the Swiss National Science Foundation (Grant Nr. 201651), the State Secretariat for Education, Research and Innovation (CANWAS), the University of Bern and Trinity College Dublin.

Additional information

Author contributions

Conceptualization: JMW, GM, LAC, ME

Methodology: JMW, GM, LAC, ME

Investigation: JMW, GM, LAC, SH, ME

Visualization: JMW, GM

Funding acquisition: JMW, ME

Project administration: JMW, ME

Writing (all drafts): JMW, GM, LAC, ME

Funding

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (SNF) (201651)

• Jamie Waterman

WBF | Staatssekretariat für Bildung, Forschung und Innovation (SBFI) (CANWAS)

• Matthias Erb

Appendix 1. 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)

Python

Define the Model and Fitting Function

This section defines the components needed to run the model:

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.

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:

• 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

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.

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.

Output: Python

The printed results show the four fitted parameters and the derived quantities:

• R_peak : fitted peak height

• t_onset, t_peak, t_mean : timing of the main phases

• Shape, Duration, Integral : derived measures of the response

The accompanying plot confirms whether the fitted curve follows the raw data.

R: Model and Fitting Functions

This section provides a complete workflow for fitting the response model in R. It includes three components:

• The model function, which computes the predicted response at each time point.

• Helper functions that calculate the derived quantities (shape, duration, total integral).

• A fitting function that estimates the four parameters from experimental data.

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:

• time : numeric vector of time values

• response : numeric vector of observed responses

• param_config : named list giving lower and upper bounds for each parameter

These bounds can be chosen by eye:

• t_onset : near the initial rise

• t_peak : near the maximum

• t_mean : in the decay region

• R_peak : near the observed peak height

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.

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.

Output: R

The fitted parameters and derived quantities are printed in a straightforward format:

• 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

The plot confirms visually whether the fitted curve captures the key features of the observed data.

Excel: Model and Fitting Functions

Excel provides a simple, visual way to use the model. The full workflow is:

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)

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:

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

Open Data → Solver and configure the following:

Objective

• Set Objective :

• To : Min

Variable cells

• By Changing Variable Cells : G2:G5

Constraints

Add both parameter bounds and temporal ordering:

Subject to the Constraints :

• G2 >= H2 and G2 <= I2

• G3 >= H3 and G3 <= I3

• G4 >= H4 and G4 <= I4

• G5 >= H5 and G5 <= I5

Temporal ordering

• G4 < G3 (onset before peak)

• G3 < G5 (peak before mean)

Method

• Select a Solving Method : GRG Nonlinear

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.

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.