A mathematical model for understanding synergistic regulations and paradoxical feedbacks in the shoot apical meristem.

The shoot apical meristem (SAM) is the primary stem cell niche in plant shoots. Stem cells in the SAM are controlled by an intricate regulatory network, including negative feedback between WUSCHEL (WUS) and CLAVATA3 (CLV3). Recently, we identified a group of signals, Epidermal Patterning Factor-Like (EPFL) proteins, that are produced at the peripheral region and are important for SAM homeostasis. Here, we present a mathematical model for the SAM regulatory network. The model revealed that the SAM uses EPFL and signals such as HAIRY MERISTEM from the middle in a synergistic manner to constrain both WUS and CLV3. We found that interconnected negative and positive feedbacks between WUS and CLV3 ensure stable WUS expression in the SAM when facing perturbations, and the positive feedback loop also maintains distinct cell populations containing WUS on and CLV3 on cells in the apical-basal direction. Furthermore, systematic perturbations of the parameters revealed a tradeoff between optimizations of multiple patterning features. Our results provide a holistic view of the regulation of SAM patterning in multiple dimensions. They give insights into how Arabidopsis integrates signals from lateral and apical-basal axes to control the SAM patterning, and they shed light into design principles that may be widely useful for understanding regulatory networks of stem cell niche.

Introduction

Spatiotemporal regulation of cell populations in stem cell niches is critical for plants and animals during developmental and homeostatic conditions [1], [2], [3], [4], [5]. The shoot apical meristem (SAM) of Arabidopsis thaliana is one of the most extensively studied model systems for understanding stem cell niches in plants. The dome-shaped SAM harbors stem cells that are located in the epidermal and subepidermal layers. Under normal conditions these cells are characterized by the high expression of CLAVATA3 (CLV3), which encodes a secreted peptide that can diffuse to neighboring cells and regulate their gene expression [6]. The patterning and size of the SAM remain stable in mature plants, suggesting tight regulation of stem cell populations. Previous studies have revealed crucial mechanisms for regulating these cells in the apical-basal axis of the SAM. Particularly, the number of stem cells is maintained by WUSCHEL (WUS), which is expressed in cells of the organizing center below the subepidermal layers of the SAM [7], [8]. WUS protein moves to the upper layers and activates CLV3 expression, whereas the diffusible CLV3 peptide suppresses the expression of WUS [9], [10], [11]. This negative feedback loop is considered the core circuit controlling the size of the SAM and the number of stem cells in this region [9], [10], [12]. More recent discoveries show that the separation of WUS and CLV3 expressing regions is driven by additional signals originating from the basal zones of the SAM, such as high concentrations of WUS and/or HAIRY MERISTEM (HAM) [13], [14], [15], and these signals prevent CLV3 expression below the subepidermal layer. Mathematical models have been used in numerous studies to describe the dynamics of gene activities in the SAM quantitatively, and have improved our understanding of the complex molecular and cellular networks controlling SAM patterning [11], [13], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25]. These models have revealed the key role of negative feedback in maintaining homeostasis, as well as the mechanisms underlying the separation between CLV3 and WUS expressing domains. While these experimental and computational studies provided ample information about regulatory networks along the apical-basal axis, regulatory mechanisms in the lateral direction of the SAM remain unclear. For example, it was proposed that some ‘pocket-like’ signals restrict CLV3 expression in both the apical-basal and lateral directions [22], but currently there is no experimental evidence showing the existence of any inhibitory factor with such an expression pattern. Although WUS has been shown to play a role in restricting the CLV3 expression laterally [13], [15], [26], [27], theoretical reconstruction of correct SAM patterning still requires some additional hypothetical ‘positional cues’ that constrain CLV3 expression [27]. Therefore, our understanding of mechanisms underlying lateral boundary of the CLV3 and/or WUS expressing domains in the SAM is incomplete.

Recently, we have identified a group of potent signals, Epidermal Patterning Factor-Like (EPFL) proteins, that originate from the peripheral zone of SAM and inhibit CLV3 and WUS expression [28], [29]. EPFL signals are communicated by four diffusive ligands: EPFL1, EPFL2, EPFL4, and EPFL6 [28]. The loss of these ligands caused the upregulation of both WUS and CLV3 and increased the size of the SAM, suggesting that these signals play an essential role in maintaining the stem cell population in the SAM, and their perturbation can override the negative feedback loop [29]. However, it is unclear how the SAM integrates information from the apical-basal direction with that from the lateral direction, and how multiple signals synergistically control SAM patterning.

To obtain insights into these questions, we built a mathematical model describing the intracellular and intercellular regulatory networks in the meristem. We perturbed the model to examine the key regulatory elements (e.g. EPFL and HAM) controlling the gene expression patterns in the middle and peripheral regions of the SAM, and to unravel the roles of interconnected negative and positive feedbacks between WUS and CLV3 in SAM patterning. We further used three new metrics of patterning to explore functional objectives that may constraint the kinetic rate constants underlying regulations of key genes: 1) downregulation of WUS expression in the lateral region of the SAM; 2) stability of total WUS concentration in the SAM in the presence of perturbations; and 3) diverse cell populations that highly express either WUS or CLV3 in the middle of the SAM. These analyses provide a holistic view of the principles governing the SAM patterning.

Results

A mathematical model recapitulates known phenotypes of SAM patterning under normal and perturbed conditions

Previous mathematical models of SAM patterning have been focused on the negative feedback between WUS and CLV3, as well as other factors controlling their expressions along the apical-basal axis [11], [13], [16], [17], [18], [19], [20], [22], [23], [24], [25]. Here, we built a two-dimensional model incorporating elements of these previous models with the potent lateral regulator EPFL that was recently characterized (Fig. 1A) [28], [29]. We considered 51 cells (a minimum number of cells that reflects multilayer patterning of the SAM in terms of key factors [13]) organized in a dome-shaped structure. This structure contains six layers of cells (L1-L6) for which gene expression patterning is critical for development [13], [14]. In each cell, concentrations of six interacting molecules were described with ordinary differential equations (ODEs) based on the gene regulatory network shown in Fig. 1A (see Methods for details). Key assumptions about the regions producing the regulatory molecules are illustrated in Fig. 1B (See Methods and Table 1 for details). Our assumptions also include diffusion/movement of EPFL, WUS and CLV3 in all directions. We sought to obtain a comprehensive view of gene regulation in SAM patterning by simulating and analyzing this model under various conditions.

