Centering the Organizing Center in the Arabidopsis thaliana Shoot Apical Meristem by a Combination of Cytokinin Signaling and Self-Organization

Plants have the ability to continously generate new organs by maintaining populations of stem cells throught their lives. The shoot apical meristem (SAM) provides a stable environment for the maintenance of stem cells. All cells inside the SAM divide, yet boundaries and patterns are maintained. Experimental evidence indicates that patterning is independent of cell lineage, thus a dynamic self-regulatory mechanism is required. A pivotal role in the organization of the SAM is played by the WUSCHEL gene (WUS). An important question in this regard is that how WUS expression is positioned in the SAM via a cell-lineage independent signaling mechanism. In this study we demonstrate via mathematical modeling that a combination of an inhibitor of the Cytokinin (CK) receptor, Arabidopsis histidine kinase 4 (AHK4) and two morphogens originating from the top cell layer, can plausibly account for the cell lineage-independent centering of WUS expression within SAM. Furthermore, our laser ablation and microsurgical experiments support the hypothesis that patterning in SAM occurs at the level of CK reception and signaling. The model suggests that the interplay between CK signaling, WUS/CLV feedback loop and boundary signals can account for positioning of the WUS expression, and provides directions for further experimental investigation.


Introduction
All the aerial plant parts are generated by the shoot apical meristem (SAM) situated at the plant apex. The SAM is formed during embryogenesis and in dicotyledonous angiosperms, such as the model plant Arabidopsis thaliana, it contains three layers of stem cells in the three outermost cell layers [1,2]. Clonal studies indicate that each layer contains about three long lived stem cells [3] in the central zone (CZ), which is marked by a lower cell division rate.
Daughter cells of the stem cells that stay in the CZ replenish the stem cell pool, whereas daughter cells that are placed towards the peripheral zone (PZ), which is marked by a higher cell division rate, enter differentiation and form organ primordia. The shape and the domain structure of the SAM are kept unchanged, although all cells continuously divide and differentiating stem cell daughters leave the meristem. Cell tracking and ablation experiments demonstrate that the fate of each cell is determined by its current position and not by lineage specific heritage, highlighting the importance of cell-cell communication [4]. Due to its changing cellular context, pattern formation of the shoot meristem does not rely on a stable point of reference, but rather occurs in a self organized manner [1].
Genetic studies mainly in Arabidopsis reveal that the WUSCHEL (WUS) and CLAVATA3 (CLV3) feedback loop is a pivotal regulator of stem cell number [1,2,5,6]. A small cell group underneath the stem cells named organizing center (OC) expresses the transcription factor WUS that maintains the stem cell in two ways. First, WUS protein moves into the stem cells, presumably through intercellular plasmatic bridges, called plasmodesmata [7]. In the stem cells, WUS directly binds to the promoter of CLV3 and promotes transcription, in addition to maintaining pluripotency through a yet unidentified mechanism [1]. CLV3 encode a small extracellular signal peptide that binds to receptor kinase complexes, including CLV1, and triggers an intracellular signal cascade that downregulates WUS transcription [1,8]. This negative feedback loop between OC and stem cells provides a mechanistic framework to keep the number of stem cells constant [1], see Fig 1A. Second, in the OC cells, WUS directly represses transcription of ARABIDOPSIS RESPONSE REGULATOR 7 (ARR7) and 15 (ARR15) genes [9], which encodes intracellular inhibitors of response to the plant hormone Cytokinin (CK), thereby promoting cellular CK response [10]. Hence, the question of how WUS expressions is centered and restricted within the SAM, becomes a key question in studying the stem cell homeostasis in the SAM. Several lines of evidence further indicate that CK is an important factor in shoot meristem regulation: first, mutants deficient in CK biosynthesis, reception, or overexpressing CK degrading enzymes, have a reduced SAM size [2,11,12]. Second, the CK receptor, Arabidopsis histidine kinase 4 (AHK4), is expressed in the meristem center, overlapping with the OC in its distal part [13]. The receptor is involved in the upregulation of WUS expression via exogenously supplied CK at relatively high levels, and it has been assumed to confer WUS regulation also at endogenous CK levels. CK response, measured by the reporter pTCS, peaks at the OC [14] in agreement with WUS enhancing CK response in these cells. Based on the expression pattern of the transcription factor SHOOTMERISTEMLESS (STM) that promotes expression of the CK synthesis gene ISOPENTENYL TRANSFERASE 7 (IPT7), CK is probably produced broadly throughout the shoot meristem, although direct evidence is still missing [15]. Furthermore, immunodetection of CKs suggest a rather broad and uniform distribution throughout the shoot meristem in Sinapsis alba [16]. Recent findings in rice indicate that activation of CK (clipping of a ribose residue) by the LONELY GUY (LOG) enzyme is confined to the 2-3 outermost cell layers of the shoot meristem, and it has been discussed whether active CK is locally produced in the shoot meristem and moves from the top downwards [2,17]. In Arabidopsis, there are eight LOG homologs. One of them, LOG4 is specifically expressed in the L1 layer, but the expression patterns of the other LOG genes are unknown and at least some of the other LOG genes seem to be also expressed in the shoot meristem [13,18,19].
peaks of concentration. One theoretical approach that was successfully applied to explain selfregulated pattern formation in developmental biology is the reaction-diffusion scheme first introduced by Turing in 1952 [20], and has since been applied to various fields of developmental biology [21][22][23]. Most of the applications of the reaction-diffusion scheme in pattern formation in biology have been in the form of activator-inhibitor systems. In its simplest form an activator-inhibitor system consists of two interacting diffusing molecule [24]. Modeling the self-organized pattern formation in the SAM has been subjected to various modeling approaches among which, activator-inhibitor models have been the main approach. Jönsson et al. [25] were the first to model the stem cell regulation in the SAM using an activator-inhibitor model. This pioneering work demonstrated that an activator-inhibitor based system can account for the observed expression of WUS in the SAM. Hohm et al. (2010), developed the first model that incorporates complete feedback between CLV3 and WUS. This model not only reproduced the expression patterns of WUS and CLV3 observed in the wildtype SAM but also some mutants and gene up and down-regulation phenotypes, further demonstrating the capability of activator-inhibitor models in accounting for SAM organization [26]. In [27] Fijuta et al implement an activator-inhibitor-based model of WUS/CLV3 regulation in a growing and dividing cellular template. Their work presents a model that is stable against perturbations caused by cellular growth and division., albeit lack of data has led to various assumptions, The activator-inhibitor based models can account for some fundamental aspects of stem cell regulation within the SAM. These models, like other spatial models of cellular development, have restrictions regarding the level of detail and the scope of the model. Often it is unavoidable to consider the input of other processes as pre-patterns or hypothetical components. Despite these limitation these models have been successful in providing an integrated view of the available experimental data regarding SAM. patterning. The hypothetical components of these models point out gaps in our biological knowledge that need to be addressed in order to obtain a mechanistic understanding of the SAM stem cell regulatory mechanism.
The aforementioned models focus on how the WUS expression pattern can emerge from the interaction of network components within the SAM. The computational models of SAM organization however, have not been limited to self-organizing systems, other models have focused on investigating the interplay between gene expression patterns rather than self-organized pattern formation [13,14,28]. These models focus on the experimentally demonstrated interactions between the WUS/CLV3 patterns and CK signaling/perception network [13,14,28]. For instance, Yadav et al. [28] investigate a model that relates CK perception via AHK4 receptor to pattern formation in the SAM. In this model CK is induced by an AHK4/CK signal, which is produced at the center of OC. The expression zones (i.e. binary expression templates) of WUS, CLV3, and KAN1 are restricted to CZ, OC, and PZ, respectively. Given these inputs the model can robustly establish the spatial patterns of WUS, CLV3 and KANDI1. Furthermore, these patterns can withstand perturbations caused by cell growth and division. In the aforementioned work the localized expression of AHK4 at the center of the OC is fundamental for correct patterning of WUS. This group of experimental and computational works, consistently propose that the patterning of OC takes place at the level of CK reception and signaling. Consequently, this implies that the self-organizing properties of the OC, can arise from the underlying CK signaling/perception pattern.

