Exchange of molecular and cellular information: a hybrid model that integrates stem cell divisions and key regulatory interactions

Stem cells give rise to the entirety of cells within an organ. Maintaining stem cell identity and coordinately regulating stem cell divisions is crucial for proper development. In plants, mobile proteins, such as WOX5 and SHR, regulate divisions in the root stem cell niche (SCN). However, how these proteins coordinately function to establish systemic behavior is not well understood. We propose a non-cell autonomous role for WOX5 in the CEI and identify a regulator, AN3/GIF1, that coordinates CEI divisions. Here we show with a multiscale hybrid model integrating ODEs and agent-based modeling that QC and CEI divisions have different dynamics. Specifically, by combining continuous models to describe regulatory networks and agent-based rules, we model systemic behavior, which led us to predict cell-type-specific expression dynamics of SHR, SCR, WOX5, AN3, and CYCD6;1, and experimentally validate CEI cell divisions. Conclusively, our results show an interdependency between CEI and QC divisions. Thumbnail image


Introduction 1
Stem cells divide to regenerate themselves and to generate all of the cell-and tissue-types in a 2 multicellular organism, such as plants. The Sarkar et al., 2007). Specifically, non-cell-autonomous WOX5 maintenance of CSCs takes place 13 through the repression of the differentiation factor CYCLING DOF FACTOR 4 (CDF4) (Pi et al., 2015). 14 wox5-1 mutants have increased QC divisions in roots and a decreased number of columella cell layers 15 (Forzani et al., 2014). In the QC cells, WOX5 controls divisions by restricting CYCD3;3 expression 16 (Forzani et al., 2014). Although the regulatory modules within the CSCs and QC are well characterized 17 (Forzani et al., 2014;Stahl et al., 2013), the molecular mechanisms by which WOX5 promotes stem cell 18 fate of CEIs remains unknown. 19 Several proteins have been shown to positively regulate WOX5, such as ANGUSTIFOLIA (AN3) / GRF-20 INTERACTING FACTOR 1 (GIF1). AN3 is expressed in the root meristem with a high peak in expression 21 in the SCN and QC and plays a role in maintaining QC identity (Ercoli et al., 2018). However, whether 22 AN3 function is dependent on WOX5 and whether AN3 has a regulatory role outside the QC in the SCN 23 is not understood. Additionally, AN3 was shown to regulate the expression of SCARECROW (SCR) ( The regulatory interactions between the different cell types of the root SCN are complex and non-30 intuitive and computational tools are essential to understanding systemic behavior. Developmental 31 processes such as auxin flow within the root and lateral shoot branching have been mathematically 32 modeled to better understand and predict system-level behavior (Canher et al., 2020;Prusinkiewicz et 33 al., 2009). Some models implement different scales of the system to simulate, understand, and predict system-level behavior as a whole. For example, a mathematical model that simulates and predicts the 35 induction of shoot branching during plant development included on a molecular scale auxin flux across 36 metamers (i.e. smaller segments of the stem) and on an organ scale the formation of metamers of the 37 stem and lateral branches (Prusinkiewicz et al., 2009). Modeling systems and allowing exchange of 38 information across different scales can also be achieved by combining agent-based models (ABM) with 39 continuous models, such as ordinary differential equations (ODEs) or partial differential equations 40 (Cilfone et al., 2015). ABMs consist of autonomous "agents" that dynamically interact and show 41 responsive behavior through a set of simple rules. ABMs have, for example, been used to simulate 42 plant-herbivore interactions (Radny & Meyer, 2018). However, within the molecular plant biology field, 43 these models are not widely used, despite their capacity to capture system-level behavior. On the 44 other hand, continuous models such as ODEs have been applied to infer gene regulatory networks 45 (Krouk et al., 2010;Yao et al., 2011) and predict dynamic gene expression patterns (Clark, Fisher, et al., 46 2020). These models are computationally intensive and lack the capability to capture system-level 47 behavior but can model complex dynamic responses over time. Hybrid models are created when, for 48 example, continuous models are used within a discrete ABM to describe a part of the system. These 49 hybrid models are usually multiscale models, given that the continuous models often describe a 50 dynamical response on a different spatiotemporal scale than the ABM (Cilfone et al., 2015). 51 In this study, we combine cell-type-specific gene expression data and experimental data with network 52 inference and parametric models to better understand how WOX5, AN3, SCR, and SHR coordinately 53 regulate CEI stem cell divisions. We transcriptionally profiled CEI cells in wild-type and wox5-1 roots, 54 as well as QC cells and non-stem cells. We found that AN3 was among the most CEI-enriched genes. 55

WOX5 regulates CEI-specific genes 66
The functional role of WOX5 in the QC and CSC has been extensively reported while its role in stem 67 cell populations remains largely unknown. WOX5 is specifically expressed in the QC cells, however, the 68 protein moves to the CSCs and the vasculature initials and has been shown to have a non-cell 69 autonomous role in these cells (Clark et al., 2019;Pi et al., 2015). To determine whether WOX5 is also 70 able to move from the QC cells to the QC-neighboring CEI cells and regulate downstream targets, we 71 used scanning fluorescence correlation spectroscopy (scanning FCS). Five-day-old 72 wox5xpWOX5:WOX5-GFP plants were analyzed with scanning FCS to evaluate the directional 73 movement of WOX5 protein between these two cell-types. Line scans were taken over time from a 74 region spanning the CEI and adjacent QC ( Fig 1A). This analysis resulted in a quantitative assessment 75 of movement and allowed us to calculate the movement index (MI). We found that WOX5 moved 76 bidirectionally between the QC and the CEI (MI = 0.90 ± 0.04 from QC to CEI, MI = 0.83 ± 0.05 from CEI 77 to QC, n = 20) (Supplemental Table 1). As a comparison, within the SCN, free GFP and immobile 3xGFP 78 have a moving index of ~0.7 and ~0.25, respectively (Clark et al., 2016). 79 To explore the potential functional role of WOX5 in CEI, we examined the expression pattern of the 80 CEI-marker pCYCD6;1:GFP in wox5. The marker showed an expression pattern that extended into the 81 cortex and endodermal cells (Fig 1B,C). This expanded expression of CYCD6;1 suggests that the 4 to 5 82 cells proximal of the CEI, referred hereafter as CEI-like cells, have gained stem cell-like characteristics 83 and also indicates that WOX5 controls CYCD6;1 expression to the CEI (Fig 1B). We then explored the 84 role of WOX5 in limiting CYCD6;1 expression and, thus, controlling CEI divisions. To this end, we 85 quantified the number of undivided and divided CEI cells in 4-, 5-, and 6-day-old wox5 and wild-type 86 roots. This quantification showed that wox5xpCYCD6;1:GUS-GFP roots had an increase of 23.43% and 87 25.33% (p = 0.0495, Wilcoxon test) divided CEI cells compared to the wild type (WT) at 4 and 6 days, 88 respectively ( Fig 1D). Taken together, these results support a functional non-cell autonomous role for 89 WOX5 in the CEI.

Network inference and node importance analysis to identify functional candidates 100
To unravel the transcriptional events regulating the extended expression pattern of CYCD6;1 in the 101 wox5 mutant background, a transcriptome analysis was performed on FACS-sorted GFP positive cells 102 from pCYCD6;1:GUS-GFP, wox5xpCYCD6;1:GUS-GFP, and pWOX5:GFP, and the meristematic cells from 103 pWOX5-GFP that do not express the marker (referred to as non-stem cells) (Supplemental Table 2). 104 Compared to the cells not expressing the pWOX5-GFP marker, 163 genes were differentially expressed 105 (FDR < 0.05) in wild-type CEI cells and 213 genes in the CEI and CEI-like cells from the wox5 mutant. In 106 total, the union of these two analyses identified 330 DEGs in CEI and CEI-like cells, of which 159 DEGs 107 (48.18%) have previously been shown to be expressed in the SCN and 53 genes were enriched in the 108 CEI (Clark et al., 2019). We hypothesized that the regulatory genes underlying CYCD6;1 expression 109 should be differentially expressed in the CEI cells (CYCD6;1 expressing cells) of the wild-type and wox5 110 roots and thus focused on the genes overlapping between these two sets of DEGs (Fig 2A). In total, 46 111 genes overlapped between the CEI and CEI-like cells, which equals an enrichment of 35.8 (p < 4.431e-112 59, Exact hypergeometric probability). To identify key regulatory proteins among these 46 genes, we 113 predicted causal relations between the TFs and downstream genes with high accuracy and constructed 114 a gene regulatory network. We inferred the causal relations by leveraging our transcriptome data with 115 a regression tree algorithm RTP-STAR ( Fig 2B)  (GIF1), which is a known regulator of cell proliferation, and an unknown TF (AT1G75710). Among the 120 inferred AN3 targets, we confirmed with TChAP data that three targets (AT1G75710, FLA10, and 121 GBSS1) were directly bound by AN3 (Vercruyssen et al., 2014). Network inference allowed us to identify 122 potential functionally important genes, however, we still needed to pinpoint the biological important 123 genes within the network. 124 To identify which genes could cause the largest impact on network stability when perturbed, we 125 performed a node importance analysis. To calculate the impact of each gene, each node received a 126 weight depending on its outdegree (i.e. number of outgoing edges), then for each node, the sum of 127 the weighted outgoing first neighbors and the sum of the weighted incoming first neighbors was taken. 128 Both sums were in turn weighted, specifically, the sum of the outgoing neighbors was weighted by 129 Average Shortest Path Length (ASPL), and the sum of the incoming neighbors was weighted according 130 to the proportion of end-nodes within the network, which is in this network 20% (see Materials and 131 Methods). We next developed an R-based Shiny application (Node Analyzer) that calculates the 132 weights and impacts of each gene within a network (Shannon et al., 2003) (see Materials and Methods) 133 (Supplemental Fig 1). Node Analyzer allowed us to rank the 20 genes in the network and select key 134 genes. The most impactful gene within our network is AN3, a transcriptional coactivator that is 135 involved in cell proliferation during leaf and flower development ( Fig 2C).

AN3 contributes to the regulation of CEI divisions 145
It was previously shown that AN3/GIF1 and its closest homologs, GIF2 and GIF3, were expressed in the 146 root stem cell niche (Ercoli et al., 2018). A triple mutant (gif1/2/3) displayed a disorganized QC and 147 increased root length as a result of an increased root meristem size (Ercoli et al., 2018). We confirmed 148 the growth repressing role of AN3 in the roots, as an3 and 35S:AN3-GFP roots showed an increased 149 and reduced root length compared to the WT, respectively (Supplemental Fig 2A). We observed a 150 disorganized stem cell niche in 56% (25/45 roots) of an3 mutant roots (Supplemental Fig 2B). 151 Additionally, an3 mutants contained starch granules in the cells that are normally CSC, suggesting that 152 AN3 plays a role in CSC maintenance (Supplemental Fig 2C). To determine whether AN3 also plays a 153 role in CEI divisions, we quantified the number of undivided and divided CEI cells in 4-, 5-, and 6-day-154 old an3 and WT roots. 6-day-old an3 roots had 19.22% fewer undivided CEI cells at compared to WT 155 (p = 0.103, Wilcoxon test), suggesting that more CEI divisions occur in the an3 mutant ( Fig 3A). 156 Additionally, when an3 is crossed with the CEI-marker pCYCD6;1:GUS-GFP, an extended expression 157 pattern is observed (Fig 3B,C). Taken together, these results support a role for AN3 in the regulation 158 of CEI divisions. 159 160 Figure

A hybrid model to dynamically simulate and predict stem cell divisions 164
If AN3 and WOX5 are indeed key regulators for CEI divisions, we would expect that their temporal 165 expression influences CEI divisions in a cell-type specific manner. To gain insight into the system-level 166 regulation of CEI stem cell divisions, we modeled the expression of CYCD6;1 and its direct and indirect 167 upstream regulators: SHR, SCR, WOX5, and AN3 ( Fig 1C, Fig 3B) (Sozzani et al., 2010). For this we 168 developed a hybrid model that combines agent-based modeling aspects with ODEs. Specifically, we 169 included four different cell types or "agents" (QC, CEI, vascular initial, and endodermal cell) and 170 constructed ODEs of the genes for each cell type that are able to recapitulate the dynamics of the 171 upstream regulatory interactions at a molecular scale. The cells/agents interact through the movement 172 of SHR and WOX5 and change state (i.e. divide) upon changes in the expression of specific proteins. 173 For example, when CYCD6;1 exceeds a certain abundance, the CEI will divide. Each time a cell divides 174 (an agent changes state), corresponding protein abundances are halved. As such, we were able to 175 exchange information bidirectionally, from molecular to cellular scale and from cellular to molecular 176 scale. To implement this hybrid model we used SimBiology to model, simulate, and analyze dynamic 177 systems that allows for rapid model optimization and provides an intuitive visualization of the model 178 (The MathWorks, 2019). 179 To analyze the temporal expression dynamics of CYCD6;1 linked to CEI divisions, and to understand 180 the regulatory role of WOX5 and AN3 in controlling the CYCD6;1 dynamics, we used ODEs to generate 181 a quantitative model that describes the dynamics of four key transcriptional regulators of CYCD6;1, 182 namely WOX5, AN3, SHR, and SCR. In our ODE systems, each ODE included a degradation term and a 183 production term that depended on its upstream regulations. The included regulations are depicted in 184 Additionally, we included ODEs that model the movement of WOX5 from the QC to the vasculature 191 initials (Supplemental Table 3), different diffusion rates of SHR from the vascular initials to the 192 endodermis and QC (Clark, Fisher, et al., 2020), the SHR/SCR complex formation, and the oligomeric 193 states of WOX5 and AN3. The oligomeric states of AN3 and WOX5 were experimentally determined 194 using scanning FCS (Supplemental Fig 3). Specifically, we performed Number and Brightness (N&B) on 195 an3 or wox5 roots expressing pAN3:AN3-GFP or pWOX5:WOX5-GFP translational fusion, respectively. 196 We found that both AN3 and WOX5 primarily exist as a monomer (98.67% and 96.01%, respectively) 197 with a very small amount of dimerization (1.33% and 3.99%, respectively) (Supplemental Fig 3). Thus,198 we fixed the oligomeric state of AN3 and WOX5 as monomers in our ODE model. As  We estimated the values for the sensitive parameters by fitting our model to computed cell-type 215 specific time course data (Supplemental Table 5 Table 5). After estimating the 219 sensitive parameters, we simulated the hybrid model to evaluate the expression dynamics within each 220 cell. For example, the hybrid model predicts high expression of SCR in the endodermal cells and a lower 221 expression in the CEI and QC. We confirmed the increased SCR expression in the endodermal cells by analyzing confocal images of the QC, CEI, and endodermal cells of pSCR:SCR-GFP for corrected total 223 cell fluorescence (CTCF) at 5 days 16 hours (Supplemental Fig 5A,B). Model simulations showed that 224 the cell-specific networks ensured robust stability of cellular behavior, such as cell division regulation 225 ( Fig 4B). The agent-based rules for cell division were set based on SHR/SCR complex and WOX5 226 expression for the QC and CYCD6;1 expression for the CEI (Supplemental Fig 6). Our hybrid model was 227 able to capture a dynamic expression pattern for the SHR/SCR complex, with high expression at 4 days 228 8 hours and 5 days 16 hours. In contrast, WOX5 shows a low expression at these time points 229 (Supplemental Fig 5C). The first peak of SHR/SCR expression at 4 days 8 hours was previously shown in 230 an ODE model, while the second peak occurred, compared to our model, earlier at 5 days 8 hours 231 (Clark, Fisher, et al., 2020). Model predictions show that the fine balance between low expression of 232 the SHR/SCR complex and WOX5 simulates a QC cell division at 5 days 5 hours. Indeed, 5-and 6-day-233 old plants show an increase in QC divisions compared to 4-day-old plants (Supplemental Fig 5D). 234 Additionally, CEI divisions are predicted to occur at 4 days 8 hours and 5 days 16 hours ( Fig 4B). We 235 observed an increased percentage of divided CEIs in 5-day-old roots compared to 4-day-old roots, 236 however, an increase was not visible in 6-day-old roots compared to 5-day-old roots ( Fig 1D, Fig 3A). 237 We found that the rate of CEI divisions within our model was influenced by the QC division. conditions. However, to further validate the model, we simulated the loss-of-function of wox5 and an3 255 and evaluated the expression patterns as well as CEI division dynamics. Based on transcriptome data of wox5 and an3, we calculated an 99.53% and 88.12% reduction of WOX5 and AN3 expression in their 257 respective loss-of-function lines (Supplemental Fig 7). As such, the initial expression levels of WOX5 258 and AN3 were set to 0.47% and 11.88% in the mutant simulation as compared to the values in a WT 259 situation, respectively. 260 Model simulations of wox5 loss-of-function predicts an additional CEI division between 4 and 5 days 261 compared to WT, which coincides with an increase in divided CEI cells at 4 days in wox5 (Supplemental 262 Fig 8A, Fig 1D). The additional division is most likely the result of the removal of WOX5 repression on 263 SHR in the vascular initials leading to an accelerated accumulation of SHR/SCR complex in the CEI. An 264 overall increase in SHR/SCR in the CEI was not predicted by the model (Supplemental Fig 9B), and 265 accordingly, CEI-specific transcriptomics and protein quantifications in the CEI of the wox5 mutant did 266 not show an increased SHR expression (Supplemental Table 2, Supplemental Fig 9A). The simulations 267 of the an3 loss-of-function predict the depletion of SCR in the QC, CEI, and endodermal cell compared 268 to WT (Supplemental Fig 8B). This decrease in SCR expression has been shown within the QC (Ercoli et 269 al., 2018). However, the CEI and endodermis still showed high levels of SCR when a repressor version 270 of AN3 is expressed in the SCR reporter line (Ercoli et al., 2018), which is in contrast to the model 271 predictions. As such, the regulation of CYCD6;1 by AN3 in the CEI may not be established via SCR but 272 another unknown mechanism. We hypothesized that AN3 is regulating an additional factor that 273 represses CYCD6;1. For this, we added an unknown factor X that is activated by AN3 and represses 274 CYCD6;1, removed the AN3 activation of SCR, updated the ODEs within the CEI agent accordingly, and 275 re-estimated 4 former and 2 new parameters (see Materials and Methods) (Supplemental Table 7 By adding competition between a repressor, transcriptionally activated by AN3, and the SHR/SCR 279 direct regulation of CYCD6;1, the model was able to accurately capture the CEI divisions in a wild-type 280 situation as well as in an an3 mutant background (Fig 5A). Notably, by adding the repressor to the 281 model, the CEI division time interval shortened to 23.3 hours (Supplemental Fig 10). To identify 282 potential candidates as a repressor downstream of AN3, we performed genome-wide expression 283 analysis on an3 meristematic root tissue (Supplemental Table 9). In total 1013 genes were differentially 284 expressed (q < 0.05) including 67 TFs of which 4 TFs were shown to interact with TOPLESS (TPL), a 285 known transcriptional corepressor (Causier et al., 2012) (Fig 5B). Of these 4 transcriptional repressors, 286 WRKY30 and MYB7 showed the highest expression correlation with the model prediction (Fig 5C).  between these cell types. WOX5 proteins can move to the neighboring vascular initials and CEI cells and SHR proteins move to the QC, CEI, and endodermal cells. Scanning FCS was used to quantify the 342 diffusion coefficient of WOX5 and SHR to include into the model (Supplemental Table 7 For imaging and root growth assays, seeds were dry sterilized using fumes produced by a solution of 377 100% bleach and 1M hydrochloric acid. The seeds were plated on square Petri dishes with solid (10 g/L 378 agar, Difco TM ) 1X MS (Murashige and Skoog) medium supplemented with 1% sucrose and stratified for 379 2 days at 4°C. The plates were grown vertically at 22°C in long-day conditions (16-hrs light/ 8-hrs dark) 380 for 4, 5, 6, or 7 days as indicated in the figures. At least three biological replicates of 10 to 20 plants 381 were performed for the root growth assays and confocal images. The different lines were always grown 382 together on one plate with the appropriate control line. For RNAseq experiments, seeds were wet 383 sterilized using 50% bleach, 100% ethanol, and water. Seeds were imbibed and stratified for 2 days at 384 4°C. Next, the seeds were plated with high density on Nitex mesh squares on top of solid 1X MS 385 medium with 1% sucrose. Seeds were plated and grown vertically at 22°C in long-day conditions. 386

Root growth assays 387
At 3, 4, 5, 6, and 7 days, the primary root length was marked. At 7 days, a picture of the marked square 388 plates was taken and the root length was measured using the software program ImageJ version 1.45 389 (National Institutes of Health; http://rsb.info.nih.gov/ij/). For the statistical analysis of the root growth 390 assays, Student's t-tests were performed on the average of each biological replicate. 488nm and 570nm lasers were used for green and red channel acquisition, respectively. Propidium 394 iodide (10μM, Calbiochem) was used to stain cell walls and mPS-PI staining was used to visualize starch 395 granules. For the N&B acquisition, 12-bit raster scans of a 256x256 pixel region of interest were 396 acquired with a pixel size of 100nm and a pixel dwell time of 12.61μs as described in ( to contribute more to the impact of a node. Because the first neighbors are weighted in regards to 451 their outdegree, genes with a lower outdegree can still have a large impact if its neighbors have a high 452 outdegree and the gene is thus centrally located. Genes with a large number of cascading targets that 453 are 2 or more nodes away will have a higher ASPL and thus a higher scaled outdegree weight, 454 accurately reflecting the hierarchical importance of the source gene itself and its first neighbors 455 targets.

Shiny app: Node Analyzer 457
To calculate necessary network statistics such as outdegree and indegree in Cytoscape® 3.8.0 (Shannon 458 et al., 2003), select Tools -> Analyze Network, check the Analyze as Directed Graph if applicable, and 459 then press OK to perform the analysis. To export node and edge files from Cytoscape, select File -> 460 Export -> Table to File, and then choose default edge or default node in the 'Select a table to export' 461 dropdown. Press OK to export each file. Import the node and edge table files into the corresponding 462 prompts (Fig 2C)

495
It was shown that the different oligomeric forms and stoichiometries of SHR, SCR, and the SCR-SHR 496 complex show a similar expression pattern (Clark, Fisher, et al., 2020). As such, the SHR and SCR 497 oligomeric forms were modeled as one variable. 498 The interaction between the different agents/cell types is modeled using mass-action kinetics. The 499 state change following division is modelled using simple agent-based rules. To simulate division of an 500 agent, the capacity of the cell doubles, subsequently halving all proteins present. 501  (Fig 1)

513
For the sensitivity analysis, the total Sobol effect index was calculated for each parameter value (Saltelli 514 et al., 2010;Sobol′, 2001). Parameter values were randomly sampled using Monte Carlo sampling to 515 obtain 150 different values for each parameter. This analysis was repeated for 10 technical replicates. 516 As such, for each parameter 170 (10 replicates x 17 ODEs) total Sobol effect indices were obtained. For 517 each ODE and replicate the sensitivities were rescaled between 0 and 1 and then averaged across the 518 17 ODEs. The obtained averaged sensitivities for each replicate were again averaged to retrieve the 519 total Sobol effect index per parameter (Supplemental Table 4). The sensitive parameters were chosen 520 as the parameters that had significantly higher Sobol indices than the lowest scoring parameter 521 (K_D2_qc) using a student's t-test (p<0.01). 522 To estimate the sensitive parameters, the model was fitted onto extrapolated cell-type specific time 523 course expression data (Supplemental Table 5 were able to extrapolate cell-type specific time course expression values (Supplemental Table 5). 529 Simulated annealing and Latin hypercube sampling as described in (Clark, Fisher, et al., 2020) produced 530 40 sets parameter estimates (Supplemental Table 6). The average of these parameter estimates was 531 used for the model simulations. The remaining sensitive parameters were set to a constant value from 532 the corresponding estimated parameter in (Clark, Fisher, et al., 2020). The value of non-sensitive 533 parameters was selected based on similar values of the model described in (Clark, Fisher, et al., 2020). 534 The production terms for WOX5 (k1_qc) and AN3 (k2_qc, k2_cei, k2_endo) were set to a constant value 535 at each time point to minimize the error between the model and the time course expression data. The 536 diffusion coefficients of SHR (a_qc, a_cei) and WOX5 (b_qc) were experimentally determined from RICS 537 experiments (Supplemental Table 3) (Clark, Fisher, et al., 2020). 538 The following changes were made in the regulatory network underlying the CEI divisions to reflect the 539 an3 loss-of-function in the hybrid model: 540 (1) Factor X; for the upstream regulation of the unknown repressor X in the CEI agent, the 541 activation by AN3 was included. 542

543
(2) CYCD6;1; for the upstream regulation of CYCD6;1 expression, we added the repression of 544 factor X in addition to the activation by the SCR-SHR complex (Sozzani et al., 2010

555
Four existing parameters (k3_endo, d3_endo, k3_cei and k5_cei) and two new parameters (k6_cei, 556 d6_cei) were re-estimated in the same manner as described above and produced 20 sets parameter 557 estimates (Supplemental Table 8). For the remaining parameters the same value as the initial hybrid 558 model was used. 559 All parameters for the initial and adjusted model are listed in supplemental table 7. To simulate the 560 hybrid models, the initial values were set as the 4D FPKM values from the extrapolated time course 561 data. For factor X, the SHR/SCR complex, and very lowly expressed genes (e.g. WOX5 in the vascular 562 initials) the initial value was zero. To simulate wox5 loss-of-function the initial value of WOX5 was set 563 to 0.47% (Supplemental Fig 7). To simulate an3 loss-of-function the initial value of AN3 in all three 564 agents, was set to 11.88% (Supplemental Fig 7). ODE45 was used as the ODE solver within SimBiology. 565

Data and Coding Availability 566
All sequencing data are available on GEO at: 567

Conflicts of Interest declarations in manuscripts 586
The authors declare no conflict of interest.