Fig. 1

A mathematical model of SAM captures roles of key regulators of SAM patterning. A. Influence diagram for key molecules regulating SAM patterning. ARR and cytokinin were implicitly modeled (WUS activates itself in the model). A single variable was used to describe the combined signaling of WUS and HAM. Arrows indicate experimentally supported regulations (Table 1). B. Key assumptions about spatial distribution of SAM regulators. Black regions highlight cells that were assumed to have the capacity to produce the denoted molecules. Maximum signaling/distribution/production refers to the distributions of denoted production or signaling in the absence of inhibitors and in the presence of activators. EPFL, CLV3 and WUS proteins were allowed to move between neighboring cells. C-J. Steady state distributions of SAM regulators under various conditions. All mutants were simulated by removing the production of the denoted molecules. All simulations were performed under initial conditions with low amount (=0) of all molecules across SAM.

Table 1

Key assumptions of the model.

Assumption	Reference for experimental evidence
Ligands EPFL are produced in the lateral region of the SAM, but are diffused broadly in the SAM	[28]
EPFL inhibits CLV3 and WUS expression	[29]
CLV3 inhibits WUS expression, and WUS activates CLV3 expression	[9], [10], [11]
WUS activates its own expression (for example: an ARR-Cytokinin-WUS loop)	[23], [39], [58], [59]
HAM inhibits CLV3 expression	[13]
HAM and WUS synergize to exert transcriptional regulation	[13], [49]
HAM is expressed in the lower middle region of the SAM	[13]
WUS may be expressed broadly in the SAM, except for L1 layer	[9], [10]
CLV3 may be expressed broadly in the SAM	[13], [56], [63]
WUS and CLV3 move/diffuse to neighboring cells in the SAM	[11], [55], [56]

* The phrase ‘may be expressed’ is used when the model allows the expression, but it is inhibited by relevant inhibitors when they are present.

We first fit our model to SAM patterning under normal conditions (wild type, WT) as well as several phenotypes in terms of distributions of key molecules due to genetic perturbations of several regulators. We identified a representative parameter set that allows the model to reproduce experimentally observed phenotypes (Fig. 1C–I, Table 2).

Table 2

Known SAM gene expression patterns captured by the model.

Genetic background	Description of patterning	Reference for experimental evidence
Wild type	CLV3 is expressed in the top middle region of the SAM	[6], [10], [31], [56], [57]
Wild type	WUS is expressed in the lower middle region of the SAM	[6], [8], [9], [10], [57]
wus mutant	No CLV3 expression	[30], [31]
clv3 mutant	Increased WUS expression in the SAM as compared to wild type, including lateral region	[9], [63]
erf mutant	Increased WUS expression in the SAM as compared to wild type, including lateral region	[32]
erf mutant	Increased and broader CLV3 expression as compared to wild type	[30], [32]
clv3-erf mutant	Highly increased WUS expression as compared to wild type	[29], [30]
ham mutant	CLV3 expression is centered at in the lower middle region of the SAM	[13], [33]
Induction of CLV3 expression by 10-fold	Stable WUS expression pattern in the SAM	[34]

In the simulated wild type, WUS and CLV3 were expressed in confined regions near the central axis of the SAM, while their proteins diffused (CLV3) or moved (WUS) broadly to other regions (Fig. 1C). Along the apical-basal axis, a moderate amount of WUS moved to L1 layer, where it maintained the expression of CLV3. Similarly, a moderate amount of CLV3 diffused to the meristem rib (Fig. 1C). Notably, both WUS and CLV3 expressing regions are limited by lateral boundaries, whereas their expressions in the apical-basal direction are anticorrelated with each other, with the L1 layer expressing CLV3 but not WUS, and bottom layers only expressing WUS (Fig. 1C). EPFL was distributed differentially within and across layers, due to the assumption of its isotropic diffusion, and this established a lateral-middle spatial gradient (Fig. 1C). In the wus mutant, no CLV3 expression was observed as we expected from the gene regulatory network (Fig. 1A and D) [30], [31]. The lateral confinement of both WUS and CLV3 was lost upon the removal of signaling triggers by EPFL (Fig. 1E): WUS expression had a moderate lateral expansion, whereas CLV3 had a dramatic expansion in the epidermal (L1) layer and the L2 layer (note that the erf mutant in this study refers to the absence of multiple ERf receptors, which requires knocking out three genes experimentally). These observations are consistent with recent experimental findings [30], [32]. Interestingly, the model showed that the anticorrelation pattern of WUS and CLV3 expressions was largely maintained in this mutant despite their own expansions (Fig. 1E). The model reproduced the well-known phenotype of the CLV3 knockout (KO): the WUS expressing region expanded both apically and laterally upon the loss of CLV3 (Fig. 1F). Moreover, the expression in the meristem rib was increased from the level in the wildtype (Fig. 1C and F), indicating the inhibitory effect of CLV3 diffused to this region. The model showed that in the absence of both EPFL and CLV3 signals, WUS had a dramatic increase in expression across the SAM (Fig. 1G). The model also showed that the anticorrelated WUS-CLV3 expression pattern was lost in the absence of HAM (Fig. 1H), and this co-expression of WUS and CLV3 in the organizing center with HAM KO was observed previously [13], [33]. Interestingly, there was a moderate decrease of CLV3 expression in the central zone with the loss of HAM compared to the wild type (Fig. 1H), even though the central zone CLV3 expression was not regulated by HAM directly in our model (see Methods). The decrease of CLV3 expression in the central zone was because of the negative feedback between WUS and CLV3: the moderately decreased WUS expression upon the expansion of CLV3 expression to the meristem rib reduced the availability of WUS protein in the central zone. This reduction in turn decreased the CLV3 expression in the central zone, but not in the meristem rib due to the dominant effect of the loss of HAM. The alteration of the CLV3 expression pattern in the ham mutant is consistent with a recent report (Zhou et al., 2018). Nonetheless, a significant level of WUS expression was maintained in the meristem rib (Schulze et al., 2010), and this level of WUS in the SAM was critical for the overall gene patterns in the region, as manifested in the phenotype upon the complete loss of WUS (Fig. 1D). In addition to these genetic perturbations, our model recapitulated the observed stability of the WUS expression when the CLV3 expression was increased by 10-fold [34] (Fig. S1A). We summarized the key SAM expression patterns that are captured by the model and their supporting experimental evidence in Table 2. In conclusion, our model captures the effects of perturbations of key genes that regulate SAM patterning in both apical-basal and lateral directions.

