Dynamic modeling of transcriptional gene regulatory network uncovers distinct pathways during the onset of Arabidopsis leaf senescence

Age-dependent senescence is a multifaceted and highly coordinated developmental phase in the life of plants that is manifested with genetic, biochemical and phenotypic continuum. Thus, elucidating the dynamic network modeling and simulation of molecular events, in particular gene regulatory network during the onset of senescence is essential. Here, we constructed a computational pipeline that integrates senescence-related co-expression networks with transcription factor (TF)-promoter relationships and microRNA (miR)-target interactions. Network structural and functional analyses revealed important nodes within each module of these co-expression networks. Subsequently, we inferred significant dynamic transcriptional regulatory models in leaf senescence using time-course gene expression datasets. Dynamic simulations and predictive network perturbation analyses followed by experimental dataset illustrated the kinetic relationships among TFs and their downstream targets. In conclusion, our network science framework discovers cohorts of TFs and their paths with previously unrecognized roles in leaf senescence and provides a comprehensive landscape of dynamic transcriptional circuitry.


Gene expression analyses reveals important emergent DEGs
The primary data source of this study is a reproducible high-resolution Hiseq2000 time course transcriptomic expression data generated by Woo et al. 1 , retrieved from NCBI's Gene Expression Omnibus 2 with the accession number GSE43616. The FPKM dataset was pre-processed, log2 transformed and fold changed with respect to day 4. We performed differentially expressed genes (DEGs) analysis as described by Woo et al. 2016. These log2 transformed expression profiles generated with expression parameters; >1 for up-regulated genes and <-1 for down-regulated genes. This produces a total of 11,453 differentially expressed genes (DEGs) for senescence. Cumulative and emerging DEGs are identified at each time point. We detected that cumulatively 475 genes are differentially expressed irrespective of time course from day 6 to day 30 ( Supplementary Fig. 1a; Supplementary Table S1). While, several DEGs are emerging at each time point ranging from 475 DEGs to 1,317 DEGs at day 6 and day 30, respectively. We observed that emergent DEGs count follow a steep increase from day 6 to day 14, and day 16 to day 30 in two distinct patterns. However, there is a decrease in number of emergent DEGs at day 16 (291 DEGs) compared to day 14 (1,217 DEGs) and day 18 (618 DEGs) (Supplementary Fig. 1a; Supplementary Table  S1). For antisense transcripts, we found approximately 600 DEGs during each sampling time point irrespective of time course from day 6 to day 30, We observed that emergent genes denoting the number of DEGs decreasing from day 6 until day 18 and begin to increase with maximum number of emergent genes at day 30 ( Supplementary Fig. 5a; Supplementary Table S1). This uncertain decrease in emergent DEGs at day 16 might be due to transition of developmental stage from growth to maturation (G to M) and maturation to senescence (M to S). Importantly, we discovered that at least 4.17% of the emerging DEGs are present in the cumulative list of DEGs at any time point ( Supplementary Fig. S1a, Supplementary Fig. 5a, and Supplementary Table S1).