Aim of this study
The capability of activator-inhibitor networks in accounting for SAM patterning has been already demonstrated. In general, the experimental identification of network components has been a major challenge in application of reaction-diffusion models in biology. Particularly in the plant field, it remains a major challenge to demonstrate the existence of reaction diffusion networks experimentally.
In our context this highlights the importance of motivating the pre-patterns of a model by known biological knowledge as much as possible; when pre-patterns are abstract and cannot be directly linked to the known biological mechanisms, the task of experimental identification of network components is complicated. In contrast, when these assumptions are motivated by experimental observations, they can be more readily investigated via experimentation. As discussed earlier, theoretical and experimental data point towards CK signaling and perception as a fundamental factor in patterning of the SAM [13,14,28]. Here we aim to expand upon the current state of research and avoid incorporation of abstract assumptions in our model, by utilizing the available data as much as possible. Our model links WUS/CLV3 feedback loop to an activator-inhibitor system based on CK signaling. We demonstrate that these components function together to position WUS expression at the OC.

The Model
In order to investigate the apical-basal position and the lateral extension of the OC within the shoot meristem, we chose a two dimensional model of a longitudinal section. In our model, mobile signals are free to diffuse out of the SAM and into the surrounding cells (Fig 1). The system is described by a set of coupled non-linear ordinary differential equations on a discrete grid, where the grid points represent the individual cells. Hence, cells are assumed to be spatially uniform and intracellular concentration gradients are not taken into account, which is considered a reasonable simplification due to the difference in timescale between cytosolic mobility (fast) and the actual pattern formation process, i.e., gene expression (slow). Therefore, we use the term diffusion not in its actual physical but rather in an effective sense, meaning unbiased bi-directional spread of molecules between cells through special openings termed plasmodesmata or via the apoplast. Moreover, for simplicity, we assume that the dynamics of the WUS/CLV3 regulatory system arising from the assumed reaction-diffusion system are sufficiently faster than the cell division rate in the SAM. Therefore, the essential features of the patterning process of the aforementioned system can be well described by a static model that does not incorporate cell division or growth.

