1. Introduction
Climate change poses significant challenges to global food security due to its negative impact on crop growth and productivity (Suter & Widmer, Reference Suter and Widmer2013). With a growing world population, especially in climate-vulnerable areas (Dai, Reference Dai2011), it becomes essential to develop crop varieties that are resilient to stresses, such as pests, heat and drought, particularly when these occur simultaneously or in rapid succession. To understand these plant environmental response mechanisms, as well as their potential trade-offs, we aim to understand the underlying gene regulatory networks (GRNs) (Quint et al., Reference Quint, Delker, Franklin, Wigge, Halliday and van Zanten2016), i.e., sets of genes which influence each other’s expression. The temporal dynamics of these GRNs are commonly modelled using ordinary differential equations (ODEs) (Cao et al., Reference Cao, Qi, Zhao, Wang, Tan and Tian2012; de Luis Balaguer et al., Reference de Luis Balaguer, Fisher, Clark, Fernandez-Espinosa, Möller, Weijers, Lohmann, Williams, Lorenzo and Sozzani2017; Valentim et al., Reference Valentim, van Mourik, Posé, Kim, Schmid, van Ham, Busscher, Sanchez-Perez, Molenaar, Angenent, Immink and van Dijk2015) fit to time-series datasets. Importantly, such ODE models can reflect biological prior knowledge and enable causal insights through formal mathematical model analysis and human interpretation. They have successfully been used to describe the influence of a small group of genes on flowering time (Chew et al., Reference Chew, Seaton, Mengin, Flis, Mugford, George, Moulin, Hume, Zeeman, Fitzpatrick, Smith, Stitt and Millar2022) or hypocotyl elongation (Nieto et al., Reference Nieto, Catalán, Luengo, Legris, López-Salmerón, Davière, Casal, Ares and Prat2022) in response to different temperature and light regimes, and enhanced our mechanistic understanding of these processes.
However, these models typically only include a small number of genes (<100) (Fröhlich et al., Reference Fröhlich, Kessler, Weindl, Shadrin, Schmiester, Hache, Muradyan, Schütte, Lim, Heinig, Theis, Lehrach, Wierling, Lange and Hasenauer2018) because they cannot (yet) be constructed for genome-scale plant stress response networks due to a lack of extensive, mechanistic, biological knowledge. On top of that, a large number of input genes makes parameter estimation prohibitively complex, as this requires significant amounts of data (Cao et al., Reference Cao, Qi, Zhao, Wang, Tan and Tian2012; Segal et al., Reference Segal, Shapira, Regev, Pe’er, Botstein, Koller and Friedman2003; Shaik & Ramakrishna, Reference Shaik and Ramakrishna2013) and results in exponentially more possible gene network topologies. To overcome these issues and enable ODE modelling of large-scale regulatory networks from time-series data, we focus on gene modules (also known as gene clusters or gene sets): groups of co-expressed and/or co-regulated genes (Saelens et al., Reference Saelens, Cannoodt and Saeys2018; Segal et al., Reference Segal, Shapira, Regev, Pe’er, Botstein, Koller and Friedman2003). By modelling temporal behaviour of gene modules instead of individual genes, we significantly reduce the number of interactions that require parameter estimates. Additionally, gene modules can help bridge the gap between single-gene expression levels and high-level phenotypes by linking these modules to broader biological processes or well-studied plant hormones. Finally, modules improve our ability to convert novel knowledge from model plants, such as A. thaliana to agronomically relevant crops, since module-level insights are easier to translate between species than individual genes (Ficklin & Feltus, Reference Ficklin and Feltus2011; Julca et al., Reference Julca, Ferrari, Flores-Tornero, Proost, Lindner, Hackenberg, Steinbachová, Michaelidis, Gomes Pereira, Misra, Kawashima, Borg, Berger, Goldberg, Johnson, Honys, Twell, Sprunck, Dresselhaus and Mutwil2021; Shaik & Ramakrishna, Reference Shaik and Ramakrishna2013; Sreedasyam et al., Reference Sreedasyam, Plott, Hossain, Lovell, Grimwood, Jenkins, Daum, Barry, Carlson, Shu, Phillips, Amirebrahimi, Zane, Wang, Goodstein, Haas, Hiss, Perroud, Jawdy and Schmutz2023). Machine learning (ML) is a promising approach to extract gene modules from expression measurements, as it has demonstrated the ability to infer complex patterns from large amounts of (biological) data (Alber et al., Reference Alber, Buganza Tepole, Cannon, De, Dura-Bernal, Garikipati, Karniadakis, Lytton, Perdikaris, Petzold and Kuhl2019). Nonetheless, ML approaches do not elucidate interactions between biological mechanisms that underlie a plant’s response to environmental influences, and therefore fall short in delivering the detailed understanding that is needed to improve plant stress resilience. Thus, recently there has been an increasing interest in combining mechanistic models with ML approaches (Alber et al., Reference Alber, Buganza Tepole, Cannon, De, Dura-Bernal, Garikipati, Karniadakis, Lytton, Perdikaris, Petzold and Kuhl2019; Baker et al., Reference Baker, Peña, Jayamohan and Jérusalem2018; Noordijk et al., Reference Noordijk, Garcia Gomez, ten Tusscher, de Ridder, van Dijk and Smith2024). In our context, this involves combining ML-derived modules, which extract relevant patterns from large amounts of biological data, with mechanistic models to provide interpretability and causal insights.
For insightful ODE modelling, we want to find modules that are both coherent and interpretable. Here, coherent means that genes within the same module share an expression pattern in a given dataset, allowing us to summarise their expression as a single variable (e.g., the mean) in our ODE model. Interpretable means that they can be linked to a clear biological process, providing insight into how plants respond to stress and deepening our mechanistic understanding. There is a trade-off between these two desired characteristics, which is related to the way in which co-expression is calculated in the context of gene clustering. Firstly, ‘local’ co-expression, calculated from gene expression measured within only one experiment (i.e., over various time points), might yield clusters that are coherent but overfit the data (Mao et al., Reference Mao, Zaslavsky, Hartmann, Sealfon and Chikina2019), resulting in an ill-defined biological function. On the other hand, ‘global’ co-expression is calculated from a gene expression compendium which covers a wide range of different conditions (e.g., developmental stages, tissues or stresses). Global co-expression should reflect gene-gene relationships more generally, providing high biological functional similarity within modules, but potentially overlooking groups of genes that are only co-expressed in a particular context. In BADDADAN, we address this coherence-interpretability trade-off by employing a module finding approach that uses a combination of global and local co-expressions.
Module detection has been extensively researched and reviewed (Saelens et al., Reference Saelens, Cannoodt and Saeys2018), and has been used to infer functionally related groups of genes in various organisms from gene expression data. For example, modules were found by converting gene expression to an association network first, disregarding the temporal dimension (Forster et al., Reference Forster, Li, Yashiroda, Yoshimura, Li, Isuhuaylas, Itto-Nakama, Yamanaka, Ohya, Osada, Wang, Bader and Boone2022; Langfelder & Horvath, Reference Langfelder and Horvath2008; Wang et al., Reference Wang, Dalkic, Wu and Chan2008). Moreover, Janky et al. (Reference Janky, Verfaillie, Imrichová, de Sande, Standaert, Christiaens, Hulselmans, Herten, Sanchez, Potier, Svetlichnyy, Atak, Fiers, Marine and Aerts2014) and Mao et al. (Reference Mao, Zaslavsky, Hartmann, Sealfon and Chikina2019) have extended such methods to incorporate prior biological knowledge. However, none of these methods explicitly leverage temporal data. In contrast, methods that exploit the temporal nature of the data (Han et al., Reference Han, Liu, Zhang, Wang, Bu, Chang, Yu, Li, Tian, Yang, Zhu and He2021; Heard et al., Reference Heard, Holmes, Stephens, Hand and Dimopoulos2005; McDowell et al., Reference McDowell, Manandhar, Vockley, Schmid, Reddy and Engelhardt2018; Schulz et al., Reference Schulz, Devanny, Gitter, Zhong, Ernst and Bar-Joseph2012) do not build upon prior biological data, and suffer from sparse, noisy measurements, common in biology. Additionally, data integration has been employed to find modules (Depuydt & Vandepoele, Reference Depuydt and Vandepoele2021; Forster et al., Reference Forster, Li, Yashiroda, Yoshimura, Li, Isuhuaylas, Itto-Nakama, Yamanaka, Ohya, Osada, Wang, Bader and Boone2022; Heyndrickx & Vandepoele, Reference Heyndrickx and Vandepoele2012; Sanchez-Munoz et al., Reference Sanchez-Munoz, Depaepe, Samalova, Hejatko, Zaplana and Van Der Straeten2024), but these methods cannot easily be applied to new transcriptomics datasets of interest, hampering their ability to improve gene module coherence. Finally, Lu et al. (Reference Lu, Liang, Li and Wu2011), Soltanalizadeh et al. (Reference Soltanalizadeh, Gonzalez Rodriguez, Maroufy, Zheng and Wu2020), and Wu et al. (Reference Wu, Liu, Qiu and Wu2014) performed ODE modelling of gene modules in yeast, mice or human, but did not build upon existing data, potentially reducing the interpretability of their modules and failing to leverage prior knowledge available for these species. Thus, many distinct elements of module finding, biological data integration, and ODE modelling are available, but no method yet integrates all these steps.
Here, we present BADDADAN (a Bioinformatics Approach to Describe Dynamical Activations of a Dimension-reduced A-priori-informed Network) which integrates existing transcriptomics and transcription factor (TF) binding data with experiment-specific expression data to infer a dynamic network model of gene modules that are coherent and reflect clear biological functions. As a proof of concept, we show that we can quantitatively model gene module temporal expression through a system of ODEs with reasonable accuracy in two real-world datasets, and that these models recover known biological insight as well as provide novel hypotheses for drought or heat resilience mechanisms. Taken together, BADDADAN offers a quantitative and interpretable framework for studying temporal gene expression responses to external perturbations, applicable to plants and other organisms.
2. Results
To demonstrate BADDADAN (Figure 1), we apply it to two time-series transcriptomics datasets in A. thaliana Col-0: 1) a progressive drought stress dataset of five-week-old plants grown in increasingly dry soil over a period of two weeks, with transcriptomics of leaf 7 sampled each day (Bechtold et al., Reference Bechtold, Penfold, Jenkins, Legaie, Moore, Lawson, Matthews, Vialet-Chabrand, Baxter, Subramaniam, Hickman, Florance, Sambles, Salmon, Feil, Bowden, Hill, Baker, Lunn and Mullineaux2016) and 2) a dataset where six-week-old plants were subjected to elevated ambient temperature, with transcriptomic samples of rosette leaves collected unevenly over a 24-hour period (Caldana et al., Reference Caldana, Degenkolbe, Cuadros-Inostroza, Klie, Sulpice, Leisse, Steinhauser, Fernie, Willmitzer and Hannah2011). We chose these datasets because they contain transcriptomics measurements over a large number of time points (>14), making them suitable for ODE modelling. Furthermore, they represent distinct stress response pathways, which allows us to test our approach on multiple types of stresses. Below, we will go over each step in our approach and its results on both datasets in more detail.
BADDADAN overview. (a) BADDADAN starts with candidate gene module finding based on a combination of experiment-specific (local) expression and compendium-wide (global) co-expression data (Obayashi et al., Reference Obayashi, Hibara, Kagaya, Aoki and Kinoshita2022). (b) To show this leads to improved coherence and interpretability, we compare these modules to modules created from either local or global data alone. (c) Our pipeline then selects a subset of the modules based on four criteria reflecting suitability for ODE modelling. (d) Next, it connects the modules using TFBS enrichment, and (e) creates an ODE model from this intermodular network which is fit to experimental (local) data and allows for biological insights.
2.1. Module finding
As a first step of gene module finding, we filter our dataset for differentially expressed (DE) genes to reduce the chance of fitting modules to noise. DE filtering is performed using limma (Ritchie et al., Reference Ritchie, Phipson, Wu, Hu, Law, Shi and Smyth2015) (see Methods), and we find 2,363 and 4,444 DE genes (FDR
$\leq $
0.05) for drought and heat, respectively, with which we proceed for clustering. The set of DE genes is subsequently clustered. Hereto, we use either the (local) co-expression between the DE genes in the time-series data or (global) co-expression across a large compendium of datasets as compiled in ATTED (Obayashi et al., Reference Obayashi, Hibara, Kagaya, Aoki and Kinoshita2022), and also combine the local and global distances by summing them to express the similarities between the DE genes (see Methods). To investigate if clusters derived from the combined distance can balance coherence and interpretability, we compare them to clusters derived from either global or local distances alone. To assess coherence, we calculate the percentage of variance in gene expression over time explained by the first principal component for each module. To measure interpretability, we identify the enriched biological process gene ontology (GO) terms for each module (see Methods). Coherence is based on expression, which is used explicitly in the clustering process; in contrast, GO annotations are never explicitly used during the clustering, and thus provides an independent assessment of clustering performance. To prevent bias (e.g., modules of size 1 are 100% coherent), we also ensure that the number of modules and their sizes are approximately equal between different clusterings (see Table S1 in the Supplementary Material, Figure S1 in the Supplementary Material and Methods).
In both datasets, we see a similar pattern (Figure 2): local distance-based clusters are more coherent than global distance-based ones, though fewer are enriched for at least one GO term. We furthermore find that combined distance-based modules balance coherence and interpretability, supporting our hypothesis that combining them can address the trade-off between the two. However, most differences between the use of local, combined or global distances for the drought data are not significant (e.g., combined versus local GO enrichment,
$p=0.07$
, Mann-Whitney U test). This could be due to the smaller number of modules in drought (
$\approx 25$
) compared to heat (
$\approx 45$
), which, in turn, reflects the fewer DE genes in the drought experiment. The total number of GO terms per module follows a trend similar to the fraction of modules with at least one annotated GO term (Figure S3 in the Supplementary Material). We also assess GO term functional relatedness per module using mean pairwise semantic similarity (Wang et al., Reference Wang, Du, Payattakool, Yu and Chen2007) (Figure S4 in the Supplementary Material and Supp. Methods in the Supplementary Material), showing that clustering on combined distances increases interpretability without reducing functional relatedness.
Module finding. (a and b) Per-module coherence, i.e., explained variance of the first principal component of modules in drought and heat stress datasets, respectively. (c and d) Fraction of gene modules with at least one enriched biological process GO term. The x-axis represents the distances used to form modules: global (derived from compendium (Obayashi et al., Reference Obayashi, Hibara, Kagaya, Aoki and Kinoshita2022)), local (based solely on the experiment of interest), combined (a sum of local and global) and random (modules formed by random gene assignment). Error bars indicate the 95% c.i. estimated by bootstrapping. Different letters indicate
$p<0.05$
, based on a two-sided pairwise Mann-Whitney U test between distributions.
2.2. Module selecting
Not all found candidate modules may be relevant to include in our dynamic model. We select modules most suitable for ODE modelling based on four characteristics, expressed as z-scores (see Methods): 1) activity, i.e., a high average expression over all samples; 2) coherence, as this facilitates representation of the expression pattern by a single variable; 3) relevance, i.e., differing in expression between control and treatment conditions; and 4) regulatory relevance, hypothesising that modules with a high number of enriched TF binding sites (TFBSs) are more likely to be involved in stress response GRNs (Heyndrickx et al., Reference Heyndrickx, Van de Velde, Wang, Weigel and Vandepoele2014). There is no strong correlation between these four scores (Figures S5 and S6 in the Supplementary Material), suggesting their complementary nature. We use a combined Stouffer z-score of these four characteristics to proceed with the top 50% of all candidate modules (see Methods), resulting in 12 modules for drought and 26 modules for the heat stress datasets. A breakdown of the scores per module shows that the highest scoring modules score well on coherence, relevance and activity, but not always on the regulatory relevance (Figure S7 in the Supplementary Material).
2.3. Module connecting
To incorporate a set of modules into an ODE model of intermodular regulation, we need to find regulatory interactions, i.e., determine which modules influence which other(s). Since such interactions cannot trivially be inferred from co-expression patterns and methods based on, for example, mutual information (Huynh-Thu & Geurts, Reference Huynh-Thu and Geurts2018) do not always perform well (Aibar et al., Reference Aibar, González-Blas, Moerman, Huynh-Thu, Imrichova, Hulselmans, Rambow, Marine, Geurts, Aerts, van den Oord, Atak, Wouters and Aerts2017; Saelens et al., Reference Saelens, Cannoodt and Saeys2018), we choose to find these interactions based on existing knowledge of TFs and their TFBSs in A. thaliana. To this end, we used TF2Network (Kulkarni et al., Reference Kulkarni, Vaneechoutte, Van de Velde and Vandepoele2018), which is based on extensive prior research (e.g., yeast one-hybrid and ChIP-Seq studies). However, connecting modules solely based on TFBSs can yield false positives, because TFs might not actually bind to a potential binding site due to, e.g., low chromatin accessibility (De Clercq et al., Reference De Clercq, Van de Velde, Luo, Liu, Storme, Van Bel, Pottie, Vaneechoutte, Van Breusegem and Vandepoele2021; Khamis et al., Reference Khamis, Motwalli, Oliva, Jankovic, Medvedeva, Ashoor, Essack, Gao and Bajic2018; Kulkarni et al., Reference Kulkarni, Vaneechoutte, Van de Velde and Vandepoele2018). Thus, we additionally require that a module can only be connected to another module when the TF that has a potential binding site to the target module is part of the source module and that this TF has an absolute correlation
$>0.75$
with the target module expression average (Figure S8 in the Supplementary Material and Methods). The sign of the correlation is then used to infer whether the TF activates or inhibits the target module. Moreover, we only consider a TF when the correlation between its expression and the average expression of its source module is above
$0.3$
(Figure S9 in the Supplementary Material and Methods). Finally, we only keep modules that have at least one ingoing or outgoing connection to other modules. The resulting intermodular networks are shown in Figure 3b and d for the drought and the heat stress network, respectively.
Intermodular networks and ODE fits. (a and c) Best ODE fit for the drought and heat stress experiments, respectively. Error bars indicate the 95% confidence interval of the mean module expression. (b and d) Intermodular network for drought and heat stress dataset, respectively. Arrows indicate activation, ‘T’-shaped ends indicate inhibition. Each edge can represent more than one TF. Module numbering is discontinuous due to the module selection step.
2.4. Model fitting
Next, to evaluate whether BADDADAN is able to capture module expression dynamics, we construct ODEs based on the inferred module network, save them in SBML format (Hucka et al., Reference Hucka, Finney, Sauro, Bolouri, Doyle, Kitano, Arkin, Bornstein, Bray, Cornish-Bowden, Cuellar, Dronov, Gilles, Ginkel, Gor, Goryanin, Hedley, Hodgman, Hofmeyr and Wang2003), and fit these using multi-start optimisation (see Methods). For the heat stress dataset, we explicitly model the circadian clock as input for the modules (see Methods) because the time points cover 24 hours and without the circadian clock as input our model was not able to attain an adequately accurate fit. For both datasets, the optimisation problem is hard; i.e., starting from 5,000 initial parameter sets, we rarely achieved the best objective value (Figure S10 in the Supplementary Material).
Nevertheless, the resulting ODE models match the dynamics of the experimental data, with a mean squared error of 0.15 and 0.05 for drought stress (Figure 3a) and heat stress (Figure 3c), respectively. However, they do not capture some of the early dynamics in the treatment condition (e.g., the peaks in modules 3, 4 and 6 for drought stress or the dip of module 27 for heat stress). In the drought stress experiment, the conditions are expected to be highly similar between control and treatment for the first few days, so what causes these peaks is not a priori clear. The fact that early peaks are not fitted well could be due to the way the external input is modelled. For drought, this is a linear increase from 0 to 1 over two weeks. Consequently, the drought condition does not have a strong signal initially in the model. For the heat stress dataset, the external input into our model is a constant 1, even though the influence of the heat stress on certain modules might change over time (Kappel et al., Reference Kappel, Friedrich, Oberkofler, Jiang, Crawford, Lenhard and Bäurle2023). Hence, further exploring the representation of the external input might enhance these fits. Alternatively, the fit may be improved by including additional module interactions or by adopting a more flexible ODE formulation (Eq. (7) in Methods).
2.5. Mechanistic insights
We end up with ODEs for networks of eight and nine modules that represent 1,203 and 1,335 genes in total for drought and heat stress, respectively. Edges between modules are derived from one TF in most cases (12 edges for both drought and heat), with some edges representing two to six TFs (five and two edges for drought and heat, respectively). While comprehensive investigation of all interactions in these models is beyond the scope of this article, below we select some hypotheses from our models and relate these to earlier work (Datasets S1 and S2 in the Supplementary Material).
2.5.1. Drought
For the drought stress experiment, the original publication (Bechtold et al., Reference Bechtold, Penfold, Jenkins, Legaie, Moore, Lawson, Matthews, Vialet-Chabrand, Baxter, Subramaniam, Hickman, Florance, Sambles, Salmon, Feil, Bowden, Hill, Baker, Lunn and Mullineaux2016) found that many genes only become DE from day eight onwards, as this likely marks a shift from mild to severe drought stress. We see a similar pattern in our data (Figure 3a); many modules behave similarly in drought and control for approximately the first half of the experiment, and only diverge after that. Only two modules (
$25\%$
) in our network could not be functionally annotated using GO terms (see Methods and Dataset S1 in the Supplementary Material), a far lower fraction than found in the original paper (Bechtold et al., Reference Bechtold, Penfold, Jenkins, Legaie, Moore, Lawson, Matthews, Vialet-Chabrand, Baxter, Subramaniam, Hickman, Florance, Sambles, Salmon, Feil, Bowden, Hill, Baker, Lunn and Mullineaux2016) where roughly
$50 \%$
of the modules could not be annotated.
The interpretability of our approach helps generating hypotheses. For example, module 12, the only one with decreased expression levels in drought compared to control, is enriched for GO terms related to organismal-level ion homeostasis, iron ion starvation and nutrient availability, suggesting that in the late phases of drought stress these can play a significant role, which has been found before (Mukarram et al., Reference Mukarram, Choudhary, Kurjak, Petek and Khan2021; Mulet et al., Reference Mulet, Campos and Yenush2020, Reference Mulet, Porcel and Yenush2023). According to BADDADAN, genes related to this homeostasis are inhibited by module 4 through the GBF3 (AT2G46270) TF, a bZIP G-box binding protein that has been found to be induced by water deprivation and was originally found in a finger millet drought transcriptome analysis (Ramegowda et al., Reference Ramegowda, Gill, Sivalingam, Gupta, Gupta, Govind, Nataraja, Pereira, Udayakumar, Mysore and Senthil-Kumar2017). Additionally, overexpression of this gene in A. thaliana has been shown to increase drought resistance (Ramegowda et al., Reference Ramegowda, Gill, Sivalingam, Gupta, Gupta, Govind, Nataraja, Pereira, Udayakumar, Mysore and Senthil-Kumar2017). As a potential mechanism, the authors suggest altered ABA signalling, but BADDADAN raises the possibility that GBF3 expression could be associated with ion homeostasis in drought resistance, which would be interesting to validate experimentally.
Module 4 and module 15 are enriched for terms related to oxidative stress and salicylic acid (SA) production or response, which are known to be involved in drought stress (Ortega et al., Reference Ortega, Celoy, Chacon, Yuan, Xue, Pandey, Drowns, Kvitko and Tsai2024). Analysis of the original paper’s data confirms that SA abundance correlates with the average expression of these modules in our ODE, with Pearson correlations of
$0.61$
and
$0.66$
for modules 4 and 15, respectively (Figure S14 in the Supplementary Material). This demonstrates how modules can be successfully linked to higher-level phenotypic traits, such as plant hormone levels. SA has been shown to induce stomatal closure (Okuma et al., Reference Okuma, Nozawa, Murata and Miura2014), and there seems to be a negative correlation between stomatal conductance (Bechtold et al., Reference Bechtold, Penfold, Jenkins, Legaie, Moore, Lawson, Matthews, Vialet-Chabrand, Baxter, Subramaniam, Hickman, Florance, Sambles, Salmon, Feil, Bowden, Hill, Baker, Lunn and Mullineaux2016) and the expression of modules 4 and 15. Unfortunately, given the absence of the original raw measurement data, this is hard to quantify. Module 4 is also enriched for terms associated with biotic stress, in line with a role for stomatal closure in preventing water loss as well as bacterial entry (Okuma et al., Reference Okuma, Nozawa, Murata and Miura2014), and SA plays important roles in both biotic and abiotic stress (Ortega et al., Reference Ortega, Celoy, Chacon, Yuan, Xue, Pandey, Drowns, Kvitko and Tsai2024). Okuma et al. (Reference Okuma, Nozawa, Murata and Miura2014) notes that stomatal closure is light dependent, which matches our model in that module 4 is regulated by module 2, which contains GO terms associated with light response. Module 4 has the highest closeness centrality of 1 in the network, which suggests it is more likely to play a key role in the overall response mechanism. Finally, module 16 is regulated by both modules 2 and 4 and enriched for starch, glucan and maltose metabolism. This suggests that these metabolic processes could be controlled through, e.g., circadian rhythm and light—all from module 2—but also stress-related processes from module 4, all of which have been reported (Lu et al., Reference Lu, Gehan and Sharkey2005; MacNeill et al., Reference MacNeill, Mehrpouyan, Minow, Patterson, Tetlow and Emes2017; Skryhan et al., Reference Skryhan, Gurrieri, Sparla, Trost and Blennow2018). Knocking out in silico (see Methods) the CBF3 (AT4G25480) TF that connects module 2 to module 4 results in altered expression of, e.g., module 16 (Figure S15 in the Supplementary Material), and in literature, it is indeed reported that changed CBF3 expression can impact sugar metabolism (Gilmour et al., Reference Gilmour, Sebolt, Salazar, Everard and Thomashow2000). All in all, these findings show that BADDADAN recovers known mechanisms and TFs, and can link module expression to other measurements such as plant hormone abundance.
Note that the expression profiles of some modules appear relatively similar (e.g., those of modules 2 and 4). Recall that the modules are created by clustering both the local distances as well as the global distances that came from a compendium of data. The latter reflects that the two modules will have different biological functions. This is indeed true as the GO enrichments of modules 2 and 4 are different; module 2 is related to light response, circadian rhythm and organic hydroxy compound biosynthetic process, whereas module 4 is related to biotic and abiotic stress. This shows the potential benefit of clustering based on a combination of local and global expression measurements.
2.5.2. Heat
The heat stress dataset shows a large variety in up- and down-expression of modules over time, which is reflected in a large number of edges representing inhibition in the intermodular network (six inhibiting edges and eight activating ones versus one inhibiting and 16 activating in drought).
Module 4 displays interesting dynamics, initially lower in heat than control but higher at the last time point. Module 4 is linked to heat and temperature stress responses as well as protein folding. Potentially, module 4 might be involved in longer-term thermotolerance, as found recently for one of its constituent genes (HSFA2) (Pan et al., Reference Pan, Yu, Wei, Qiu, Mao, Liu, Yan, Linghu, Li, Guo and Tang2024). BADDADAN suggests this is potentially through protein-folding-related processes, evidence for which can be found (Wang et al., Reference Wang, Zhu, Tang, Wang, Sun and Deng2024). The expression pattern of module 4 seems to result from the type-2 coherent feed-forward loop (Alon, Reference Alon2007) from module 9 through module 11. Module 9 inhibits module 4 through the LCL1 (AT5G02840) TF, which has been found to be involved in heat shock response (Kidokoro et al., Reference Kidokoro, Konoura, Soma, Suzuki, Miyakawa, Tanokura, Shinozaki and Yamaguchi-Shinozaki2023). Module 11 stimulates modules 4, 5, 14, and 29 through the ABF1 (AT1G49720) TF, which is linked to ABA response as well as heat (Song et al., Reference Song, Kim, Chung and Lim2017), highlighting it as a potential key regulator. Unfortunately, ABA was not measured in this experiment, so we cannot study how its abundance correlates to module expression. Nonetheless, knocking out ABF1 in silico (see Methods) reveals a stronger influence on modules 4 and 29 compared to modules 5 and 14 (Figure S15 in the Supplementary Material). Specifically, the parameters
$\beta _{11,4} \approx 195$
and
$\beta _{11,29} \approx 20$
are much higher than
$\beta _{11,5} \approx 7 \cdot 10^{-4}$
and
$\beta _{11,14} \approx 4 \cdot 10^{-4}$
, highlighting that not all inferred intermodular edges represent equally strong relationships. This underscores the value of the interpretable parameters provided by our approach.
Module 5 is initially downregulated in both conditions, with a slightly higher expression in heat in the first 10 hours, after which we observe a higher expression in the control condition. This module is associated with defence response, stomatal movement and response to biotic (bacterial) stimulus. As we also found in the drought dataset, the biotic stimulus is likely related to the stomatal movement (Okuma et al., Reference Okuma, Nozawa, Murata and Miura2014). Module 5 is also enriched for programmed cell death, glucosinalate catabolism and callose localisation, which have been linked to (a)biotic stress responses (Janda et al., Reference Janda, Lamparová, Zubíková, Burketová, Martinec and Krčková2019; Hopkins et al., Reference Hopkins, van Dam and van Loon2009). Overall, this shows how BADDADAN leverages prior knowledge to identify modules potentially related to combined stresses, and suggests that module 5 is an interesting target for studies on interactions between them.
Finally, Module 29 exhibits a similar response in control and heat during the initial phase of the heat shock, but at later time points is more active in heat. This could again be through a type-2 coherent feed-forward loop (Alon, Reference Alon2007), this time from module 5. Module 29 is related to cellular homeostasis, glucosinolate transport and sulphur compound transport, all related to heat stress (Ihsan et al., Reference Ihsan, Daur, Alghabari, Alzamanan, Rizwan, Ahmad, Waqas and Shafqat2019; Sugio et al., Reference Sugio, Dreos, Aparicio and Maule2009). An increase in catabolism of glucosinalate (which contains sulphur) due to module 5 will lead to a need for transport. Module 5 inhibits module 29 through the LHY1 (AT1G01060) TF, linked to the circadian clock; BADDADAN suggests that it could also be involved in controlling sulphur transport, through module 29. This has been observed in Medicago truncatula (Achom et al., Reference Achom, Roy, Lagunas, Picot, Richards, Bonyadi-Pour, Pardal, Baxter, Richmond, Aschauer, Fletcher, Rowson, Blackwell, Rich-Griffin, Mysore, Wen, Ott, Carré and Gifford2022), and LHY1 has been linked to sulphur in A. thaliana (Van De Mortel et al., Reference Van De Mortel, Schat, Moerland, Van Themaat, Van Der Ent, Blankestijn, Ghandilyan, Tsiatsiani and Aarts2008). As the underlying mechanism is still unknown, this is a potentially interesting direction for future research.
Of the modules that initially display a large difference between the control and heat conditions (4, 9, 25, 27 and 45), only three could be annotated for function (4, 9 and 25). Module 9 has interesting annotations such as fluid transport—which suggests the early phase response of heat is causing the plant to transport more water—and hydrogen peroxide transmembrane transport, but because it has no ingoing connections, we cannot infer what processes or TFs are involved in regulating this. This shows that BADDADAN only recovers a part of the intermodular network and thus does not necessarily provide insight into the full network; however, any model based on a subset of genes will be similarly limited.
3. Discussion
Here, we introduce BADDADAN, a framework for constructing an ODE model of gene modules encompassing more than 1,000 genes in total. We apply it to two real-world datasets covering different stresses and time scales, demonstrating its advancements in gene module expression modelling. First, by combining TFBS enrichment with expression data, it identifies directional, signed relationships between modules. Second, by estimating parameters from data that capture these relationships, it quantifies interaction strengths. Third, by combining experiment-specific expression data with compendium-wide co-expression data, it yields coherent and interpretable modules. Together, these features enable accurate predictions of gene module expression over time and in silico predictions of responses to perturbations, such as knockouts, that would not be apparent from network structure alone. Moreover, BADDADAN uncovers known biological processes and their interactions while also generating novel hypotheses that can inform the design of future experiments.
BADDADAN enables out-of-the-box integration of experiment-specific data with prior experimental data (ATTED (Obayashi et al., Reference Obayashi, Hibara, Kagaya, Aoki and Kinoshita2022)), addressing a limitation of the widely-used WGCNA (Langfelder & Horvath, Reference Langfelder and Horvath2008), which does not directly support such integration and assumes a scale-free network topology. In contrast to recent black-box ML approaches (Baker et al., Reference Baker, Peña, Jayamohan and Jérusalem2018; Sapoval et al., Reference Sapoval, Aghazadeh, Nute, Antunes, Balaji, Baraniuk, Barberan, Dannenfelser, Dun, Edrisi, Elworth, Kille, Kyrillidis, Nakhleh, Wolfe, Yan, Yao and Treangen2022; Shen et al., Reference Shen, Liu, Tu and Tang2021), BADDADAN is composed of intermediate steps that are inspectable; e.g., module selection provides z-scores, connections between modules can be visualised as a directed graph, and the ODE model contains human-interpretable parameters (e.g.,
$\beta $
can be linked to interaction strength). Moreover, existing ODE modelling approaches comparable to BADDADAN (Lu et al., Reference Lu, Liang, Li and Wu2011; Wu et al., Reference Wu, Liu, Qiu and Wu2014) infer connections between modules solely in a data-driven way—similar to LASSO regression—which has been shown to underperform for GRN inference (Saelens et al., Reference Saelens, Cannoodt and Saeys2018). However, as potential future research, such connections could be integrated with our TFBS-based ones, a hybrid approach similar to Aibar et al. (Reference Aibar, González-Blas, Moerman, Huynh-Thu, Imrichova, Hulselmans, Rambow, Marine, Geurts, Aerts, van den Oord, Atak, Wouters and Aerts2017) where connections are inferred if they are important for predicting dynamics and/or underpinned by the presence of a TFBS. Also, incorporation of additional data, e.g., in the form of chromatin accessibility, could further inform finding of regulatory connections (Staut et al., Reference Staut, Manosalva Pérez, Matres Ferrando, Dissanayake and Vandepoele2024).
By addressing essential challenges in gene module network modelling—such as edge directionality determination, activation/inhibition inference and associated weights estimation within an ODE framework—BADDADAN opens up extensive possibilities for further analysis. For example, looking deeper into confidence interval estimation of parameters, applying an identifiability analysis to our model or performing topological sensitivity analysis (Babtie et al., Reference Babtie, Kirk and Stumpf2014), where multiple different network topologies are tested for their ability to describe the data. This will tell us if we have found a network topology which would uniquely describe the resulting behaviour. Related to this, approaches for optimisation of a given network topology by adding or removing a small number of nodes (modules in our case) could be applied (Astola et al., Reference Astola, Stigter, van Dijk, van Daelen and Molenaar2014). Also, modelling choices (e.g., the value of the Hill coefficient, linear environmental effects or additive module interactions) can significantly affect the fitted parameter values—and possibly subsequent predictions of knockout effects—so care should be taken when interpreting them. Experimental validation of BADDADAN’s hypotheses will help support or refine these modelling assumptions.
Settings of each individual step are now optimised separately in the pipeline, but it would be an interesting to research if gains can be made by optimising the pipeline end-to-end through differentiable optimisation methods developed for deep learning (AlQuraishi & Sorger, Reference AlQuraishi and Sorger2021). This would also allow the use of deep-learning-based embeddings for clustering, higher dimensional representations of module expressions and neural ODEs or biologically informed neural networks for ODE fitting. On top of that, the loss function could be extended to make modules adhere to, e.g., rate limits or other existing biological knowledge.
Furthermore, to aid in integration with existing or novel methodologies, BADDADAN builds on existing standards and software. For parametrisation of the ODE, we utilise standard formats and scripts for systems biology models (Schälte et al., Reference Schälte, Fröhlich, Jost, Vanhoefer, Pathirana, Stapor, Lakrisenko, Wang, Raimúndez, Merkt, Schmiester, Städter, Grein, Dudkin, Doresic, Weindl and Hasenauer2023; Schmiester et al., Reference Schmiester, Schälte, Bergmann, Camba, Dudkin, Egert, Fröhlich, Fuhrmann, Hauber, Kemmer, Lakrisenko, Loos, Merkt, Müller, Pathirana, Raimúndez, Refisch, Rosenblatt, Stapor and Weindl2021), and we save our models as SBML format files (Hucka et al., Reference Hucka, Finney, Sauro, Bolouri, Doyle, Kitano, Arkin, Bornstein, Bray, Cornish-Bowden, Cuellar, Dronov, Gilles, Ginkel, Gor, Goryanin, Hedley, Hodgman, Hofmeyr and Wang2003), allowing interaction with the wide ecosystem that is already available for systems biology models. This would allow us, for example, to couple BADDADAN to high-quality ODE models for small GRNs that are already available (Chew et al., Reference Chew, Seaton, Mengin, Flis, Mugford, George, Moulin, Hume, Zeeman, Fitzpatrick, Smith, Stitt and Millar2022; Nieto et al., Reference Nieto, Catalán, Luengo, Legris, López-Salmerón, Davière, Casal, Ares and Prat2022; Valentim et al., Reference Valentim, van Mourik, Posé, Kim, Schmid, van Ham, Busscher, Sanchez-Perez, Molenaar, Angenent, Immink and van Dijk2015). Given that both are composed of ODEs, integration should be possible. For example, we could use a BADDADAN module related to the production of a metabolite—similar to how we linked two modules to SA in the drought dataset—and use the abundance of that metabolite predicted by BADDADAN as input for an existing mechanistic model. Vice versa, mechanistic models that predict metabolite abundance could be used as input (cf.
$c(t)$
or
$u(t)$
) to BADDADAN. Alternatively, a TF in the BADDADAN model could be taken to directly influence expression of a gene in the mechanistic model, improving its prediction accuracy.
BADDADAN may face computational challenges in parameter estimation, similar to those in GRN modelling, when applied to datasets requiring ODE modelling of large numbers of modules. However, this seems unlikely to pose a problem in practice, as we successfully apply BADDADAN here to real-world datasets representing complex biological systems. If needed, the number of modules requiring ODE modelling could be reduced by adjusting the deepsplit parameter or by being more stringent in the module selection step, at the cost of modelling accuracy, coherence and/or biological interpretability.
Both datasets contained microarray measurements, so BADDADAN’s performance on RNA-seq data remains untested. However, existing studies suggest that clustering structures are similar between microarray and RNA-seq data (Sîrbu et al., Reference Sîrbu, Kerr, Crane and Ruskin2012), which would suggest that BADDADAN’s module finding would also work on RNA-seq. However, how the specific characteristics of RNA-seq measurements influence BADDADAN’s ODE fitting of module dynamics is still unknown. Hence, the parameter limits and form of the ODEs might have to be re-evaluated. To further investigate this, we have concrete plans to apply BADDADAN to time-series RNA-seq data in the future.
Finally, in the future, BADDADAN could be applied to datasets of other (non-)plant species and/or stresses and could form the basis for translational science, using fundamental insights into model organism A. thaliana resilience mechanisms to breed relevant crops with higher resilience. Since our model can predict possible effects of knockouts under different stresses, it provides possibilities to study the effects of genetic variation on resilience. Also, a comparative analysis of modules between different individual or combined stresses could give insight into the temporal dynamics of processes important for multiple stresses specifically (Rizhsky et al., Reference Rizhsky, Liang and Mittler2002; Verslues et al., Reference Verslues, Bailey-Serres, Brodersen, Buckley, Conti, Christmann, Dinneny, Grill, Hayes, Heckman, Hsu, Juenger, Mas, Munnik, Nelissen, Sack, Schroeder, Testerink, Tyerman and Wigge2023), which would further support targeted resilient crop development. For application to non-model species, alternatives to TF2Network and ATTED should be found. This could involve inferring regulatory interactions purely from expression data or applying the ATTED methodology to existing datasets available for the species of interest.
Summarising, BADDADAN provides a proof-of-concept approach for inferring gene module-based dynamic models from time-series transcriptomics data, to help unravel the intricate gene regulatory mechanisms that plants employ to deal with stresses. Increasing knowledge in this area might eventually be used to breed more climate resilient crops which can be used to provide food to a growing world population, even under the pressure of climate change.
4. Methods
4.1. Gene expression input data
Progressive drought microarray expression data (Bechtold et al., Reference Bechtold, Penfold, Jenkins, Legaie, Moore, Lawson, Matthews, Vialet-Chabrand, Baxter, Subramaniam, Hickman, Florance, Sambles, Salmon, Feil, Bowden, Hill, Baker, Lunn and Mullineaux2016) was extracted from Gene Expression Omnibus (GEO) (Barrett et al., Reference Barrett, Wilhite, Ledoux, Evangelista, Kim, Tomashevsky, Marshall, Phillippy, Sherman, Holko, Yefanov, Lee, Zhang, Robertson, Serova, Davis and Soboleva2013) accession GSE65046, and probe names were converted to TAIR gene IDs using the Python package GEOParse (v2.0.3) (Gumienny, Reference Gumienny2024). For the heat dataset (Caldana et al., Reference Caldana, Degenkolbe, Cuadros-Inostroza, Klie, Sulpice, Leisse, Steinhauser, Fernie, Willmitzer and Hannah2011), expression data was extracted from ArrayExpress accession E-MTAB-375 and again probe names were mapped to gene identifiers using the annotations available on GEO (Barrett et al., Reference Barrett, Wilhite, Ledoux, Evangelista, Kim, Tomashevsky, Marshall, Phillippy, Sherman, Holko, Yefanov, Lee, Zhang, Robertson, Serova, Davis and Soboleva2013). Both original datasets were already
$\log _2$
-transformed to make noise additive instead of multiplicative and to stabilise variance (Archer et al., Reference Archer, Dumur and Ramakrishnan2004).
4.2. Differential gene expression
We used limma (v3.60.3) (Ritchie et al., Reference Ritchie, Phipson, Wu, Hu, Law, Shi and Smyth2015) in R (v4.4.0) for differential gene expression analysis, as it has been shown to perform well on time-course microarray datasets (Moradzadeh et al., Reference Moradzadeh, Moein, Nickaeen and Gheisari2019). For the drought stress dataset, contrasts were used to test if genes changed expression differently between subsequent time points in the drought group compared to the control group. However, such contrast testing only works if multiple biological replicates are available, which was not the case for the heat stress data. Thus, we ran limma on this dataset by fitting a natural cubic spline with five degrees of freedom through the time points using the R splines library. Next, we used limma to conduct an F-test between the conditions on the five parameters that correspond to interaction, i.e., difference between the control and heat groups. Finally, for both drought and control, only genes with a Benjamini-Hochberg false discovery rate (FDR)
$\leq 0.05$
were kept for further processing.
4.3. Module inference
To obtain global (compendium-wide) co-expression patterns, we downloaded pairwise z-scores from ATTED v11.1 (Obayashi et al., Reference Obayashi, Hibara, Kagaya, Aoki and Kinoshita2022), which were calculated based on 27,427 samples and corrected for biases due to differences in sampling conditions and sample sizes (Obayashi et al., Reference Obayashi, Hibara, Kagaya, Aoki and Kinoshita2022). The z-scores (denoted with
$W_{\textrm {global}}$
) were converted to a distance matrix
$D_{\textrm {global}}$
, which represents pairwise distances between genes, with each row and column corresponding to a gene:
$$ \begin{align} D_{\textrm{global}} = \max (W_{\textrm{global}}) - W_{\textrm{global}}, \end{align} $$
where
$\max (M) = \max _{i,j} M_{ij}$
is the maximum value over all rows and columns. This is consistent throughout the manuscript and also holds for uses of
$\min $
.
Local distances were calculated within each dataset as follows:
$$ \begin{align} D_{\textrm{local}} = 1 - C, \end{align} $$
where C is the matrix of pairwise Pearson correlations between genes over the samples (taking the average expression level if biological replicates were available). Saelens et al. (Reference Saelens, Cannoodt and Saeys2018) found that this ‘unsigned’ (Langfelder et al., Reference Langfelder, Luo, Oldham and Horvath2011) metric yielded higher module detection accuracy compared to, e.g., an absolute Pearson correlation or Spearman correlation.
To obtain a clustering based on both global and local co-expression, we combined
$D_{\textrm {global}}$
and
$D_{\textrm {local}}$
. To avoid biases due to a difference in distributions, we first separately converted both distances to a z-score
$$ \begin{align} Z = \frac{D - \mu_D} {\sigma_D}, \end{align} $$
where
$\mu _D$
is the average and
$\sigma _D$
is the standard deviation of all distances.
Next, we summed the resulting z-score matrices:
$$ \begin{align} Z_{\textrm{combined}} = w \cdot Z_{\textrm{global}} + (1-w) \cdot Z_{\textrm{local}}, \end{align} $$
where w influences the balance between global and local co-expression; i.e., a lower w would result in higher coherence if that is more desired (Figure S2 in the Supplementary Material). We set
$w=0.5$
for all experiments as we want an equal balance between interpretability and coherence.
Finally, to ensure no distances were negative (needed for clustering), we converted them to a positive distance:
$$ \begin{align} D_{\textrm{combined}} = Z_{\textrm{combined}} - \min (Z_{\textrm{combined}}). \end{align} $$
To find modules, we use hierarchical clustering. This has been shown to obtain reasonable performance in gene module detection (Saelens et al., Reference Saelens, Cannoodt and Saeys2018), where it was benchmarked on simulated and experimental data with a known ground-truth. Average linkage hierarchical clustering was performed in R on
$D_{\textrm {combined}}$
,
$D_{\textrm {global}}$
and
$D_{\textrm {local}}$
. To extract modules from the resulting dendrogram, it could be cut at a specific height, but such a global approach has been shown to be suboptimal for gene module detection (Langfelder et al., Reference Langfelder, Zhang and Horvath2008). Thus, we applied a Dynamic Tree Cut method (Langfelder et al., Reference Langfelder, Zhang and Horvath2008) designed specifically for gene module detection. Also, since module size affects both coherence (e.g., modules of size 1 are always 100% coherent) and interpretability (e.g., a module containing all genes would lack enriched GO terms), we ensured that module sizes were roughly equal between the three clusterings. To achieve this, we tuned the deepsplit parameter (Figure S1 in the Supplementary Material). Increasing deepsplit changes the internal dynamic tree cut thresholds used to define clusters, promoting splits that yield lower within-cluster and higher between-cluster dissimilarity (Langfelder et al., Reference Langfelder, Zhang and Horvath2008). This results in more, but smaller, clusters. For the drought stress data, we set deepsplit=1 (default value) for the global and combined distance matrices, but used deepsplit=2 for the local one, as this ensured cluster sizes similar to the other two clusterings (Figure S1a in the Supplementary Material). For the heat dataset, module sizes were comparable for all input distances at deepsplit=1 (default value) (Figure S1b in the Supplementary Material). All other Dynamic Tree Cut settings were kept at their default values.
4.4. Characteristics of clusterings
We compared our local, combined and global distance-based clusterings on coherence and interpretability (GO enrichment) of the resulting modules.
4.4.1. Coherence
To assess coherence, we performed a PCA on each module separately, where we first normalised the expression of each gene to have mean 0 and variance 1 (Eq. 3). Coherence was then calculated as the percentage of variance explained by the first principal component, found using the PCA function in Scikit-learn (v1.1.3) (Pedregosa et al., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion, Grisel, Blondel, Prettenhofer, Weiss, Dubourg, Vanderplas, Passos, Cournapeau, Brucher, Perrot and Duchesnay2011). We plotted the distribution of these coherence scores over all modules using Seaborn (v0.11) (Waskom, Reference Waskom2021).
4.4.2. GO enrichment
GO enrichment was performed using GOAtools (v1.2.3) (Klopfenstein et al., Reference Klopfenstein, Zhang, Pedersen, Ramírez, Warwick Vesztrocy, Naldi, Mungall, Yunes, Botvinnik, Weigel, Dampier, Dessimoz, Flick and Tang2018) with all DE genes as background. Only gene annotations of experimental evidence codes were used to prevent data leakage from existing co-expression-based electronic annotations. We focused on biological process GO terms, as we were interested in biological interpretation of the gene modules. For all modules, only enriched GO terms with a Benjamini-Hochberg FDR
$\leq 0.05$
were used for further analysis. To study the variability in the fraction of modules that had at least one enriched GO term, we used Seaborn’s bootstrapping (Waskom, Reference Waskom2021) to estimate the 95% confidence interval.
4.5. Statistical testing
Statistical tests were performed using the Python package statannotations v(0.5.0) (Charlier et al., Reference Charlier, Weber, Izak, Harkin, Magnus, Lalli, Fresnais, Chan, Markov, Amsalem, Proost, Krasoulis and Getzze Repplinger2022). Because we could not assume all data to be normally distributed, we compared coherence and interpretability using the non-parametric two-sided Mann-Whitney U test. As a null model, we measured the characteristics of random modules, where we ensured that the module sizes were distributed identically to those of the combined distances by randomly permuting these cluster assignments.
4.6. Module selection
To select candidate modules for modelling, we assigned each module a score composed of four different components. 1) Activity: to obtain modules that are composed of genes that are active, we determined the average expression over all samples of all genes in a module. 2) Coherence: to measure the coherence of a module, we performed PCA on the normalised expression of its constituent genes over time in our dataset of interest and measured the percentage of variance explained by the first principal component, as described above. Coherence is desired because it facilitates representation of the module expression by a single value (e.g., average). 3) Relevance: to assess whether a module’s expression differs between treatment and control, we first summarised the module by calculating the average expression of its genes over time for each condition. We then measured the difference between conditions using the mean squared error of these averaged expression profiles. 4) Regulatory relevance: we assumed modules with a high number of enriched TFBSs were more likely to be important for a plant’s stress response; e.g., genes related to stimulus response and gene regulation in A. thaliana are known to have a higher number of regulators (Heyndrickx et al., Reference Heyndrickx, Van de Velde, Wang, Weigel and Vandepoele2014). Moreover, modules with many enriched TFBS are more likely to successfully be linked to other modules in our model. To score this regulatory relevance, we used the count of TFBS
$N_T$
that were enriched in a module according to TF2Network (Kulkarni et al., Reference Kulkarni, Vaneechoutte, Van de Velde and Vandepoele2018). Since this distribution was highly skewed, we applied a transformation
$\log (N_T + 1)$
. We chose not to include interpretability (i.e., GO enrichment) in the scoring to ensure that it remains an independent criterion for evaluating the selected modules.
To ensure no single criterion dominated the selection, all scores were scaled to z-scores through Eq. (3). Subsequently, we combined these individual z-scores using Stouffer’s method to obtain one combined z-score per module because this allowed for subsequent statistical interpretation and filtering:
$$ \begin{align} Z = \frac{\sum_{i=1}^k Z_i}{\sqrt{k}}. \end{align} $$
Here, k is the number of different characteristics (in our case, four). Only the activity was based on non-normalised gene expression, the other scores were calculated from gene expression where each gene was z-scored over all samples. We selected the top 50% (a user parameter, here chosen arbitrarily) of all modules based on their summary score for subsequent modelling. The module selection step results in an intermodular network that has a higher density compared to the network from modules that have not been filtered (Table S3 in the Supplementary Material).
4.7. Model structure
To infer connections between modules, we used the TFBS enrichment results of TF2Network (Kulkarni et al., Reference Kulkarni, Vaneechoutte, Van de Velde and Vandepoele2018). We inferred regulatory interactions between modules based on three criteria: 1) The target module contained an enriched binding site for a TF that belonged to a source module. 2) The TF had a Pearson correlation of at least 0.3 with the average expression of its source module; any TFs with correlations to their source module lower than this were considered an outlier and excluded from regulatory interaction inference. 3) The absolute Pearson correlation between the expression of the TF and the average expression of its target module had to exceed a certain threshold. This step fulfilled two functions. Firstly, it reduced the number of false positives that could result from TFBS enrichment analysis, e.g., due to chromatin accessibility (De Clercq et al., Reference De Clercq, Van de Velde, Luo, Liu, Storme, Van Bel, Pottie, Vaneechoutte, Van Breusegem and Vandepoele2021; Khamis et al., Reference Khamis, Motwalli, Oliva, Jankovic, Medvedeva, Ashoor, Essack, Gao and Bajic2018; Kulkarni et al., Reference Kulkarni, Vaneechoutte, Van de Velde and Vandepoele2018). Secondly, as we needed to choose a different form of the Hill equation for inhibition of activation (Eq. 7), we used the sign of the Pearson correlation as a heuristic to infer up- or down-regulation. If multiple TFs from the same source module regulated the same target module, we asserted they all consistently up- or down-regulated; if there was disagreement, no valid network would be output by BADDADAN. Finally, modules that were left with no interactions with any other modules were removed from the model.
To determine the correlation cutoff between TFs and their target modules, we inferred model structures as described above for cutoffs varying from 0 to 0.95 at intervals of 0.05, and kept track of the total number of modules, interactions and whether the model graph consisted of one single weakly connected component or multiple disjoint components (Figure S8 in the Supplementary Material). We then picked 0.75 as a cutoff for both the heat and drought stress datasets, as this led to the sparsest model not split into disjoint components. Intermodular network visualisation and analysis were performed in NetworkX (v2.8.8) (Hagberg et al., Reference Hagberg, Schult, Swart, Varoquaux, Vaught and Millman2008) and Cytoscape (v3.10.3) (Shannon et al., Reference Shannon, Markiel, Ozier, Baliga, Wang, Ramage, Amin, Schwikowski and Ideker2003).
4.8. ODE model
For each module, we model its activity in an ODE representing the average gene expression levels over time after normalisation (Eq. 3) because related genes can have similar expression patterns but at different magnitudes (Wu et al., Reference Wu, Liu, Qiu and Wu2014). Since z-score normalisation of the gene expression could result in negative values, a constant of 3 (minimum value of module expression was around −2.5) was added to guarantee positive values, preventing biologically implausible scenarios where negative expressions would invert the decay term and falsely suggest gene production. For a module m that is inhibited by a set of modules I and activated by a set of modules A, we represent its rate of change as follows:
$$ \begin{align} \frac{dy_m}{dt} = \gamma_m \cdot u(t) - \delta_m \cdot y_m + c_m(t) + \sum_{i \in I} \beta_{i,m} \frac{k_{i,m}^2}{k_{i,m}^2 + y_i^2} \notag \\ + \sum_{a \in A} \beta_{a,m} \frac{y_a^2}{k_{a,m}^2 + y_a^2}. \end{align} $$
The expression of each module
$y_m$
is dictated by an external stimulus-dependent term weighted by
$\gamma _m$
, a decay term governed by the parameter
$\delta _m$
, a circadian clock input
$c_m(t)$
that we only used for the heat stress dataset (because this is sampled over a period of 24 hours) and activation/inhibition Hill equations of the second order where
$\beta _{n,m}$
represents the maximum magnitude of the activation/inhibition of module n on module m, and
$k_{n,m}$
the abundance of the regulator at which the module is 50% activated/inhibited. We chose these formulas as gene regulation is highly nonlinear and Hill kinetics are widely used in ODE modelling in systems biology (Frank, Reference Frank2013; Polynikis et al., Reference Polynikis, Hogan and di Bernardo2009; Valentim et al., Reference Valentim, van Mourik, Posé, Kim, Schmid, van Ham, Busscher, Sanchez-Perez, Molenaar, Angenent, Immink and van Dijk2015). To investigate the influence of the Hill coefficients, i.e., the strength of the nonlinearity, we also performed parameter optimisation (see below) with Hill coefficients of 1 and 3, and found that the ability to fit the data was not lost when the strength of nonlinearity was changed (Figure S11 in the Supplementary Material). However, resulting parameter values differed significantly (Figures S12 and S13 in the Supplementary Material), suggesting that the choice of Hill coefficient alters BADDADAN’s prediction of underlying dynamics.
The function
$u(t)$
captures the external input, which was different for the drought and heat experiment. For the former, we modelled the linear increase in drought stress over two weeks (Bechtold et al., Reference Bechtold, Penfold, Jenkins, Legaie, Moore, Lawson, Matthews, Vialet-Chabrand, Baxter, Subramaniam, Hickman, Florance, Sambles, Salmon, Feil, Bowden, Hill, Baker, Lunn and Mullineaux2016) as follows:
$$ \begin{align} u(t) = \begin{cases} \frac{t}{13}, & 0 \leq t \leq 13, \quad \text{if drought}, \\ 0, & \text{control}, \end{cases} \end{align} $$
to simulate a linear increase of drought where t is time in days.
For the latter, heat stress, we modelled the external input as follows:
$$ \begin{align} u(t) = \begin{cases} 1, & 0 \leq t \leq 24, \quad \text{if heat}, \\ 0, & \text{control}, \end{cases} \end{align} $$
as is also common in plant mechanistic modelling papers (Nieto et al., Reference Nieto, Catalán, Luengo, Legris, López-Salmerón, Davière, Casal, Ares and Prat2022). Here, we measured t in hours.
Initially, the model struggled to accurately capture the data in the heat experiment, likely due to the circadian clock having a strong influence on gene expression over a 24-hour period (which is not observed during daily sampling in the drought dataset). To address this, we explicitly included a circadian clock input
$c_m(t)$
for each module:
$$ \begin{align} c_m(t) = \begin{cases} a_m \sin\left(\frac{t}{24} 2\pi + \phi_m \right) + b_m, & 0 \leq t \leq 24, \quad \text{if heat}, \\ 0, & \text{if drought}. \end{cases} \end{align} $$
Here,
$a_m$
quantifies the input strength of the circadian clock on the module expression change,
$\phi _m$
determines the phase offset of the circadian clock influence (which we give a 24-hour period by definition) and
$b_m$
models the baseline influence of the circadian clock. Again, t is measured in hours. Finally, BADDADAN exported the full ODE models in the Systems Biology Markup Language (Hucka et al., Reference Hucka, Finney, Sauro, Bolouri, Doyle, Kitano, Arkin, Bornstein, Bray, Cornish-Bowden, Cuellar, Dronov, Gilles, Ginkel, Gor, Goryanin, Hedley, Hodgman, Hofmeyr and Wang2003) to ensure reproducibility, consistency and interoperability.
4.9. Parameter optimisation
We specified the parameter optimisation problem per dataset (drought or heat) in PETAB (Schmiester et al., Reference Schmiester, Schälte, Bergmann, Camba, Dudkin, Egert, Fröhlich, Fuhrmann, Hauber, Kemmer, Lakrisenko, Loos, Merkt, Müller, Pathirana, Raimúndez, Refisch, Rosenblatt, Stapor and Weindl2021) and optimised the parameters of the ODE model to match experimental measurements on both the control and treatment conditions using PyPESTO (Schälte et al., Reference Schälte, Fröhlich, Jost, Vanhoefer, Pathirana, Stapor, Lakrisenko, Wang, Raimúndez, Merkt, Schmiester, Städter, Grein, Dudkin, Doresic, Weindl and Hasenauer2023) and AMICI (Fröhlich et al., Reference Fröhlich, Weindl, Schälte, Pathirana, Paszkowski, Lines, Stapor and Hasenauer2021) for ODE integration. We started the ODE model with
$t=0$
being the measured initial module expressions and predicted the full time course data from there. Both conditions shared all parameters except for
$u(t)$
and
$c_m(t)$
. We minimised the negative log-likelihood objective function using the L-BFGS optimisation algorithm implemented in SciPy (Virtanen et al., Reference Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser, Bright, van der Walt, Brett, Wilson, Millman, Mayorov, Nelson, Jones, Kern, Larson and van Mulbregt2020) and AMICI’s adjoint sensitivity method. For fitting, PyPESTO required the noise level for each module to calculate the negative log-likelihood. We set this value to 0.1, based on the 95% confidence interval estimated through bootstrapping of the mean per module. Parameter limits were chosen empirically to ensure the parameters of a best fit never hit their upper or lower bounds (Table S2 in the Supplementary Material), and all non-negative parameters were estimated on a
$\log _{10}$
scale to aid exploring different orders of magnitude (Schälte et al., Reference Schälte, Fröhlich, Jost, Vanhoefer, Pathirana, Stapor, Lakrisenko, Wang, Raimúndez, Merkt, Schmiester, Städter, Grein, Dudkin, Doresic, Weindl and Hasenauer2023). To maximally explore the optimisation landscape, models with random initial parameters were initialised 5,000 times in parallel, and waterfall plots were created to inspect if a global minimum was found (Figure S10 in the Supplementary Material). In the heat dataset, the last two data points (
$t \approx 10$
and
$t \approx 20$
) were oversampled fourfold to mitigate the effects of uneven sampling and ensure a more balanced representation across time points. Finally, to predict TF knockout effects in silico, we set the parameter
$\beta $
of the associated edge(s) to
$10^{-6}$
. ODE fits and resulting parameters were stored using MLFLOW (Chen et al., Reference Chen, Chow, Davidson, DCunha, Ghodsi, Hong, Konwinski, Mewald, Murching, Nykodym, Ogilvie, Parkhe, Singh, Xie, Zaharia, Zang, Zheng and Zumar2020). Optimisation took approximately four days per dataset on a server equipped with two Intel(R) Xeon(R) E5-2640 v3 CPUs.
Open peer review
To view the open peer review materials for this article, please visit http://doi.org/10.1017/qpb.2025.10016.
Acknowledgements
The authors thank Tijmen den Butselaar for helping find the public data that was used to test BADDADAN, Elena Del Pup and Francesca Giaume for helpful discussions regarding the visualisation and biological interpretation of the results and Jordi Alonso Esteve for input on the methodology.
Competing interests
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Data availability statement
All source code is available at https://github.com/CropXR/BADDADAN.
Author contributions
B.N., A.D.J.vD. and D.dR. conceived the initial idea. B.N. carried out the research and wrote the initial manuscript. M.R., A.D.J.vD. and D.dR. provided critical feedback. All authors have read and approved the final manuscript.
Funding statement
This publication is part of the long-term program PlantXR: A new generation of breeding tools for extra-resilient crops (KICH3.LTP.20.005), which is financed by the Dutch Research Council (NWO), the Foundation for Food & Agriculture Research (FFAR), companies in the plant breeding and processing industry, and Dutch universities. These parties collaborate in the CropXR Institute (www.cropxr.org) that is founded through the National Growth Fund (NGF) of The Netherlands.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/qpb.2025.10017.
Open access