Weighted gene co-expression network construction identified significant modules within two different senescence networks
The availability of high-resolution, large-scale transcriptome datasets enables coexpression network analysis to identify clusters of highly correlated genes which are potentially co-regulated in response to a physiological condition within a cell 3 . Thus, coexpression networks have huge potential for identification a set of genes, which might participate in a common biological process. To determine senescence associated common gene signatures, we implemented a correlation-based R package; weighted gene co-expression network analysis (WGCNA) 4,5 on two subsets of dataset with 11,453 DEGs representing overall senescence process and 2,226 DEGs representing expression of antisense transcripts. Specified complexity of the multi-course dataset, connections to the particular network (Supplementary Table S2), while information centrality is a substitute of closeness centrality centered on definite resistance amid nodes in a specified network to determine the course of information starting one to other node built on length of largest subnetwork of a network 8 . It is found that Blue module genes have significantly higher degree than turquoise module genes ( Supplementary  Fig. S5b).
Additionally, we computed information centrality with weighted parameter (i.e. all edges are treated as per their actual weights). The information centrality IC(i) for node i in a graph G is calculated as: Here, n is the total count of nodes and = ( + − ) − 1, is a component of R matrix. D is a weighted degree diagonal matrix for each node and J is a matrix consisting of 1 for all elements. Therefore, = ( ) = [ − + ] − 1. Mathematically, is well-defined as infinite. Hence, 1 = 0.
The information centrality distribution of ASCENT reveals that the blue module has significantly higher information centrality compared to turquoise and ASCENT ( Supplementary Fig. S7a). Furthermore, we also calculated average information centrality of the largest component of the network. It is observed that blue module has highest average IC (0.0069) followed by 0.0027 for turquoise module, while average IC for ASCENT was only 0.0009 ( Supplementary Fig. S7b). Thus, our data advocates that antisense Arabidopsis leaf senescence transcripts relies on collective genes of blue and turquoise modules to spread the information through DEGs. Also, the shortest paths between blue module and turquoise module genes were calculated to determine total numbers of shortest paths navigate through a node. Shortest path distribution reveals that genes in the blue module are closer to each other than turquoise and ASCENT (F.test, P= 0.020970693 between Blue and ASCENT, F.test, P= 0.295933097 between turquoise and ASCENT) ( Supplementary Fig. S7c). Thus, any stimuli can travel faster in blue module than turquoise module.
Similarly, the information centrality distribution of SCENT reveals that the darkolivegreen, greenyellow, violet, red, grey60, brown, blue modules have significantly higher information centrality compared to turquoise and SCENT (Fig. 1c). Furthermore, we also calculated average information centrality of the largest component of the network. It is observed that darkolivegreen module has significant highest average IC (0.0022) followed by 0.0018 for greenyellow and 0.0014 for violet modules while average IC for SCENT was only 0.00049 ( Supplementary Fig. S3a). Thus, our data suggests that Arabidopsis leaf senescence relies on collective genes of darkolivegreen, greenyellow, violet, red, grey60, brown and blue modules to spread the information throughout nuclear DEGs. Also, the shortest paths between all seven modules genes were calculated to fine the total count of shortest paths navigate through a node. Shortest path distribution reveals that genes in the darkolivegreen, violet, greenyellow and red modules are more close to each other than grey60, brown, blue, turquoise and SCENT (F test; p<9.8e -10 between darkolivegreen and SCENT, p= 2.18e -08 between violet and SCENT, p= 5.136e -08 between greenyellow and SCENT, p= 1.735e -06 between grey60 and SCENT, p= 1.025e -05 between red and SCENT, p= 0.0001345 between brown and SCENT, p= 0.0001265 between blue and SCENT) ( Supplementary Fig. S3b, Supplementary  Table S2). Thus, any stimuli can travel faster in darkolivegreen module than greenyellow module.