The model makes an immediate prediction that the loss of HAM and EPFL will cause dramatic expansion of the CLV3-expressing region, as well as expansion of the WUS-expressing region in both the lateral and apical-basal axes (Fig. 1I). In the following sections, we describe specific insights and other predictions that we obtained from the model.

Differential roles of EPFLs and CLV3 in regulating SAM patterning

It has been previously shown that WUS expression in the middle SAM is essential for maintaining stem cell populations [7], [8], [26]. In contrast, WUS expression in the lateral region of the SAM is associated with abnormal development [10], [30]. These observations suggest that the two regions may require different modes of regulation. To gain deeper understanding of the roles of EPFL and CLV3 in limiting the WUS expression region, we focused on two SAM areas described in the model (Fig. 2A). Region 1 (Fig. 2A, blue) includes cells near the middle apical-basal axis of the SAM, whereas Region 2 (Fig. 2A, Indian red) includes cells in the peripheral zone. In terms of our simulation results (Fig. 1C), Region 1 was defined as cells that express either WUS or CLV3 under the wild-type condition, whereas Region 2 was defined as cells that neither express WUS nor CLV3 under the same condition.

Fig. 2

EPFL and CLV3 control patterns of WUS expression at different regions. A. An illustration of two regions that were analyzed for WUS mRNA production. Region 1 (blue) contains cells near the apical-basal axis of the SAM, whereas Region 2 (Indian red) contains cells at the lateral boundary between the middle zone and the peripheral zone. B. Total WUS mRNA in the two regions indicated in A. Simulations were performed with the model structure described in Fig. 1. Mutant genotypes were modeled by removing the production of the indicated mRNAs. The hatched bars were calculated by adding the amounts of mRNA in the erf mutant and clv3 mutant increased from wild type, and the amount of wild-type mRNA. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We quantified the amount of WUS mRNA (in an arbitrary unit, or a.u.) in these two areas based on our simulations. We found that upon the loss of EPFLs there was a moderate increase of total WUS expression near the middle apical-basal axis, but the increase was more prominent in the clv3 mutant (Fig. 2B, blue). This increase in the clv3 mutant is comparable to that observed in the erf-clv3 mutant, suggesting that the apical-basal boundary of the WUS-expressing zone is primarily controlled by CLV3. In contrast to the dominant role of CLV3 in Region 1, EPFL and CLV3 both control WUS expression in Region 2: the loss of either factor (single-mutants) resulted in a significant increase of WUS expression, and the loss of both factors produced an even more dramatic increase (Fig. 2B, Indian red). Interestingly, the sum of the effects on WUS expression from the two single-mutants was less than the effect of the erf-clv3 mutant (Fig. 2B, hatched vs. solid), suggesting a nonlinear cooperativity between EPFL and CLV3 signals. We found that this cooperativity was maintained when we perturbed the parameter values in the model (Fig. S2). Together with the observation that the loss of EPFL resulted in an increased CLV3 signal (Fig. 1E), our model further suggests that part of the lost inhibition of WUS expression in the erf mutant is compensated for by an increased level of CLV3. Together, these results show that EPFL and CLV3 signals synergistically control the lateral boundary of the WUS-expressing zone, whereas the apical-basal boundary of the WUS-expressing zone is primarily controlled by the CLV3 signal.

In addition to the synergy between EPFL and CLV3 in the regulation of WUS, we found that the spatial distributions of the cells producing EPFL and HAM (Fig. S1B and C) were critical for maintaining the expression patterns of both WUS and CLV3: the expansion of HAM production to upper layers resulted in the reduction of CLV3 expression and the lateral expression of WUS expression; the expansion of EPFL to the middle SAM reduced the overall CLV3 expression, and it dampened the WUS expression in the middle SAM while causing the lateral expansion of the WUS expression (Fig. S1B and C). These results suggest the synergy between EPFL and HAM, and the importance of their spatial distributions.

Inhibition of WUS and CLV3 expression by EPFL is important for maintaining SAM patterning

It was recently observed that EPFL inhibits the expression of both WUS and CLV3 [29]. It was, however, unclear how each of these inhibitions contributes to SAM patterning. We therefore removed these two inhibitory relationships from our model one at a time and examined the steady state patterning of the SAM (Fig. 3). We found that upon the loss of CLV3 inhibition by EPFL, there was a 30% increase of both CLV3 mRNA and CLV3 protein in the SAM. Interestingly, CLV3 was expressed around the organizing center in both the apical and lateral directions (Fig. 3A and C). This lateral expansion of CLV3 expression suggests that the inhibition of CLV3 by EPFL is critical for restricting CLV3 expression in the top layers of the middle SAM. In addition, WUS expression moderately decreased compared to the wild type (Fig. 1C) due to the expanded CLV3 expression. When we removed the inhibition of WUS by EPFL, the WUS expression region had a significant lateral expansion, and its expression in the middle region also increased compared to the wild type (Fig. 3B and C, dotted line). These analyses demonstrated the paradoxical roles of EPFL on WUS expression: it directly downregulates WUS in both the lateral and middle regions of the SAM while it indirectly upregulates the same gene by inhibiting CLV3. The direct downregulation has the dominant effect, because blocking this regulation produced a WUS pattern similar to that of the erf mutant (Fig. 3, gray line). Nonetheless, the slightly higher WUS expression in the blocked EPFL-to-WUS condition than in the erf mutant provides further evidence for the paradoxical roles of EPFL (Fig. 3, gray line). Together, these results show that EPFL inhibition of both WUS and CLV3 is critical for SAM patterning.

Fig. 3