Facts and assumptions underlying the model
The proposed model is based on the following published results: • CK binds to the AHK4 receptor, which in turn causes phosphorylation of both type-A and type-B ARRs via arabidopsis histedine phosphotranfer proteins (AHPs) [29]. In absence of CK the receptor functions as a phosphatase. [30,31].
• Type-B ARRs are transcription factors that activate transcription of CK response genes, including type-A ARRs [32].
• Type-A inhibit type-B ARR function, the precise mechanisms has yet to be determined [10,33]. In general two hypotheses exist regarding the inhibition of type-B ARRs via type-A ARRs. Evidence suggests that type-A ARRs inhibit type-B ARRs via repression of upstream CK signaling. In addtition it has been proposed that type-B repress type-A via competition for phosphate molecules [34].
• There is feedback loop between WUS and CLV3 genes, where WUS moves from the OC into the stem cells and activates the transcription of the CVL3 gene. The CLV3 peptide is mobile and inhibits the expression of WUS [1,35].
• Expression of WUS is activated by CK signaling [13], presumably via canonical type-B ARRs effector genes. Additionally WUS represses the expression of type-A ARRs [9], thus promoting CK signaling.
In addition, our model incorporates the following assumptions: • The SAM consists of equivalent cells that have the potential to express all genes included in the model. The exception is the epidermis (L1 layer), which is assumed to be different from the rest of the cells in the SAM. This means that in our model the identity of the L1 cell layer is not determined via the proposed self-organizing mechanism.
• We hypothesize a molecule (L1 signal) that is supplied by the uppermost cell layer (L1) and diffuses downwardly establishing a gradient. The presence of this molecule is necessary for the cells to be able to respond to WUS signal by producing CLV3. Such a molecule has been identified by Knauer et al (2013), who charactrized a microRNA, miR394, that is produced at the L1 layer and is required for establishment of CLV3 expression. In our model the L1 signal is necessary for cells to be able to respond to WUS and establish stem cell identity [36].
• A diffusing inhibitor and a self-activating component are essential parts of pattern forming activator-inhibitor mechanisms [24]. Currently there is no evidence of such an inhibitor involved in SAM patterning; our trials show that several molecules within the model can be assumed to act as an inhibitor or to induce an inhibitor. For example, type-A ARRs appear as a plausible candidate for the inhibitor within the model. It is known that type-A ARRs inhibit CK signaling [10]. In the model the inhibitor is assumed to be downstream of the type-A ARRs. This is because the to fulfill the role of the inhibitor within an activator-inhibitor system, type-A ARRs have to be highly mobile signals. In absence of evidence regarding the mobility of these molecules, we assumed that the inhibitory function is conveyed via a highly mobile intermediate, factor X. Thus in our model two mechanisms exist for inhibition of type-B by type-A ARRs, via phosphate competition and via factor X.
• We assume type-B ARRs promote CK signaling via direct induction of AHK4. Experimental results presented on Arabidopsis eFP browser [37], (data from AtGenExpress project [38]) show that application of Zeatin leads to significant up-regulation of AHK4 levels. For the model, this assumed interaction constitutes the autocatalytic loop of the activator-inhibitor subnetwork.
• In our reductionist approach, we do not distinguish between mRNA and protein of the genes unless it is essential in addressing the question at hand. Considering the expression pattern of CLV3 mRNA [1] and its demonstrated inhibitory effect at the OC, it becomes apparent that CLV3 elicits a signal that travels further than it its mRNA. This is reflected in the model by distinguishing CLV3 mRNA and protein. The mRNA is assumed to be immobile while the protein is able to diffuse.
• In the model we assume that phophotransfer from the receptor to the ARRs are sufficiently fast processes compared to gene expression. Therefore the phosphotransfer is implemented using a quasi-steady-state assumption. See detail in Material and methods.
Reaction-diffusion modeling of the SAM has a history of more than a decade and the model presented in this work is inspired and motivated by earlier modeling efforts. It utilizes concepts and components (experimentally verified as well as hypothetical) put forward in earlier works In particular, factor X is a universal component of activator-inhibitor models of the SAM [25][26][27]39]. As the inhibitor in an activator-inhibitor systems, models consistently predict it to be a fast diffusing molecule with an expression pattern centered around SAM. To date evidence for a molecule that fulfills the role of such inhibitor and matches its predicted expression pattern has not emerged. Similarly the concept of L1 signal was first established by Joensson et al in [25]. In later works this was utilized as a signal defining the lateral expression of WUS [40], as well as in an apical-basal setting, as a cofactor that along with WUS is required for production of CLV3, in both two-dimensional [7], and three-dimensional [28]settings. In our model the L1 signal is essentially the same as the in latter; an apical basal signal required for CLV3 inductio in response to WUS. As already pointed out in [28], the strongest evidence for the existence of such a signal comes from the observation that in pCLV3::WUS both WUS andCLV3 are expressed in the uppermost three cell layers of the SAM [41].