Module significance by gene ontology term enrichment
Gene ontology (GO) terms enrichment analyses within blue and turquoise modules of ASCENT were performed with DAVID 9 and Panther 10 to offer a genetic explanation of the constructed gene co-expression network. Both modules were augmented in GO terms associated to explicit Biological Processes (BP). Significance validation of the detected GO term enrichments for both modules was supported by Panther. The blue module having 356 genes with the largest number of enriched terms included enriched GO categories associated with transcriptional regulation, oxidation-reduction process, signal transduction, response to abscisic acid and response to salt stress GO:0006355, GO:0055114, GO:0007165, GO:0009737 and GO:0009651 respectively ( Supplementary Fig. S8a). While, the turquoise module, comprised of 632 genes obligated an over-representation of BP terms related to transcriptional regulation, oxidation-reduction process, photosynthesis, defense response, response to cold and response to ABA (e.g., GO:0006355, GO:0055114, GO:0015979, GO:0006952, GO:0009409 and GO:0009737) ( Supplementary Fig. S8b, Supplementary Table S3). Additionally, we identified the protein classes encoded by genes in both modules by Panther. It is observed that Nucleic Acid Binding, Oxidoreductase, Hydrolase, Transferase, Ligase and Transporter protein classes were coded abundantly with minimum ten genes in blue module. While, prominent proteins classes encoded by turquoise module genes are Oxidoreductase, Nucleic Acid Binding, Transferase, Hydrolase, Ligase and Transporters (Supplementary Table S3).
Similarly, we pulled out the genes of seven significant modules from SCENT for GO analysis. We formed two module groups; 1 st (darkolivegreen, greenyellow, violet) and 2 nd (red, grey60, brown and blue). It was observed that genes of the 1 st group are mostly involved in oxidation-reduction process (41 genes) followed by protein ubiquitination (25 genes), protein transport and multicellular organism development (17 genes each) and vesicle-mediated transport (11 genes) ( Supplementary Fig. S4a, Supplementary Table S3). While, the 2 nd group of module genes are involved in protein phosphorylation (101 genes) followed by response to salt stress (71 genes), response to cold (42 genes), protein folding (39 genes), response to water deprivation (38 genes) and abscisic acid-activated signaling pathway (32 genes) ( Supplementary Fig. S4b, Supplementary Table S3). Also, panther protein classification signifies hydrolase being the most encoded protein by both module groups followed by transferase, nucleic acid binding, oxidoreductase, transporter, transcription factor, enzyme modulator, cytoskeletal protein, ligase and lyase. Importantly, it must be noted that DAVID and Panther GO term enrichment classification represents almost similar classification.

A. Transcription factor associated regulatory network
Direct coordination of major switches in gene expression relies on transcription factors (TFs) and their downstream target genes in an intricate fashion. Thus, TFs can play a dominating role in controlling the onset and complete development of the senescence in leaves. TFs bind to particular target genes' regulatory sequences to activate and/or suppress their transcription. Importantly, several senescence associated genes encoding transcription factors have been known; NAC, MYB, AP2/EREBP, WRKY, C2H2 zinc-finger, GRAS, and bZIP being the largest groups of TFs family [11][12][13] . However, a systems-level study that discovers a comprehensive list of TFs in leaf senescence is still missing. To investigate the TF-targets coordination of 11,453 nuclear genes and 2,226 antisense transcripts across the Arabidopsis leaf lifespan, we examined the significance of TFs in transcriptional regulation. TFs, we performed a massive TF binding regulatory sequence scanning on 11,453 senescence nuclear genes ( Supplementary Fig. S9). First, we generated all possible combinations (forward, reverse, complement, reverse complement along with degeneracy) of nucleotides for binding of respective TF family based on text mining and family associated TF motif search. Second, we retrieved the 1,000 upstream, 5' UTR sequences and AGI transcript of 11,453 genes. Third, 111 TFs regulatory sequence were subjected to blastn 17 against 1,000 upstream, 5' UTR and AGI transcript sequences of 11,453 genes. Finally, 458,302 possible TF/target combinations (Supplementary Table S4) were identified in Arabidopsis leaf lifespan.
Similarly, among 2,226 antisense transcripts, 163 TFs are found among all annotated TFs in Arabidopsis. To determine the TFs binding to the senescence target genes on regulatory sequences we implied the respective TFs to AGRIS 14 , AthaMap 15 , and Plant Cistrome 16 databases. Out of 163 TFs, 65 TFs and their corresponding target genes were available through aforementioned online resources, while 98 senescence nuclear genes TFs and their targets were not reported on any database. To determine the TF/target pairs of remaining 98 TFs, we performed a massive TF binding regulatory sequence scanning on 2,226 senescence antisense transcripts ( Supplementary Fig. S9). First, we generated all possible combinations (forward, reverse, complement, reverse complement along with degeneracy) of nucleotides for binding of respective TF family based on text mining. Second, we retrieved the 1,000 upstream, 5' UTR sequences and AGI transcript of 2,226 antisense transcripts. Third, 98 TFs regulatory sequence were subjected to blastn 17 against 1,000 upstream, 5' UTR sequences and AGI transcript of 2,226 genes. Finally, 204,761 possible TF/target combinations were identified for 163 TFs in 2,226 senescence antisense transcripts (Supplementary Table S4).
We combined co-expression networks (SCENT and ASCENT) and integrated TF/target association data to represent overall TF/target regulation in senescence. Interestingly, we report 220 TFs in Arabidopsis leaf senescence. While, Woo et al. 1 reported only 94 TFs in senescence, we captured more than double TFs of which 190 TFs are identified by our network biology framework (Fig. 2a, Supplementary Table S4). Moreover, our pipeline could be used to simulate any TF-target relationships in leaf senescence since we are providing all the raw data in the supplementary information. Finally, this computational pipeline can be used for other plant-related processes.
Furthermore, we focused our TF/target identification to blue and turquoise modules of ASCENT in order to find the TF affecting network in Arabidopsis leaf senescence. It has been observed that 24 TFs and 48 TFs were present in blue and turquoise modules respectively (Supplementary Table S4). The expression pattern of these TFs in senescence data revealed that only nine TFs are down-regulated at early 'G to M' phase, eight TFs down regulated at transition phase ('G to M' to 'M to S'), while 12 TFs are up-regulated and four TFs are down-regulated at late 'M to S' phase in blue module. Interestingly, 12 TFs are part of hub 75 genes pertaining high connections. Also, TFs of turquoise module followed a similar pattern; 24 TFs are down-regulated at early 'G to M' phase, 18 TFs down regulated and three TFs are up-regulated at transition phase ('G to M' to 'M to S'), while 11 TFs are up-regulated and 18 TFs are down-regulated at late 'M to S' phase in senescence life span. While, two TFs are part of hub 75 genes pertaining high connections. Interestingly, significant greater quantity of differentially expressed TFs at 'M-to-S' phase advocates that Arabidopsis leaves employ the alterations in TFs more dynamically at the onset or late stage of senescence. Likewise, we focused our TF/target identification to seven modules of SCENT in order to find the TF affecting network in Arabidopsis leaf senescence. It has been observed that 21 TFs are present in red module followed by 12 TFs each in blue and brown modules. While, 11 TFs are present in greenyellow followed by 4 TFs each in darkolivegreen and grey60. Violet modules represent 3 TFs only. Though all TFs are scattered throughout SCENT, interestingly, four TFs are part of hub 500 genes pertaining high connections (ANAC046, ANAC019, ANAC032 and ANAC055) regulating most of the genes (Supplementary Table S4).