The effects of WUS and CLV3 inhibition by EPFL. A. A simulation with the removal of CLV3 inhibition by EPFL. B. A simulation with the removal of WUS inhibition by EPFL. In both A and B, the indicated transcriptional inhibition was removed from the model by increasing the threshold to a large number (1000). Simulations were performed with the procedure as those described in Fig. 1. Steady state SAM patterns were analyzed. Wild type and erf mutant patterning are shown for comparison. C. Left: WUS and CLV3 mRNA in the middle region (Fig. 2, blue). Right: WUS and CLV3 mRNA in the lateral region (Fig. 2, Indian red). Dotted lines are visual guides for wild-type WUS mRNA. Gray lines are visual guides for erf mutant WUS mRNA. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We found that the alterations of the gene expression upon perturbations were not always intuitive. For example, blocking the EPFL-to-CLV3 inhibition resulted in a reduction of CLV3 mRNA in the middle of the SAM (Fig. 3C), suggesting the importance of the negative feedback loop between CLV3 and WUS. In the next section, we describe our analyses of the interconnected feedback loops in the SAM.

Paradoxical feedbacks between WUS and CLV3 maintain robust patterning in apical-basal axis of SAM

We noticed that the SAM regulatory network contains two feedback loops between WUS and CLV3 (Fig. 1A and Fig. 4A). In the first network motif, WUS activates CLV3 expression and CLV3 inhibits WUS expression, forming a negative feedback loop. The second motif consists of WUS inhibition by CLV3, and CLV3 inhibition by HAM and high concentration of WUS [13], [14], [15]. This motif is a double-negative (positive) feedback loop (Fig. 4A, top panel). These two feedback loops are both interconnected (i.e. they share a regulation) and paradoxical (i.e. they are negative and positive respectively). While the first loop is a well-known negative feedback that has been proposed to be critical for maintaining the homeostasis of the SAM, the function of the interconnected and paradoxical feedbacks between WUS and CLV3 is unclear. To examine the roles of these two feedbacks on SAM patterning, we constructed two alternative models with disrupted negative and positive feedbacks respectively (NFL-KO, PFL-KO models. Fig. 4A). For NFL-KO model, we introduced a hypothetical activation signal for CLV3 in the L1 and L2 layers of the central zone. This signal avoids the trivial outcome of the complete absence of CLV3 in the SAM, and the addition of the signal gave rise to a SAM pattern similar to the wild type (Fig. 4B). For the PFL-KO model, we removed the downregulation of CLV3 by the WUS signal. In order to compare these alternative models to the basal (wild-type) model in terms of WUS expression patterning, we further adjusted the synthesis rate constants of CLV3 mRNA in the alternative models such that the overall WUS expression patterns obtained with the three models are comparable (Fig. 4B).

Fig. 4

Roles of paradoxical feedbacks between WUS and CLV3. A. Top: two motifs showing the negative and the positive feedback between WUS and CLV3. Middle: network structure of negative feedback knockout model (NFL-KO) which does not have a WUS-CLV3 negative feedback loop. Dashed line indicates a hypothetical activation signal for maintaining CLV3. Bottom: network structure of positive feedback knockout model (PFL-KO) which does not have the WUS-CLV3 positive feedback loop. B. Steady state distributions of WUS and CLV3 mRNA in SAM with circuits of wild-type, NFL-KO and PFL-KO models. C. Steady state distributions of WUS and CLV3 mRNA in SAM with circuits of wild-type, NFL-KO and PFL-KO models with 50% reduction of WUS mRNA production rate constant. D. Total WUS mRNA in SAM with respect to reduction of WUS mRNA production rate constant. E. Distributions of 17 middle SAM cells in the phase space of WUS and CLV3 mRNA with three models under normal condition. Each dot indicates a cell and in the middle region of the SAM (Fig. 2A Indian red), and its coordinates show the levels of WUS mRNA and CLV3 mRNA, respectively. Contour maps show the density of these cells in the WUS and CLV3 mRNA space. Numbers in quadrants indicate cell numbers. The cell numbers in the top left and lower right quadrants reflect the diversity of the WUSon and CLV3on cells. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

We next perturbed the expression of WUS by reducing its mRNA synthesis rate constant. This perturbation reflects possible aberrant environmental or genetic changes. Because the regulation of WUS itself was not altered in the PFL-KO or NFL-KO models, this perturbation was not biased towards any model in a trivial way. However, the wild type model was robust with respect to the perturbation in terms of total WUS production in the SAM (Fig. 4C): a 50% decrease of the rate constant resulted in less than 10% loss of the WUS in the wild type SAM (Fig. 4C, orange). In contrast, the loss of either negative or positive feedback decreased the stability of WUS production when challenged with the same perturbations. In particular, the NFL-KO model had a more prominent decrease of WUS production (>3-fold from unperturbed condition) than the PFL-KO model did (Fig. 4C, purple), suggesting that the negative feedback has a dominant role in protecting the stability of WUS production in the SAM. In earlier analysis, we showed that blocking the inhibition of CLV3 by EPFL gave rise to a WUS expression pattern comparable to the wild type (Fig. 3A), a feature found in the NFL-KO and PFL-KO models. We observed that the inhibition of CLV3 by EPFL also facilitates the maintenance of WUS expression stability, but its effect is less than that of the feedback loops (Fig. 4C, yellow). We found that when the WUS mRNA production rate constant was reduced in the wild type, the level of CLV3 expression was significantly decreased (Fig. 4D), which in turn compensated for the loss of WUS mRNA production, thereby maintaining the stability of WUS expression in the SAM via a negative feedback. This compensation effect via downregulation of CLV3 was not observed with the NFL-KO model, and was only moderate in the PFL-KO model as compared to the wild type model (Fig. 4D). We further perturbed the system with increasing WUS production rate constant, and we used additional metrics to examine the response in terms of the SAM patterning. We found that the wild-type model consistently performed equally well as or better than the alternative models (Fig. S3).