Model equations
Integrating the above stated experimental observations and additional assumptions in a mathematical model, we arrive at the following system of non-dimensionalized coupled ordinary differential equations: where we defined: The subscript ij denotes the position x = x ij on the grid, where i and j refer to row and column indices respectively. L ij and Ck ij refer to the L1 signal and the CK expression profiles derived in material and methods. The function Γ is the transfer function for the two step phosphorelay from CK binding to the phosporylation of the ARRs. It comprises the kinase and phosphatase activity of the AHK4 receptor, the phosphorelay via the AHPs and the competition between the type-A and type-B ARRs for the AHPs, for the derivation see material and methods. The spatial coupling between the grid points is achieved by the diffusion operatorD (see material and methods). We use reflecting boundary conditions for the apical side. The basal and lateral boundary is not well defined; we use boundary conditions which are, for simulation purposes, equivalent with using an infinite domain for the apical-basal dimension; in the numerical simulations the grid is extended in basal and lateral directions until the concentrations decay to almost zero; this makes the boundary condition at the basal and lateral end of the grid irrelevant and provides a good approximation for the in vivo SAM. We close the domain basally and laterally using reflecting boundary conditions. All simulations, unless stated otherwise, were carried out in a cell grid where the meristemic dome is represented by three cell layers that are 9, 11 and 13 cells wide. The uppermost cell layer and the laterally outermost cells of the other two layers make up the L1 layer ( Fig 1B). We always checked that the grid is large enough to approximate an infinite domain in the described manner. In the following figures, the area of the grid that contains no information has been cropped.
Mobility of molecules in the model. In a model, the 'assigned' mobility of molecules can occur at the level of any intermediate components that are not explicitly present in the model. Moreover, in such a case, often, mRNA and protein of a gene are considered a single identity, hence in reality, the assigned mobility can occur at the level of either mRNA or the gene.
The correct patterning of the model depends on mobility of WUS, CLV3 peptide, L1 signal and factor X. in the model WUS needs to me mobile to reach L1 layer and trigger the expression of CLV3. The mobility of WUS protein has been demonstrated previously and WUS protein is detected at L1 layer [35]. In order to inhibit WUS at the OC, CLV3 is required to be mobile in the model. Similarly the intercellular movement of CLV3 peptide has been established [8]. As an inhibitor in an activator-inhibitor system, the mobility of factor X is required for model functionality. As mentioned earlier, a feasible candidate for the role of L1 signal is miRNA394, which has been shown to act as a mobile signal. [36] Model analysis We tested whether the proposed model can account for the observed patterns of CLV3 and WUS in the SAM and whether it can reproduce known experimental results, which are relevant to the patterning process. To this end, we performed numerical experiments: we examined whether the model is capable of reproducing the observed phenotypes wildtype SAM as well as various non-wildtype scenarios including mutants and overexpression lines. Because we use a mechanistic model, we can map experimental manipulations directly to the parameters of the model. Therefore the non-wildtype scenarios can be implemented by changing the model component that corresponds to the specific mutation, overexpression, etc. For instance, a knock out mutation is implemented by setting the production rate of the affected gene to zero. To simulate the ablation scenarios, the appropriate changes are applied to the wildtype system at the steady sate. Once the system reaches a steady state again, the resulting expression patterns are compared against the experimental observations.
Model subnetwork structure. The model in essence consists of two coupled subnetworks: WUS/CLV3 (Fig 2A, green shaded) and the CK signaling (Fig 2A, blue shaded). In addition boundary information is supplied by CK and the L1 signal (Fig 2A, red arrows; also see S1 Fig for the profiles). Parts of the CK signaling network correspond to components of a classical activator/ inhibitor system (Fig 2A). The AHK4/B/B p part of the network acts as an autocatalytic activator (Fig 2C), while the pathway leading from B p to X, fulfills the role of induction of the inhibitor by the activator (Fig 2D). The CK gradient confines the domain of pattern formation to the upper part of the meristem. The activator-inhibitor network is coupled to the WUS/CLV3 subnetwork via an incoherent feed-forward loop (B p /WUS/A) which specifies WUS expression by type-B ARR. The WUS/CLV3 subnetwork generates the expression domain of CLV3. Boundary information supplied by the L1 signal determines the correct orientation of the CLV3 expression in apicalbasal direction. S2 Fig contains examples of the model out put in the uncropped simulation grid.
Robust reproduction of the wildtype expression patterns of the genes included in the model. The sine qua non for the model is of course whether the observed wildtype pattern can be established and maintained. The simulated wildtype pattern is shown in Fig 3. WUS is present in a high concentration in a small region at the center of SAM, which in both lateral and apical basal directions corresponds to the observed experimental pattern [1]. In our model the lateral position of a single WUS peak is always at the center, whereas the apical basal position depends on the region defined by CK. The size of the WUS domain depends on the dynamics of the activator-inhibitor subnetwork as well as inhibition from CLV3.
We  [35,42]. In contrast, when WUS mobility is considered in the model the resulting expression extends to the L1 layer, S3B and S3D Fig. The predicted patterns of AHK4 and WUS by the model, overlap in the OC. This has been observed experimentally and reproduced by previous models of WUS/CLV3 interactions [14,43]. Furthermore the model predicts the WUS expression domain to constitute a sub-section of the broader AHK4 domain. This is in agreement with the observed distribution of WUS::DsRed-N7 in the apical half of the AHK4 receptor domain marked by AHK4::GFP, in inflorescence meristem [13]. Patterns of type-A and type-B ARRs, see Fig 3C and 3D, are comparable to the patterns reported by [14]. It should be noted that the pattern of type-A ARR expression in the model refer to the phohporylated portion of these proteins, while the relevant experimental data primarily consists of GUS reporter and transcriptional marker gene expressions [9,14,44]. This complicates the comparison of model output in terms of type-A ARR expression against experimental data. Nevertheless the model predicts that WUS expression domain and the domains associated with CK signaling (AHK4, type-A and type-B ARR expression domains), largely overlap. This is in agreement with experimental observations of ARR5 [13] and AHK4 [13,14]transcriptional reporters in the SAM, as well as with previous models of mutual interactions between CK signaling and WUS in the SAM [13,14].
CLV3 (mRNA) expression is limited to the tip of the meristem, with an expression zone that is wider at the apical end and narrows towards basal limit of CLV3 expression, Fig 3E. This is in agreement with the CLV3 (mRNA) patterns observed experimentally.
We investigated the effects of our choice of representation of meristem and L1 layer on model output. We observed that the model is not dependent on this particular representation; model simulations using several other representations of meristemic geometry and L1 layer, produce the correct output, see S4  Positioning the Organizing Center in Arabidopsis Meristem tested parameters for which the model performed in sufficient agreement with the experimental data we analyzed their sensitivity (see Material and Methods).
The model displays little sensitivity to variations of most of the parameters, demonstrating robustness within the defined parameter sub-space, S6 Fig. The model shows high sensitivity to parameters k 1 , k 6 , k 14 and d 2 , which correspond to production rate of type-B ARRs (k 1 ), phosphorelay inhibition by X (k 6 ), production of the inhibitor X (k 14 ) and the diffusion rate of WUS (d 2 ). The first three parameters (k 1 ,k 6 ,k 14 ) are essential for the correct functioning of the reaction-diffusion system and correspond to the activity of the autocatalytic loop (k 1 ), inhibitory effect of the inhibitor (k 6 ) and the production rate of the inhibitor (k 14 ). The model shows the highest sensitivity to the diffusion rate of WUS d 2 which represents the ratio of WUS diffusion to type-B ARR diffusion. The direct effect of d 2 is to influence the width of the WUS expression peak. Additionally WUS diffusion along with L1 signal determines the expression of CLV3. This double effect of the WUS diffusion rate d 2 on both WUS and CLV3 expression domains explains why it is the most sensitive parameter. Other diffusion/degradation rates in the reaction-diffusion subnetwork do not directly affect the expression pattern of CLV3.
The L1 signal has to be confined to a few cell layers. By altering the hypothetical L1 signal we can identify some properties of this signal, which are essential for the correct behavior of the model. By this, we can assess the model hypothesis regarding this directional signal in the SAM, when a candidate for such a signal is identified. For simulations of the wildtype shown in Fig 3, the L1 signal extends to only a few cell layers below the L1 layer as shown in Fig  3G. We examined the scenario where the L1 signals extend farther down the meristem. This can be achieved either by increasing the diffusion rate of the signal and/or by decreasing its degradation rate. The extension of L1 signal results in enlargement of CLV3 domain and presence of CLV3 in the OC (compare Fig 4A and 4B), which is never reported in wildtype. Therefore, the model predicts that in wildtype conditions the directional signal originating in the L1 layer, is confined to the 3-4 uppermost cell layers, which is in good agreement with the spread detected for miR394 [36]. Limitation of the CK response profile. For the proposed model it is only important that CK is limited to the upper 20-30 cell layers of the meristem, the actual process by which this is achieved is not important. As there is no experimental evidence that a diffusion or transport barrier-such as a Casperian strip-in this region exists, we analyze the consequences of the assumption that the CK profile is not limited by a physical barrier. Because no evidence for a directed transport of CK in the SAM exists, we consider the mobility of CK as a non-directional diffusion-like process. In this case the CK profile is governed by three parameters: the size of the synthesis zone n 0 , the average lifetime τ of a CK molecule and the effective diffusion rate D eff . Unfortunately, for none of these parameters estimates are available. CK profile was experimentally measured to cover the first 25 cell layers of the meristem [16], i.e., the synthesis regime does not extend beyond this. It seems reasonable to limit it further to the actual meristem [3]. From this follows that 1 n 0 < % 7. The long-distance translocation of CKs is mediated by the xylem and the phloem and has been experimentally investigated [45]. However, for this study the local short-distance mobility of CK across the plasma membrane and the cell wall is important, for which the mechanisms are not well understood [46]. The purine permease family and the equilibrative nucleoside transporter family have been proposed as candidates for CK transporters. While the first can transport free-base CKs in a proton-coupled manner the latter facilitate diffusion of nucleosides along a concentration gradient [47]. In any case, the mobility of CK in the SAM is determined by its diffusion in the cytoplasm and the transport/diffusion across the cell boundaries, where the latter is likely to be the limiting process. In order to obtain an estimate for the upper limit for the effective diffusion rate in the SAM tissue we consider the diffusion of a molecule in the cytoplasm. Based on measurements in E. coli, we find as a rough estimate of the diffusion constant in the cytoplasm D > eff % 241mm 2 s À1 [48,49]. An estimate for the lower limit can be obtained by considering the diffusion of molecules within the cell wall [47]. We obtain for CK as lower limit D < eff % 42mm 2 s À1 . The degradation of CK is catalyzed by CK oxidase/dehydrogenase [3]. It appears that degradation is a highly regulated process, which makes it difficult to say something about the rate. To date the degradation rate of CK in the SAM has not been measured. Therefore, no further information is available. However, we can use these considerations to obtain an idea about the average lifetime of a CK molecule. For a given diffusion rate D eff and a range n 0 of the synthesis zone the average degradation rate or life-time τ of CK follows from limiting the profile to the upper 25 cell layers (see material and methods). The resulting τ as a function of D eff and n 0 is shown in Fig 5. Due to the rough estimates available for the other parameters, there is of course a range for τ, but interestingly the deliberations above, point towards a lifetime of CK of the order of a few minutes.
An interesting question in this context is, how the phloem might alter the CK profile. The phloem starts several cell layers below the meristem [50,51]. A simple physical picture for the CK transport inside the phloem is that of mass transport due to the bulk motion of a fluid (advection). The results of the analysis described in the material and methods suggests that advection via phloem, at a distance of roughly 220μm, hardly affects the CK profile in the meristemic zone, which is rather defined by the diffusion length scale. However, in young plants the situation is quite different. Unless the diffusion length is modified during the growth process, the meristemic zone defined by the CK profile would be unrealistically large. The close proximity of the phloem to the meristemic region (from 0−10μm for a mature embryo [52] to about 80μm for a 11 days old seedling [50]) suggests that advection via the phloem can dictate the length scale of the CK profile and hence the size of the meristemic region in a young plant. In mature mgol-4 plants, the meristem is enlarged and becomes fragmented into multiple apices, each containing a separate stem cell niche. We simulated this effect by laterally extending the width of the meristem, as shown in Fig 7. In vivo, the mutant exhibits other developmental defects and the enlarged SAM does not possess a smooth and uniform edge, but forms a rather jagged and disorganized structure [53]. Our aim was to investigate wether the model, in general, is capable of generating multiple WUS centers in a larger domain. When the width of the domain is doubled, two WUS centers appear, the expression of WUS and CLV3 can be seen in Fig 7A and 7B. This multiplication of the pattern in a larger domain, is a known characteristic of reaction-diffusion systems [54], and further demonstrates the competence of a reaction-diffusion system in modeling the WUS regulation within the SAM.
The model exhibits regenerative ability of the meristem following laser ablation. Experimental observations show that after the removal of the OC and stem cell domain (SCD) in the SAM of tomato via laser ablation, two new WUS centers form at the opposite sides of the ablation [55]. Starting from a wildtype expression pattern, Fig 8A and 8B), we eliminated the WUS and CLV3 expressing cells. When the system again reaches a steady state, two new OCs and SCDs form at the either side of the ablation site in a very similar manner to the experimental observation, (Fig 8C and 8D) as well as a previous model [26]. Such regenerative ability is Positioning the Organizing Center in Arabidopsis Meristem an essential property of the SAM. The model predicts the presence of CK signaling and AHK4 at newly formed WUS expression centers after ablation. To our knowledge the presence of AHK4 expression patterns and CK activity have not been analyzed in the SAM after laser ablation. We therefore performed laser ablation experiments and tested for recovery on the level of CK signaling.