B. miRNA associated regulatory network
MicroRNAs (miRNAs;miRs) play significant regulatory role in multiple cellular activities throughout plant developmental progression 18 . It has been reported that miRNAs regulate the genic regulation by the RNA interference process 19,20 . To identify the role of miRNA regulation in Arabidopsis leaf senescence, we retrieved all 325 high confidence Arabidopsis mature miRs sequences from miRBase miRNA repository 21 and 2,226 antisense transcripts from TAIR online repository 22 in fasta formats. The miR-target identification was performed by psRNATarget online tool 23 with default parameters for small RNA sequences. It is observed that out of 325 miRs, only 60 miRs interact with 115 target genes generating 164 miR-target interactions to establish the functional dynamic miR-based regulatory networks in ASCENT (Supplementary  Table S4). Finally, 281 miR-target interactions were predicted on Arabidopsis leaf senescence (SCENT) by aforementioned method (Supplementary Table S4).

Reconstructing senescence responsive dynamic regulatory events
Scalable Models for the Analysis of Regulation from Time Series (SMARTS), a method which incorporates static and time series expression data to reconstruct condition-specific reaction network in an unsupervised manner 24 . First, Time Series Expression Experiments (TSEEs) is associated for datasets representing a collective biological time scale. Next, clustering is utilized to acquire a preliminary assembly for the TSEEs using associated time series. Finally, several iterations are executed to create regulatory models for sets of TSEEs and iteratively conveying TSEEs to models until convergence is achieved. Additionally, the regulatory model identifies specific stimulated pathways and genes, which uses statistical analysis to recognize TFs that vary in activity among models. Also, SMARTS can integrate miRs Dynamic Regulatory Events Miner (mirDREM) 25 , a unique miR time series expression measurement and miR-target genes for dynamic miR-regulated interaction networks. We implemented SMARTS 24 on 11,453 nuclear genes expression pattern with log2 normalization for dynamic regulatory event mining with all possible 458,302 TF/target and 281 miRtarget interactions (SCENT, Fig. 1d, Supplementary Table S5).The miR expression data was retrieved from 1 . The dynamic activated pathways regulated by TFs and miR was generated by TAIR GO ontology function. Also, we executed SMARTS 24 on 204,761 possible TF/target combinations and 164 miR-target interactions for 2,226 AGI antisense transcripts (ASCENT) and complete leaf senescence expression profile with log2 normalization (Supplementary Fig. S10a, Supplementary Table S5). In ASCENT, this analysis modeled 22 different paths that are governed by 11 different groups of TFs and several miRs ( Supplementary Fig. S10A, Supplementary Table S5, significance threshold = 0.040 for TFs and 0.251 for miRs).