Positive feedback (or double-negative feedback) is widely used in developmental systems for maintaining tissue boundaries [35], [36]. In the wild type SAM, cells along the middle apical-basal axis are either WUSoffCLV3on or WUSonCLV3off except for L2 cells which had significant expression of both WUS and CLV3 (Fig. 4B and E). The robustness of the anticorrelated expression pattern of WUS and CLV3 in most cells of the SAM was also reflected in mutant simulations when EPFL-driven inhibition was partially or completely lost (Fig. 1E, 3, and 4B). However, in the absence of the positive feedback loop (Fig. 4A, PFL-KO model), the number of cells co-expressing WUS and CLV3 increased significantly, and the anticorrelation pattern of WUS and CLV3 expressions was lost (Fig. 4B and E).

Together, these results show that the paradoxical feedbacks between WUS and CLV3 are critical both for stabilizing WUS expression and for generating the separation between WUS-expressing and CLV3-expressing regions. Stable WUS expression may be important for maintaining the size of the SAM, while separation between WUS-expressing and CLV3-expressing cells may be important for maintaining diverse cellular properties, such as proliferation rate, in different regions of the SAM [23]. The specialization of cellular phenotypes along the apical basal axis may in turn support the robust structure of the SAM.

Tradeoff in maintaining SAM patterning

Due to the lack of quantitative measurement in terms of the SAM patterning, we did not attempt to obtain a single parameter set that gives the best fit to experimental data in this study. However, the emergence of the simulated SAM patterns that are qualitatively consistent with many experimental observations allows us to interrogate the influence of the parameter values on the SAM patterning in a collective manner. In the previous sections we discussed three key features of the SAM regulatory network: 1) the network uses two signals to fully repress WUS expression in the lateral region of the SAM; 2) the network maintains the stability of WUS production in the SAM using two feedback loops; and 3) the network generates a distinct populations containing WUSon cells and CLV3on cells in the middle region of the SAM. If these three features are critical for plant physiology, how would the kinetic rate constants of the system be optimized to achieve them? To gain insight into this question, we used three metrics describing three performance objectives: lowering WUS expression in the lateral region, lowering the variability of WUS production in the presence of perturbations, and increasing the number of WUSon cells and CLV3on cells in the middle region (Fig. 5A). For easier interpretation, all three metrics were designed such that lower scores mean ‘better’ performance and higher scores mean ‘poorer’ performance (see Methods, Fig. 5B top). We perturbed 23 model parameters by decreasing each of them 2-fold, and then increasing each 2-fold, scored the performance of the perturbed models using the three performance metrics, and then normalized them to the scores obtained from wild type model (Fig. 5B) by subtracting the wild-type scores from the scores of the perturbed models (see Methods). Among the three metrics, the scores of the diversity metric (Metric III) were most robust with respect to the perturbations among the three metrics, possibility because it is based on a discrete measurement in terms of cell numbers rather than a continuous one. Changes in each parameter had some influence on at least one performance score, and changes in 13 out of the 23 parameters had one performance score better than the basal one. However, none of the changes of individual parameters gave rise to better performance in two or three metrics (Fig. 5, heatmap, Fig. S4).

Fig. 5

Influence of parametric changes on three performance objectives of SAM patterning. A. Illustrations of the three performance metrics of SAM patterning. See Methods for details of these metrics. B. The performance scores upon perturbations of 23 parameters. Each parameter was decreased and increased by 2-fold and simulation was performed in the same procedure as that with the wild-type (basal) model. Steady state SAM patterning was scored based on the three metrics, and the scores were compared with those obtained with basal model. Blue: perturbed model has better score than basal model does. Red: perturbed model has poorer score than basal model does. White: perturbed model and basal model have the same score. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

This analysis suggests a possible tradeoff in optimizing the kinetics of the SAM network: achieving better performance in one metric is likely accompanied with jeopardized performance in at least one other metric. However, this phenomenon might be simply due to the fact the very few perturbations of individual parameters had improved performance in any metric, so the absence of better performance in two or three metrics may just occur by chance. To further examine the relationship among these three performance objectives, we generated 104 perturbed models with all of their parameter values deviated from the basal set. These parameter values were randomly selected from intervals bounded by values with 50% changes from the basal set. Among these 104 randomly perturbed models, 1964 of them had improved performance score (less than 0, or better-than-basal) for Metric I, 202 of them had improved score for Metric II and 423 of them had improved score for Metric III (Fig. 6A). However, only 2 of them had improved scores for all three metrics. We performed a permutation test and confirmed that this number is significantly low (Fig. 6B, upper panel, p < 10−4). Furthermore, when a model had a better-than-basal score in one or more metric, the probability that it also had a lower-than-basal score in one or more metric was 99.8%, a number that is significantly higher than expected (Fig. 6B, lower panel, p = 0.006). These results show that there exists a tradeoff in optimizing the three performance objectives of SAM patterning. We then asked which pairs of metrics have a significant tradeoff problem when multiple kinetic rate constants are varied, and we found that if a model has a better score for Metric I or Metric II (lateral inhibition and stability), there is a significantly high probability that it has a poorer score for Metric III (divergence), suggesting that there is a pairwise tradeoff between optimizing for WUS-CLV3 diversity and other features (Fig. 6C). Other pairs of metrics did not show significantly high probability of having such opposite trends of performance change (Fig. 6C). Overall, these results suggest that there may exist a tradeoff among the three objectives of the SAM patterning when the system varies its rate constants to achieve those goals simultaneously.

Fig. 6

Tradeoffs in optimizing parameters to gain better performance of SAM patterning. 10,000 models were generated with random parameter values around those of the basal (wild-type) model. These values were chosen from a uniform distribution bounded by ±50% from the corresponding basal values in the wild-type model. Each model was simulated in the same way as with the basal model, and then the steady state SAM patterning was scored with the three metrics mentioned in Fig. 5 (see Methods for details). Positive scores mean poorer performance than wild-type model. Negative scores mean better performance than wild-type model. A. Score distributions in phase space of all pairs of metrics. Red numbers show cases in quadrants, e.g. numbers for the lower left quadrants show the numbers of models with improved scores from basal model in two metrics. Models with scores on the x or y axis were not included in red numbers. Spearman correlation coefficient and its p-value are shown in each panel. B. The scores of the perturbed models were arranged in a 10,000 × 3 matrix. Each of the three columns were randomly permutated 10,000 times independently, and the resulting 10,000 matrices were analyzed based on statistics of their rows. Blue: histogram of scores or probabilities with permuted scores. Red: Observed scores or probabilities with the perturbed models. Top panel shows the cases with three better (improved) scores (there are two observed cases out of the 10,000). Lower panels show the probability of obtaining at least one poorer score under the condition of obtaining at least one better score. C. Pairwise conditional probabilities of obtaining opposite performance in two metrics. The right-tail p-value is >0.05 unless otherwise indicated. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Discussion