Microsurgical and laser ablation experiments
The model assumptions presented here suggest that the WUS expression in the OC is maintained via a reaction-diffusion network at the CK expression level. Laser ablation and microsurgical studies have shown that upon removing the WUS expressing cells the WUS expression is regenerated [55]. Our model further predicts that this recovery takes place at the level of CK signaling. In order to test this hypothesis, we performed microsurgical and laser ablation experiments, where we removed the cells within meristem that express the GFP under the TCS marker gene. We also performed the same experiment with GFP expressed under WUS promoter. In our initial trials, we observed that the expression of TCS fades away following the dissection of the meristem. This could be due to the lack of CK supply through the stem to the meristem. In order to compensate for the lack of supply of CK via the stem and to aid the visualization of the expression of TCS in the days following the dissection, we cultured the meristems after dissection in a CK containing medium. Our results demonstrate that the CK signaling domain within the meristem regenerates within 1-2 days following microsurgical ablation as shown in Fig 9G-9I, in a similar manner and time-frame as WUS expression, see   that contribute towards regulation of the SAM. Despite these findings, and several modeling efforts, it is still unclear how the WUS expression domain is restricted and centered within the SAM. We argue that the patterning and regulation of WUS within the SAM cannot be well understood without addressing the cell lineage-independent nature of it. The current knowledge, puts CK forward as a major factor in positioning and patterning of the SAM. We developed a model strictly based on the known mechanism of CK reception and signaling via the AHK4 receptor. Our experimental results show that, upon laser ablation of OC and SCD, the meristem as a whole is able to regenerate the WUS expressing cells as well as stem cells. In addition we demonstrate that the CK signaling domain within the SAM shows similar Positioning the Organizing Center in Arabidopsis Meristem regenerative capabilities. Considering the time-frame of the recovery of WUS and TCS expression after ablation and the current understanding of the role of CK signaling in regulating the SAM activity, the experimental results suggest that the CK signaling could be the basis of the regenerative ability of the SAM as a whole. The time-frame of the regeneration of CLV3 compared to WUS and TCS, demonstrated that the recovery of CK signaling and subsequently the WUS expression is sufficient for the recovery of stem cell population.
Concepts such as L1 and factor X have been consistent features of activator-inhibitor-based modeling of SAM and their role in the work presented here is in principle same as earlier works. The main contribution of this work lies in the observation that some know components of CK signaling network have the capability of functioning as an activator-inhibitor system. By incorporating our assumptions of a diffusing inhibitor, a feedback loop involving the type-B ARRs and AHK4 and L1 signal, we demonstrate the potential of the CK signaling network in generating patterns within the SAM, in close agreement with the experimental observations. In addition our experimental results suggest that specification of OC including its self-organizing properties can arise, at least in part, via CK signaling. If type-A ARRs are assumed to be highly mobile, the intermediate factor X is not required and the model can function with type-A ARRs directly inhibiting the phosphorylation of type-B ARRs. We explicitly tested this scenario and observed, with adjustment of parameters the model is capable of producing the same patterns. In case the type-A ARRs do not fulfill the requirements, the proposed intermediate factor X is necessary. To our knowledge the studies focused on genetic regulation of the SAM do not put forward a candidate for factor X. While inhibiting type-B ARRs, factor X is predicted to have an overlapping maxima and expression domain with type-B ARRs. At first sight, AHP6, a well-known inhibitor of CK signaling [43] appears as a likely candidate. However its expression pattern does not match the predicted pattern for factor X [56,57]. In the model factor X has a higher diffusion rate than other molecules, including the small CLV3 peptide. This hints at factor X being an even smaller molecule, perhaps a micro-RNA (miRNA). miRNA165/166, are mobile signaling molecules that suppress CK signaling via inhibition of CK production [58,59]. However the expression pattern of these molecules is very different from the predicted expression of factor X [60]. It remains to be seen whether the new and active area of research on the role of miRNAs in plant genetic regulation would identify an experimental counterpart for this hypothetical molecule.
The inhibition of upstream CK signaling via type-A ARRs is essential for correct model output. This is because the aforementioned interaction constitutes a part of the core activatorinhibitor motif. In contrast the type-A ARR inhibiton of type-Bs via competition for phosphate is not necessary for correct model functionality, see S7 Fig. Our simulations show that the model can produce the corrects output in presence and absence of phosphate competition. While verification of the mechanism of type-A ARRs inhibitory effects are out of the scope of the model, the results suggest that the upstream inhibition is the main mechanism in meristem patterning.
One important ingredient of the model is the observation that the CK concentration profile is limited to the upper 20-25 cell layers. The exact cause of this is not important for the model to work, but as there is no evidence of a diffusion barrier for CK, we explored the consequences of a diffusion-like transport of CK within the tissue. Because the determining parameters for the CK profile are unknown, we cannot limit the synthesis regime of CK, besides the plausible assumption that CK synthesis is confined to the SAM. However, based on these considerations and using estimates for the effective diffusion rate of CK we conclude that the average lifetime of a CK molecules within the SAM is of the order of a few minutes. A further consequence of the model is that the size of the meristem might be determined by two distinct physical mechanisms. The model suggests that in adult plants the size of the SAM is governed by the length scale of CK diffusion, while in young plants it is determined by the distance from the L1 to the phloem.
We have shown that a combination of boundary driven patterns-CK and L1 signal acting as morphogenes-and a reaction diffusion system including AHK4 and its assumed inhibitor can account for a variety of observed phenomena regarding the SAM. The reorganization of the SAM after laser ablation and the appearance of multiple OCs upon lateral extension of the SAM closely resemble the properties of a reaction diffusion system. Our results suggest that both a short and a long range morphogene are required for establishment and regulation of WUS and CLV3 expression patterns within the SAM. The long ranged (on the scale of the SAM) CK confines the WUS expression peak to the SAM. The short ranged L1 signal is required in order to restrict the expression of CLV3. We show that L1 signal, originating from the L1 layer and diffusing downwards, can adequately explain the induction of CLV3 expression in a specific location at the tip of the SAM. The model predicts that the signal does not diffuse past the first few cell layers beneath the L1. The recently characterized miR394, produced at the L1 layer and necessary for establishment and regulation of stem cells by WUS, provides a suitable candidate for the role of L1 signal in the SAM. The proposed minimal model focuses on specific aspects in order to understand the concepts and is not expected to capture the complex biological system in its full detail. This is specifically true for redundancy, which is a common feature of many biological systems. A survey of literature reveals a high degree of redundancy within the CK sub-network. There are several types of CKs in plants. Many ARRs have similar expression patterns and are thought to be at least partially redundant [10]. Single and even multiple mutants deficient in CK biosynthesis do not show significant SAM phenotypes. The same is true for many type-A and type-B ARR genes. Therefore, single mutant phenotypes predicted by the reductionist model presented here, cannot be expected to correspond to the observed single mutant phenotypes. However, the model is expected to exhibit systemic behavior that could be used to assess the hypothesis under study. Furthermore, the model makes specific predictions that can be utilized to design experiments to test the model hypothesis and to further clarify underlying mechanism of gene expression patterning within the SAM. In summary, we show that regulation of the CK receptor AHK4, through a reaction-diffusion mechanism can plausibly account for an array of observed phenomena regarding WUS patterning and thus providing one possible answer to the question of how the organizing center is centered.