Gene regulatory network simulation
Standardized Qualitative Dynamical systems (SQUAD) 26 executes steady state studies and dynamic simulations of genes and their interacting TFs by highlighting system equilibrium and their fluctuations upon regulatory stimuli (TF). The activity significance of each node was fixed to zero (initial steady state) representing the preliminary point for continual dynamic simulations to achieve the equilibrium state for network. The reciprocated impression of a stimuli (TF) is studied by modulating the initial activation to 1. Additionally, the partial initiation for input stimuli was adjusted between 0 and 1. SQUAD alter the outcome of these regulation into a time-fixed sigmoid activity curve for all connecting nodes of the network which can be visualized by trajectory plots. We implemented SQUAD simulation on integrated TF-target directional regulatory network ASCENT consisting of 237 nodes with 282 connections and SCENT containing 2,637 nodes and 5,020 connections; a total of 5,295 TF/target co-regulation (Fig. 2a). Here we display the activation states of two representative TFs i.e. ANAC046 and AT1G71520. The activation state of TF (AT1G71520) was changed from 0 to 1 to simulate the resulting trajectory plots of 22 putative target genes according to logical connectivity including partial activation of nodes ( Supplementary Fig. S10b, Supplementary Table S6). ANAC046 co-regulates with 137 genes in SCENT among which, 49 nodes exhibit hub 800 property in SCENT, whereas AT1G71520 co-regulates with 22 genes in ASCENT of which 18 nodes exhibit hub 75 property. 18 of these 49 targets of ANAC046 ( Supplementary Fig. S10c, Supplementary Table S6), and 10 of 18 targets of AT1G71520 have been reported in autophagy or senescence-associated processes. The activation of ANAC046 and AT1G71520 over time was used as an index to measure the behavior of its putative downstream target genes ( Supplementary Fig. S10b and c; Supplementary Table  S6). This analysis displayed varied levels of activation for these 137 and 22 genes for ANAC046 and AT1G71520 respectively. Schematic representation of target genes that are regulated by only ANAC046 and AT1G71520 or in combination with ANAC046, AT1G71520 and other TFs estimated through Regulatory Network Simulation by SQUAD (Fig. 2B). Similar to these two TFs, SQUAD can be applied to all the 210 TFs and 163 TFs for SCENT and ASCENT, respectively.

Plant material and growth conditions
Arabidopsis thaliana seeds (Col-0) were sown on soil and stratified for 3 days at 4°C. The seeds were then transferred to growth conditions (12h light / 12h dark at 22°C). Dark-induced senescence was conducted as previously described in 27 . Briefly, the fifth or sixth leaves of four-week old plants were detached in placed on two filter papers soaked in 10mL of distilled water in a petri dish. The petri dishes were incubated in the dark at 22°C.

Protein extraction and Western blotting
For each reported day, leaves were frozen in liquid nitrogen and ground using a mortar and pestle. Protein was extracted from leaf tissue in a buffer containing 50mM HEPES, 10mM EDTA, 50mM NaCl, 0.5% (v/v) Triton X-100, 5mM DTT, 5% glycerol (v/v), protease inhibitor (1x; P8340; Sigma), 1mM PMSF. The subsequent slurry was centrifuged at 13,000 rpm for 20 minutes at 4°C. 75µL of the supernatant was mixed with 25µL of 4x sample buffer (NP0008; Thermo Fisher Scientific) and boiled at 100°C for 5 minutes. After 2 minutes of ice the samples were centrifuged at 13,000 rpm for 2 minutes. The resulting supernatant was loaded onto a 10% SDS-PAGE gel. SAG12 was visualized after transfer using a SAG12 specific antibody (AS14 2771; Agrisera) ( Fig.  2c).