Two-axis control of WUS expression in the SAM

Patterning of WUS and CLV3 expression in the SAM is considered to be a crucial system for stem cell maintenance in plants. In this study, we built a multicellular model describing the key gene regulatory network controlling SAM patterning. In particular, we considered signals that regulate WUS expression from both the middle SAM and peripheral regions. This differs from previous models which focused on the middle SAM region around the apical-basal axis [11], [13], [16], [17], [18], [22], [23], [24], [25]. Our model shows how the system combines the peripheral signal EPFL with the CLV3 signal originating from the top central SAM to restrict WUS expression to the organizing center in the meristem rib. In addition, this peripheral signal synergizes with the WUS and HAM signals from the meristem rib to restrict CLV3 expression to the top central SAM. Notably, this localization of CLV3 expression can be achieved without assuming pocket-like expression regions of any inhibitory factor [13], [22]. The prominent roles of EPFL in regulating SAM patterning by modulating WUS and CLV3 expression is consistent with its key roles in plant development that were identified previously [28], [29], [30]. The model reveals a remarkable cooperativity of the peripheral and middle signals for shaping the patterning of WUS and CLV3 expression. This cooperativity is particularly manifested in our observation that WUS expression in the lateral SAM region is controlled by both CLV3 and EPFL, each of which has partial but significant impact on WUS inhibition. Our results suggest that stem cell maintenance in Arabidopsis requires a highly regulated crosstalk between the middle and peripheral regions of the SAM rather than signals in the middle apical-basal axis alone.

Patterning and size regulation of the SAM and model limitations

Our modeling study focuses on gene expression patterning, which has a major influence on the size of the SAM. Because of the central role of WUS in maintaining the stem cell population and its positive correlation with cell proliferation rates, it is reasonable to use changes in the amount of WUS in the SAM as an indicator for changes in SAM size. For example, upregulation of WUS that was observed in erf, clv3 and other mutants is expected to correlate with an increase of stem cell populations, which in turn gives rise to an increase in SAM size. This is consistent with previous experimental observations under such conditions [29], [32], [37], [38], [39]. However, the interplay between cell proliferation and gene expression patterning is bidirectional, dynamic, and more complex than this simple inference assumes. For example, lateral expansion of the SAM during development could in theory reduce the effect of EPFL on the middle region unless the production rate of EPFL scales with the expansion. Understanding these forms of interplay will be important for quantitative characterization of SAM patterning and size control. Furthermore, it has been shown that the size and the shape of the SAM undergo significant changes in processes such as ontogeny [40], suggesting the complexity of the size-patterning interplay. Finally, cytokinin and auxin hormones, transcription factor SHOOTMERISTEMLESS (STM), microRNAs and other CLE peptides may influence meristem patterning directly, by regulating the rate of cell proliferation, or by interacting with the WUS-CLV3 loop [41], [42], [43], [44], [45], [46], [47], [48]. Ultimately, these regulators also need to be incorporated into the model of SAM regulation once their role in the SAM and the interplay with the WUS-CLV3 loop are clearly defined. Nonetheless, our work provides critical insights into the middle-peripheral interplay of SAM patterning, which can serve as a foundation for future development of more complex and more realistic models of SAM growth.

The roles of paradoxical feedbacks

It has been long considered that the negative feedback between WUS and CLV3 plays an essential role in maintaining the stem cell population in the SAM [9], [10]. The self-limiting property of negative feedback ensures the stable expression of both genes. However, recent data suggest that the interactions between WUS and CLV3 might be more complex than a negative feedback loop: WUS and HAM from the meristem rib can inhibit CLV3, which is prevented from being expressed below the subepidermal layer [13], [14], [15], [49]. As such, WUS and CLV3 are involved in a negative feedback loop as well as a positive (or double-negative) feedback loop. The latter network motif is known for its function in generating switch-like behaviors and formation of tissue boundaries during development [50], [51]. Consistent with these features, we found that this interconnected, paradoxical feedback WUS-CLV3 network supports the stability of WUS expression against alteration of kinetic rates, and it facilitates the diversification of SAM cells, in particular cells in the middle region, into WUSon and CLV3on cells. Although the two types of cells do not have a clear boundary along the apical-basal axis at single cell-layer resolution, the formation of a heterogeneous population with stable phenotypic composition is consistent with established theories about positive feedback loops [52].

Since the homeostatic function of a negative feedback loop does not necessarily depend on the spatial separation of two types of molecules, why does the SAM need to separate WUSon and CLV3on cells? One plausible reason is that the stem cell transcriptional program depends on a high concentration of CLV3 and a low concentration of WUS. Another possibility is that the differential cell proliferation rates along the apical-basal axis [23], which can be governed by heterogeneous expression patterns of WUS and CLV3, may play essential roles in maintaining the structure of the SAM. Future work is warranted to connect these molecular and cellular features of SAM patterning to the physiology of plants.

Tradeoffs in optimizing the SAM gene regulatory network

Tradeoffs in designing biological circuits have been extensively studied previously; e.g., the tradeoff between noise attenuation and speed of response, and that between maximizing information content in individual cells and cell populations [53], [54]. However, the roles of tradeoffs in complex traits such as tissue patterning remain elusive. Although the exact physiological function of SAM patterning is only partially understood, it is clear that the structure of the SAM gene regulatory network and its underlying biochemical properties serve multiple purposes, because diverse alterations of SAM patterning are often associated with abnormal development ranging from organ formation failure to nonoptimal size of organisms [6], [7], [28], [33]. Our modeling study revealed that the kinetic rates in the SAM gene regulatory network have a general tradeoff among achieving multiple ‘desired’ features of patterning. For example, parameter sets that can achieve more diverse populations of WUSon and CLV3on cells are likely associated with less stability in terms of WUS expression, or less regulation of WUS expression in the lateral region. These nontrivial tradeoffs may act as a complex selection pressure similar to multi-objective optimization problems, and they can play crucial roles in shaping the structure and kinetic rates of the SAM regulatory network through evolution. Furthermore, while there are a wide variety of SAM regulations in plants and a single optimal set of biochemical kinetic rate constants may not exist even within a single species, the tradeoff among the goals may serve as a general principle in shaping these kinetic rate constants.