Plant material and shoot meristem culture
Plants and cultured apices were grown under the long photoperiod (16 h light). The following lines of Arabidopsis have been described previously: Arabidopsis TCS-GFP containing an enhanced version of the published construct [61] and WUS-GFP [62] are in the Col-0 background. CLV3-GFP [19], is in the Landsberg erecta (Ler) background. For in vitro Arabidopsis shoot meristem culture, inflorescence meristems of Arabidopsis plants were dissected and transferred to MS medium containing 0.7% phytagel. For CK treatment, benzyladenine (Sigma-Aldrich) was added to the medium at final concentration of 500 μM.
In laser ablation experiments, a total of 6 WUS-GFP seedling were ablated and all subsequently recovered. Out of the 9 TCS-GFP seedling that were ablated and subsequently placed in BA medium 7 recovered. All of 5 CLV3-GFP seedling that were ablated subsequently recovered. In microsurgical experiments 2/2, 2/2 and 3/3 recovered for TCS-GFP, WUS-GFP and CLV3-GFP respectively.

Microscopy and laser ablation
Confocal analysis was carried out using a Leica upright confocal laser-scanning microscope (Leica TCS SP5) with long-working distance water immersion objectives. The cell wall was stained with 0.2% propidium iodide (PI; Sigma-Aldrich) for 1min. Apices were observed in the 3% agarose medium. Following laser settings are used for the observation; GFP (Argon laser, excitation 488nm, emission, 500-530nm), PI-staining (Argon laser, excitation 488nm, emission 600-700nm). Poking was carried out by fine tungsten needles (World Precision Instruments). Laser ablation was carried out by diode laser at the wavelength of 405 nm. Target cells were chosen using the Leica bleach point function and submitted to UV laser irradiation (90% laser power for 25 seconds). Confocal z-stacks were 3D reconstructed by MorphographX software [63].