Fv/Fm quantification
Fv/Fm was used as an indicator of senescence by accessing the efficiency of photosystem II (PSII). Using genetics approach, it was shown that ANAC046 participates in leaf senescence 28 . Transgenic plants overexpressing ANAC046 (ANAC046-OX) and ANAC046 containing a SRDX domain (ANAC046-SRDX) were previously described 28 . Eight detached leaves from Col-0, ANAC046-OX, and ANAC046-SRDX were dark adapted for 15 minutes and measurements were acquired using a chlorophyll fluorometer (OS-30p+ Opti-Sciences). Measurements were taken immediately after detachment of eight leaves from four different plants of above described genotypes on days 2, 3, and 6. A two-tailed Students t-test was applied to determine statistical significance. The ANAC046-SRDX line had a p-value of 0.03 on day 6 when compared to day 6 Col-0. We determine the efficiency of the PSII at each time point using a chlorophyll fluorometer to measure Fv/Fm (Fig. 2c). At day 6, we show the ANAC046-SRDX plants is significantly higher than Col-0 (p-value 0.03), while ANAC046-OX showed reduced efficiency of PSII, although not statistically significant. These data together indicate that ANAC046 is a positive regulator of leaf senescence.

RNA extraction and quantitative real time PCR (qRT-PCR)
Total RNA from leave samples was extracted with RiboZol TM RNA Extraction Reagent (AMRESCO), followed with DNase I treatment to remove any Genomic DNA contamination. Complementary DNA (cDNA) was synthesized with the SuperScript III first-strand RT-PCR kit (Invitrogen). qRT-PCR ( Supplementary Fig. S11) was performed with gene specific primers (Supplementary Table S7) and GoTaq qPCR Master Mix (Promega) in a RealPlex S MasterCycler (Eppendorf).

Mathematical model of qRT-PCR data by C0 method
C0 is the first method in Real-time PCR analysis that does not require the assumption of equal efficiency between unknowns and the standard curve. The point of intersection between the x-axis and tangent through the inflection point of the PCR data plot is C0. Where, x is the PCR cycle number, Fx is fluorescence at cycle x, Fmax is the maximum fluorescence, d is the Richards coefficient, and Fbrf is the background reaction fluorescence.
The inflection point (Flex) was calculated as follows: And the Cγ0 was calculated as follows: C0 value is inversely proportional to amount of cDNA in the reaction and reflects the reaction kinetics since it is derived through slope of the inflection point of the raw fluorescence data. For example, C0 for SAG12 in ANAC046-OX (overexpressor) is lower than both Col-0 (wild type) and ANAC046-SRDX (knockout) at day 6, which illustrates that ANAC046-OX is kinetically more active than the Col-0 or ANAC046-SRDX. Similar C0 kinetics was observed for AT5G40690 transcript accumulation in the abovementioned genotypes ( Fig. 2e and f).

Raw datasets was fitted to the mathematical model (Equation 1), and C0 (Equation 3)
was calculated if R 2 was greater than 0.997. All calculations were done using in house python scripts.