Overall, our model includes a few key elements for controlling SAM patterning that were not considered previously. It captures many patterning phenotypes observed under normal conditions and with genetic perturbations. It offers insights into the interplay between peripheral and middle SAM signals, the intricate feedback regulations, and the principles governing the network design of this system.

Resource availability

Computer code for reproducing all key results and figures is available at: https://github.com/lfsc507/sam

Methods

Construction of mathematical model

To model the spatiotemporal dynamics of gene regulation in the shoot apical meristem (SAM), we considered 51 cells that are organized in a dome-like structure. We estimated this number of cells in the SAM from Chen et al. [38], and a similar number was used in a recent SAM model [13]. We described this 2D cellular network with 51 points within a half-circle with a radius of 25 [38]. In the model, the diffusive EPFL ligands are synthesized in the peripheral regions distant from the middle apical-basal axis and inhibit the expression of both WUS and CLV3 through binding to their receptors which are assumed to be broadly expressed in the SAM [28], [29]. CLV3 is a diffusive peptide, and WUS is a transcription factor that moves across cells [11], [55], [56]. In addition to the WUS-CLV3 negative feedback and these lateral regulators, our model describes a HAIRY MERISTEM (HAM) signal that originates from the rib zone and inhibits CLV3 expression in the organizing center [13]. In this model, we do not assume that HAM can move across cells due to the lack of experimental evidence. The spatial distribution of HAM expression is established by a diffusive microRNA that is not considered in the model explicitly [46]. It has been shown that HAM and WUS control gene expression synergistically [49], and that high concentration of WUS may also contribute to CLV3 downregulation [14], [15]. We therefore assumed that the inhibition of CLV3 expression by HAM signal depends on WUS. As such, WUS has a paradoxical role (both activation and inhibition) in regulating CLV3 in the presence of HAM. Finally, we considered a CLV3 independent positive feedback involving WUS. This feedback may be supported by a WUS-cytokinin mutual activation loop: it was previously shown that cytokinin activates WUS expression [23], [39], [57], whereas WUS derepresses cytokinin signal by inhibiting Type A ARABIDOPSIS RESPONSE REGULATOR (ARR) genes which act as inhibitors of cytokinin [14], [58], [59]. In addition, the WUS autoactivation loop may be supported by other factors [27]. Based on these assumptions, dynamics of six interacting species representing concentrations of regulatory molecules is described with nonlinear ordinary differential equations (ODEs) in each cell (point) of the model (additional spatial constraints are shown in Fig. 1B):

here, state variables , , , , and represent the concentrations (or strengths) of WUS mRNA, WUS protein, CLV3 mRNA, CLV3 protein, EPFL, and HAM respectively. represents a meristem rib signal that combines both HAM and WUS, the latter of which inhibits CLV3 expression at high concentration. A full list of parameter descriptions and their numerical values is available in Table S1. Briefly, is the production rate constant of molecule ; is the degradation rate constant of molecule ; is the threshold of activation or inhibition of by . is the cooperativity of activation or inhibition of by ; is the rate constant of passive diffusion‐like transport of molecule is the Laplace operator describing gradients of concentrations, which govern passive diffusion‐like transport; has a unit of concentration per unit area. was adjusted by multiplying with a scaling factor , where represents the distance between the centers of the two cells [60]; and neighboring cells are defined as cells that are located within a radius of 10 . We neglected the subcellular geometry of the cells, their contact areas and the influence of mechanics in this study (the effected contact area for WUS transport cannot be directly inferred from total contact area of plasma membrane) [13]. The diffusion of EPFL and CLV3, and the diffusion-like symplastic transport of WUS are responsible for the intercellular communication in the model. We used Hill function to describe nonlinearity in the gene regulation. Previous models of the SAM and other complex systems have used similar nonlinear functions [17], [19], [22], [61]. is a constant for us perturb the negative feedback regulation (see next section). When a molecule is controlled by multiple factors, we assumed a multiplicative form of Hill functions (AND-gate-like) (for example, nonlinear interaction of WUS and HAM may arise from their physical interactions [49]), except for the inhibitions of WUS by CLV3 and EPFL, which were assumed to be additive (OR-gate-like). Because the absolute concentrations of these molecules have not been measured experimentally, we used an arbitrary unit (a.u.) to describe concentration (or strength) of each molecule. Once these measurements become available, one can easily scale these variables to fit to specific concentrations. We used no-flux boundary condition for the model, and this is similar to a recently published SAM model [13].

We fit the parameters to known patterning phenotypes of the SAM under normal and genetically perturbed conditions (see details of mutant models below). These phenotypes are listed in Table 1. Because only qualitative information is available from the experimental data, we performed the fitting manually. Here, we do not attempt to obtain an optimal set of parameters that give the best fit to experimental data. We instead discuss the general trends of the influence of each parameter on multiple features of the SAM patterning (see Performance Metrics). To perform a simulation for a SAM system, we solved the system of ODEs numerically using the Tellurium package [62]. The initial concentrations for all variables were set to zero. An example of the time course solution of the wild type SAM is shown in Movie S1. For all our analyses, steady state solutions (at Day 100) were used to determine the patterning of the SAM.

Simulations for mutants

For each mutant SAM model, we simulated the genetic perturbation by setting the production rate constant for the knocked-out gene(s) to zero. These parameters include for the wus mutant, for the clv3 mutant, for the erf mutant, and for the ham mutant.

To examine the roles of individual interactions in the gene regulatory network, e.g. the activation of CLV3 by EPFL, we set the parameters describing the activation/inhibition thresholds () of the interactions to 1000 a.u., which exceeds the maximum concentrations of the activators/inhibitors.