Simulation details
Parameters were chosen from the biological and physical relevant ranges and adjusted to maximally approximate the available data. All simulation where carried out with arbitrary initial concentrations of all the molecules in the model, within the [0,0.5] range, unless otherwise stated. Simulations were continued until the steady state was reached.

The diffusion operator
The diffusion operatorD operating on the square grid index is defined by:

L1 signal and CK profile
For both the L1 signal and CK the problem can be described as diffusion in one-dimensional semi-infinite space with finite production regime. We approach the problem by dividing the space into two section, a region where the signaling molecule is produced, and a region where the production of the signaling molecule does not occur.
For simplicity we treat space as continuous; the corresponding partial differential equation for the concentration ϕ of a diffusing molecule in steady state reads: In above equations λ is the degradation rate, Γ is the production rate of the signaling molecule in the production domain and D is the diffusion rate. Rescaling the spatial dimension with the typical length scale L for a cell in the SAM tissue, i.e.x ¼ x=L, the solution to these equations unoccupied receptors: R f = R(1+αCK) −1 , where CK denotes the CK concentration (we neglect the depletion of the free CK by binding to AHK4), R is the concentration of AHK4, and α is the corresponding inverse K d value of the binding reaction. Assuming further a surplus of AHP phosphotransfer proteins compared to the amount amount of AHK4 receptors, the abundance H p of phosphorylated AHPs is given by: H p = γR b (1+βR f ) −1 , where γ and β describe the kinase and phosphatase activity of the receptor, respectively. The phosphorylated AHP (H p ) can bind either to type-A (A) or type-B ARR (B), hence A and B compete with each other for H p . Using again the quasi steady-state assumption we find for the fraction of phosphorylated A: Here σ A and σ B are the inverse K d values for the binding reactions for A and B, respectively. Inserting the above results for H p , R b , and R f we arrive at: To obtain non-dimensional quantities we rescale all concentrations with the half-maximal concentration W1 2 of the CLV3 activation by WUS, Eq (6). Defining: where Γ is the transfer function for describing the two-step phosphorelay given in Eq (8). Following the suggestion in the literature [65] that the phosphotransfer to the type-B ARR is inhibited by phosporylated type-A ARR, we assume that a non cell-autonomous inhibitor X, which is activated by A p , inhibits the phosphorelay to type-B ARR non-competitively: Rescaling time with the inverse degradation rate of WUS, the other terms in Eqs (1)-(7) follow directly from the assumptions described in Results section.