Supplementary Discussion
Senescence is a multilayered and extremely organized aging process in plants' life cycle 1 . Plants employ a wide range of regulators that reprogram regulatory components in a highly synchronized manner to properly execute this complex developmental stage of plant life. On the contrary, the static regulatory network models demonstrate the direct transcriptional regulation without any information of time/stage specific regulation or molecular events/rewiring followed through regulators 24 . Therefore, revealing the dynamic transcriptional regulation during the onset of senescence is essential for a compressive understanding of temporal relationships among TFs and their targets. Towards this, we developed a computationally intensive pipeline to predict the regulatory networks that integrates publically available tools to study the regulation dynamics as well as the behavior of master regulators and their downstream target genes by computational simulations. With limited availability of TFs binding site information in the public domain, our sequence and co-expression based TF binding site prediction method can enrich the TF-target associations. By employing this method, we enriched the TFs-target prediction by over two folds in this study. Furthermore, the implementation of predicted TF-target association, miRs-terget interactions and transcriptome expression can reveal the stage specific dynamics of TFs and miRs regulation to predict the changes in functional association of genes. Finally, the simulated response of downstream targets and master regulators predictions based on network topology measures and ODE models can be useful for laboratory experimental validations, which will eventually reduce the time and cost of experimentation. Our computational pipeline can be an ideal model to the scientific community for predicting the TFs binding motifs and application of transcriptional regulatory networks to determine the dynamics of regulation to explore diverse plant biology challenges like abiotic and biotic stresses in crop plants.

Supplementary Fig. S1
Gene Expression analyses for nuclear encoded genes in senescence. (a) Average Information Centrality (IC) measure of darkolivegreen, greenyellow, violet, red, grey60, brown, blue and turquoise modules genes and SCENT. Darkolivegreen, greenyellow, violet, red, and grey60 modules have a higher average IC than SCENT. (b) Shortest path distribution reveals that genes in the darkolivegreen, violet, greenyellow and red modules are more close to each other than grey60, brown, blue, turquoise and SCENT (F test; p<9.8e -10 between darkolivegreen and SCENT, p= 2.18e -08 between violet and SCENT, p= 5.136e -08 between greenyellow and SCENT, p= 1.735e -06 between grey60 and SCENT, p= 1.025e -05 between red and SCENT, p= 0.0001345 between brown and SCENT, p= 0.0001265 between blue and SCENT).   To determine the TF/target pairs of TFs whose promoter binding sites are unknown, we performed a massive TF binding regulatory sequence scanning on senescence nuclear genes. First, we generated all possible combinations (forward, reverse, complement, reverse complement along with degeneracy) of nucleotides for binding of respective TF family based on text mining and family associated TF motif search. Second, we retrieved the 1000bp upstream and 5' UTR and AGI sequences of senescence nuclear genes. Third, TFs regulatory sequence were subjected to standalone blastn against 1,000 upstream 5' UTR and AGI sequences of senescence nuclear genes. 458,302 possible TF/target combinations were identified, finally we integrated coexpression networks to get the high confidence TF/targets for 220 TFs senescence.  ATG18a  AT4G25690  AT4G27720  AT4G30630  AT4G31410  AT5G05110  AT5G09630  NIK1  EMB1241  HUA2  AT5G54860  AT5G61530  TGA1  CYP86A2  AT4G39270  FTSHI1  AT4G16155 CAX7 AT5G40690

Time (arbitrary units)
Activation b c TF TF gene clusters. The green nodes denote a set of genes, mutually regulated till that interval split. The responsive group of TFs and miRNAs are listed in Supplementary  Table S5. The dynamic activated pathways regulated by TFs and miRNA was generated by TAIR GO ontology function, represented by a particular color with path numbers. This analysis modeled 22 different paths that are governed by 11 different groups of TFs and several miRNAs ( Supplementary Fig. S10A, Supplementary Table  S5, significance threshold = 0.040 for TFs and 0.251 for miRNAs). (b) Senescence Regulatory Network Simulation by SQUAD (Standardized Qualitative Dynamical systems). The simulated steady state regulatory network trajectory plots were modeled according to logical connectivity including partial activation of nodes. The change in activation state of TF (AT1G71520) from 0 to 1 stimulated the resulting trajectory plots of 22 putative target genes according to the connectivity of partial activated nodes. (c) Senescence Regulatory Network Simulation by SQUAD. The simulated steady state regulatory network trajectory plots were modeled according to logical connectivity including partial activation of nodes. The change in activation state of TF (ANAC046) from 0 to 1 stimulated the resulting trajectory plots of 27 putative target genes out of 127 (for clear visualization) according to the connectivity of partial activated nodes.