All other parameters were kept the same as those in the wild-type model. Steady state distributions of all modeled molecules were obtained with the same procedure as the simulation for the wild-type SAM.

Normalization of concentrations for visualization

All analyses of molecular abundance (concentrations) were based on raw values obtained from the simulations. However, because visualization of the SAM patterning with multiple molecules would involve using the same color scale (jet colormap) for different ranges of concentrations, we normalized the concentrations of all molecules by dividing all values by their own maximum concentration across all genotypes before visualization. Therefore, the ranges of all visualized abundance are [0, 1], where the value 1 effectively represents the maximum concentrations of individual molecules across all simulations.

Models without feedback loops

To examine the roles of feedback loops on SAM patterning, we created two alternative models, each of which has a feedback loop (negative or positive) removed from the basal (wild type) model. To ‘knockout’ a feedback loop, we first removed the regulation of CLV3 by either the WUS or HAM signal by setting threshold constants () to 1000 a.u.. The removal of these interactions generates prominent changes of SAM patterning. For example, if WUS-to-CLV3 inhibition is removed from the basal model, then CLV3 is not expressed in the SAM and WUS is highly expressed compared to the basal model. However, these changes may not reflect the fitness advantage of the feedback loop because the other kinetic rate constants can be altered to compensate for the effect of the removal of the interactions. We therefore further adjusted the parameters to mimic this compensation. For the negative feedback knockout model, we introduced an activation signal for CLV3 by letting , reducing by 15%, and reducing by 50%. With this adjustment, the WUS expression pattern became comparable to the wild type. For the positive feedback knockout model, we reduced by 85% and reduced by 50%. The WUS expression pattern with this model also became comparable to the wild type. After these changes, we compared the feedback knockout models with the basal model in terms of WUS expression patterning in similar dynamic ranges.

Performance metrics for SAM patterning

To quantify the effect of changes in parameter values on SAM patterning, we used three performance metrics. These metrics are used to describe traits that might be critical for normal plant physiology rather than fit to experimental data, which are not available in a quantitative manner in most studies of SAM patterning. We do not claim that plants optimize their biochemical rate constants to achieve these three goals in general, because there are many other objectives that plants must achieve to gain better fitness. Our focus is rather on the relationship among these three goals.

Metric I (lateral inhibition) describes the ability of the system to inhibit WUS expression in the lateral region of the SAM. Specifically, it is the total WUS mRNA at steady state in 34 cells that are closest to the peripheral boundaries, i.e.

where is the total number of cells in the lateral region (34), and is the steady state WUS mRNA concentration in cell .

Metric II (WUS stability) measures the deviations of total WUS in the SAM from unperturbed models to perturbed models when the mRNA production rate constant is reduced by 0–50%. We chose 10 levels of such reduction in the parameter in the interval [0, 0.5], and obtained the deviation score given by:where is the number of perturbations (10), n is the number of cells the SAM (51), and is the steady state level of WUS protein in cell with perturbation , and is the steady state level of WUS protein in cell with the unperturbed parameter set.

Metric III (WUS-CLV3 heterogeneity) describes the size and heterogeneity of the cell population consisting of WUSonCLV3off and WUSoffCLV3on cells in the middle region of the SAM. First, the WUS and CLV3 expression was binarized with a threshold of 0.5 a.u.. Next, we calculated the heterogeneity score with the following function:

(4)where is the total number of WUSonCLV3off cells in the middle region (17 cells closest to the central axis) of the SAM, and is the total number of WUSoffCLV3on cells in the same region.

The scores of the basal (wild type) model are 5.3, 0.48 and −8 for , and respectively. We perturbed the parameters in the basal model in two ways and then examined the performance of the perturbed models (note that this perturbation is different from that in Metric II). In the first analysis, we decreased each parameter of the model by 2-fold, and then increased it by 2-fold. These two perturbations were performed for each parameter and six scores were obtained. In the second analysis, we perturbed all parameters in the model by randomly selecting their values from the intervals [, ], where is the basal value of each parameter. These parametric changes represent the alterations of biochemical rate constants that may occur through evolution to achieve desired fitness goals.

Since we are interested in how these performance scores (, and ) change when the parameters of the model are systematically perturbed from the basal set, we further normalized the raw values of these scores with the performance score obtained from the basal parameter, i.e.where is the performance score of the perturbed model, and is the performance score of the basal model. The same normalization was used to scale the scores for all three metrics. As such, if the perturbation gives rise to a performance better than that obtained with the basal set, the normalized score is negative. If the performance is poorer than that with that obtained with the basal set, the normalized score is positive.

Permutation test

The parametric perturbations and performance scoring described in the previous section generate performance scores that can be organized in an matrix, where is the number of perturbations. In our random perturbation of parameters, the value of is 104. We are interested in whether there exists a significant tradeoff in obtaining better scores of all three metrics when multiple parameters are perturbed. We first described a possible tradeoff in a small number of cases in which all three performance scores are improved () in these perturbations. Out of the 104 parametric perturbations, only 2 of them gave rise to three negative scores. To test whether this number is significantly low, given that the total numbers of improved (better than basal) scores for the three metrics are 1964, 2020 and 423, respectively, we permutated each column of the matrix independently for 104 times. We next compared the distribution of the numbers of cases in which all three performance scores are improved in these 104 matrices with the observed number of cases (2), and we calculated the empirical p-value with the area of the left-tail of the distribution with the right bound of 2, i.e. the probability of observing 2 or less cases.

The same strategy was used to quantify the significance of other descriptions of tradeoffs. We tested whether there is a significantly high probably of obtaining one or more poorer-than-basal score when an improved score (in any metric) is obtained with a parameter set. This analysis was then performed for each pair of metrics, i.e. permutation tests for six conditional probabilities were conducted. In these analyses, empirical right-tail p-values were obtained.

CRediT authorship contribution statement

Ziyi Liu: Investigation, Methodology, Visualization, Writing - original draft, Writing - review & editing. Elena D. Shpak: Conceptualization, Supervision, Funding acquisition, Writing - original draft, Writing - review & editing. Tian Hong: Investigation, Methodology, Visualization, Conceptualization, Supervision, Funding acquisition, Writing - original draft, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.