Model evaluation
In order to evaluate the resulting simulation patterns one needs to define appropriate quantitative measures. In order to evaluate the output of the model, quantitative measures must be defined based on experimental observations of WUS and CLV3 pattern in the SAM. In evaluating the model output we are concerned with general patterning capabilities rather than reproduction of experimentally observed patterns in detail. Therefore we focus on essential features that define the existence of the correct pattern. This allows for variability in model output by not imposing a too strict of a criteria for the correct pattern. The output of the model can be assessed using the WUS and CLV3 concentration distributions. We define the following marginalized distributions:

Cost function construction
A set of reference marginalized and normalized concentrations of WUS and CLV3 can be obtained from experimental data as described in [66]; GFP intensity is used as a proxy for the concentration of WUS and CLV3 in an apical-basal cross-section of the 3D confocal stacks of SAM in pWUS-n3GFP and pCLV3-n3GFP respectively. Using the two marginalized distributions given in Eqs (11) and (12), we define two objectives as the distances between the marginalized profiles: Model parameters S1 Table lists the non-dimensionalized parameters used in simulations presented in the main text unless otherwise stated.

Parameter survey
Focusing on parameter set,p 1 , described in S1 Table, we define a hypercube O on the logarithmic scale, where we extended each parameter one order of magnitude in each direction: O ¼ P j ½p j 1 =10; 10p j 1 . Within this hypercube a total of N = 2.5×10 5 parameter sets were generated. For each parameter setp 2 fp 1 ;p 2 ; :::;p N g we calculated a score LðpÞ, as described above. The subset defined by o ¼ fp 2 fp 1 ;p 2 ; :::;p N gjLðpÞ < 0:1g was obtained. Subset ω consists of parameter sets whose outputs are consistent with experimental observations. The threshold of 0.1 allows for variation in model output, while it insures the existence of the correct patterns.

Sensitivity analysis
To characterize the effect of perturbations on model output a local sensitivity analysis was carried out. The normalized sensitivity of parameter p j i belonging to the parameter pointp i is defined by: Lðp 1 i ; Á Á Á ; p jÀ1 i ; p j i þ d; p jþ1 i ; Á Á ÁÞ À Lðp i Þ d : From the set of S j i we calculated the quartiles as shown in the box-plot S6 Fig. [For the sensitivity analysis δ = 0.01 was used.
Supporting Information S1 Fig. Boundary layer profiles. Re-scaled (ϕ/ϕ 0 , ϕ 0 = ϕ(0)) concentration profiles as given by Eqs (9) and (10)  when phosphate competition between ARRs is absent from the model. The expression profile of WUS and CLV3 in absence of phosphate competition were subtracted from the corresponding wildtype expression profile. The expression profiles were normalizes to maximum levels in each case and the absolute difference between the two was calculated. To simulate the absence of phosphate competition, parameters k 4 and k 5 are set to zero. Model output is not significantly changed in absence of the phosphate competition. (EPS) S1 Table. The Parameters used in the main text simulations. All parameters are dimensionless. (TIFF)