Putting Theory to the Test: Which Regulatory Mechanisms Can Drive Realistic Growth of a Root?

In recent years there has been a strong development of computational approaches to mechanistically understand organ growth regulation in plants. In this study, simulation methods were used to explore which regulatory mechanisms can lead to realistic output at the cell and whole organ scale and which other possibilities must be discarded as they result in cellular patterns and kinematic characteristics that are not consistent with experimental observations for the Arabidopsis thaliana primary root. To aid in this analysis, a ‘Uniform Longitudinal Strain Rule’ (ULSR) was formulated as a necessary condition for stable, unidirectional, symplastic growth. Our simulations indicate that symplastic structures are robust to differences in longitudinal strain rates along the growth axis only if these differences are small and short-lived. Whereas simple cell-autonomous regulatory rules based on counters and timers can produce stable growth, it was found that steady developmental zones and smooth transitions in cell lengths are not feasible. By introducing spatial cues into growth regulation, those inadequacies could be avoided and experimental data could be faithfully reproduced. Nevertheless, a root growth model based on previous polar auxin-transport mechanisms violates the proposed ULSR due to the presence of lateral gradients. Models with layer-specific regulation or layer-driven growth offer potential solutions. Alternatively, a model representing the known cross-talk between auxin, as the cell proliferation promoting factor, and cytokinin, as the cell differentiation promoting factor, predicts the effect of hormone-perturbations on meristem size. By down-regulating PIN-mediated transport through the transcription factor SHY2, cytokinin effectively flattens the lateral auxin gradient, at the basal boundary of the division zone, (thereby imposing the ULSR) to signal the exit of proliferation and start of elongation. This model exploration underlines the value of generating virtual root growth kinematics to dissect and understand the mechanisms controlling this biological system.


Introduction
Regulation of plant growth and development has been the subject of intensive research for over a century and this will likely continue, given the increasing need for crop production to sustain a growing world population [1]. By virtue of new experimental approaches Arabidopsis thaliana has surpassed classical crops like wheat, tobacco and maize as the main model system to study the underlying molecular processes, thereby considerably advancing the field [2]. The principal processes that determine the growth of plant organs, including the primary root, are cell division and cell growth. Detailed analyses of these processes at high spatial and temporal resolution can be obtained by means of kinematic analyses providing a link with whole organ growth measurements [3,4]. Under consistent growth conditions, typically around six days after germination the Arabidopsis seedling root has established a steady root growth rate. At the distal end of the root a stable zone of cell proliferation (division or proliferation zone: DZ or PZ, respectively) precedes a delineated zone of accelerated growth (elongation zone, EZ). Cells in the elongation zone grow anisotropically reaching 20-fold length increases in as little as six hours, by massive vacuolar expansion, before reaching their final length (in the mature zone, MZ), where root hairs in specific epidermal cell files (trichoblasts) indicate a final differentiation [5,6]. This relatively simple outline of the growing root apex originates from a more complex underlying organisation. Arranged concentrically around the longitudinal axis, the root apex contains well-defined cell layers originating from stem celllike initials organized in a tiered arrangement around a zone with low mitotic activity (the quiescent centre, QC; [7]). These initials give rise to 6 distinct cell layers in the proximal direction: stele, pericycle, endodermis and cortex (from initials proximal to the QC), epidermis and lateral root cap cells (from initials lateral to the QC), and one in the distal direction: columella cells (from initials distal to the QC). How the underlying tissue structure is established and maintained remains the subject of extensive investigation [8].
As ample evidence demonstrates, plant hormones play a central role in growth and development and they enable plants to respond to changing circumstances [9]. The wealth of molecular knowl-edge has inspired various attempts to capture aspects of those complex phenomena in mathematical models (reviewed in [10,11]). Some modelling studies have focused on root growth. The landmark study of Grieneisen et al. [12] used experimental findings on PIN auxin-exporter distributions to reproduce realistic auxin patterns in the root apex that govern cell division and growth. Laskowski and coworkers [13] extended this model to account for lateral root initiation based on auxin maxima induced by root bending. Swarup et al. [14] and Band et al. [15] used modelling approaches to link asymmetric auxin redistribution at the root tip to the gravitropic response. Mironova et al. [16] presented the so called 'reflected flow' mechanism as an alternative to the 'reverse fountain' mechanism [12] to explain auxin selforganisation, especially during the very early stages of root development. Cruz-Ramírez et al. [17] proposed models that integrate radial patterns of the transcription factors SHR and SCR with longitudinal information from an auxin gradient to control asymmetric cell divisions required for differentiation within the Arabidopsis root. With a recent spatial model of the root tip Santuari et al. [18] implicate differential endocytosis in explaining how the auxin gradient is interpreted throughout the meristem. Band et al. [19] have modelled root growth as a single file of cells with growth-driven dilution of gibberellins providing a mechanism for growth arrest toward the end of the elongation zone.
Although some of those models comprise cell growth and division, these processes are typically implemented in an arbitrary and artificial way since their relationship with hormone action has not been sufficiently elucidated. Moreover, plant cells are connected mechanically through their cell walls, which by their stiffness counteract the hydrostatic turgor pressure driving cell growth. The collective framework of growing cell walls (apoplast), imposes major constraints on organ growth, preventing neighbouring cells to become separated (precluding for instance 'cell sliding'). This type of growth, designated 'symplastic', [20], has only recently been represented in a vertex-based computational model of root growth [21]. By creating mechanical stresses, in principle the apoplast operates as an integrative function in coordinating organ growth [22][23][24]. In this study we explore the prowess of previously proposed mechanisms for regulation and coordination of growth and division in a symplastic framework by building and testing computational models within the simulation platform VirtualLeaf [25]. The simulation results are discussed in relation to classical kinematic studies of Arabidopsis thaliana.

Model construction and evaluation
Kinematic growth analyses aim to infer properties of cell growth and division from microscopic time series of the growing root, yielding detailed profiles of cell expansion (extension) and cell division (partitioning) along the principal axis of growth [3]. Among others Green [26] developed and advocated this method. Interestingly, Green also derived how proposed elemental growth models would be manifested in a kinematic framework. Based on our current concepts of growth regulation we have extended this analysis that relates changes in cell size along the root axis to relative rates of cell division and expansion. To this end we used the vertex-based plant modelling software VirtualLeaf [25] to directly simulate a selection of distinct growth models (see Methods). These models evolve a two-dimensional cellular grid that represents an axial bisection of the growing root apex ( Figure 1A and 1B). Cells are defined as polygons and cell wall segments correspond to the edges acting as linear springs. Cells and cell walls are endowed with biochemical properties represented by reaction and transport equations and logical rules. Specifying regulatory mechanisms ( = input) amounts to specifying rules that minimally determine (i) cellular growth rates (via changing target areas), (ii) cell division rates and orientation (here strictly horizontal), (iii) the transition at the border of the QC and the proximal meristem, (iv) the transition between division and elongation zone (DZ and EZ, respectively), and (v) the transition to mature (differentiated) cells. By modifying those rules and analysing the resulting virtual phenotypes in terms of microscopic and kinematic characteristics we then delineate the requirements for cell growth and division needed to produce a realistic and stably growing root. We have opted here to look at proposed elementary mechanisms and gradually increase the degree of realism compared to real data.
Cell division in plants is spatially confined to meristematic zones, which has inspired discussions on whether its regulation takes place at the supra-cellular ('organismal') level instead of the level of individual cells ('cellular') [27]. Although as often the truth probably lies somewhere in between [28][29][30], we here classify regulatory mechanisms as cell-autonomous or non-cell-autonomous [31]. In the first case cells behave as pre-programmed automata and there is no role for spatial signalling. In the second case regulation relies on spatial cues that are a function of the organisation of the symplast and that feed back into the cellular decision centres. Here we will define and evaluate the outcome of simulations of models corresponding to those two classes of mechanisms.

Regulatory mechanisms must conform to the 'Uniform Longitudinal Strain Rule'
A consequence of the symplastic growth of the root is that at a given distance from the tip all cells have the same relative expansion rate [32]. As stated by Ivanov [33], any observed difference in cell lengths between tissues must therefore reflect differences in cell proliferation (see also [26]). Inversely, any form of growth regulation that results in different elongation rates for cells at the same distance from the tip would disrupt symplastic growth ( Figure 1A). For instance, suppose all cells at the same

Author Summary
The growth of a plant root is driven by cell division and cell expansion occurring in spatially distinct developmental zones. Although these zones are in principle stable, depending on the conditions, their size and properties can be modulated. This has been meticulously described by kinematic studies, which have led to the proposal of mechanisms underpinning those observations. At the same time, much knowledge of the identities and interactions of molecules involved in these mechanisms has accumulated, in particular from the model species Arabidopsis thaliana. Here we attempt to resolve the longstanding question whether observed growth patterns can be explained by autonomous decision-making at the level of individual cells or if the aid of some external signal is required. We then ask, building on the accumulated molecular information, which minimal models can provide for stable growth while keeping sufficient flexibility to regulate growth. Therefore, we constructed computational models for different growth mechanisms operating in a virtual two-dimensional Arabidopsis root and compared their behaviour with biological experiments. The simulations provide strong indications that spatial signals are required for realistic and flexible root growth, likely orchestrated by the plant hormones auxin and cytokinin.
(vertical) position in a downward growing root have the same absolute (areal) expansion rate, irrespective of their size (Model 1, Tables 1 and S1). With inner cell files narrower than outer cell files (similar to the real root) this fixed size increment results in consistently larger relative elongation rates for the inner tissue layers leading to tissue distortion and unbalanced distribution of mechanical stresses ( Figure 1B and C). Note that the same situation would occur when adjacent files contain cells of similar width, but different lengths growing at the same absolute rates. Hence, non-uniform relative strain rates at some position along the principal growth axis eventually lead to malformations. Disruption can be manifested in different forms ranging from small local cell or tissue deformations up to changes of whole organ growth. In fact, 'disruption' may be an overly negative term as it could be argued that a carefully-coordinated breach of that principle, as for instance in the root gravitropic response, can be beneficial to the plant. Furthermore, the stated necessary condition might be too stringent since small and short-lived random perturbations are likely to yield no significant distortions since they can cancel each other out to some extend as demonstrated in the next section. The apoplast, by its ability to transmit mechanical stress, may effectively function as a buffer to those small and random perturbations. Nevertheless, it seems evident that systematic differences in strain rates will eventually result in geometrical changes to the organ structure. In order to evaluate whether growth mechanisms allow for stable root growth we therefore reformulate the previously mentioned findings of Ivanov [33] and Green [26]: If in a uni-directionally growing root at least two points at the same axial position (with respect to the growth direction) have a persistently different longitudinal strain rate (relative elongation rate) then the symplastic structure will be distorted.
'Persistent' should here be interpreted as present during a minimum time interval sufficient to produce an arbitrary distortion based on that local strain rate difference. From here on, we will refer to this formulation as the 'Uniform Longitudinal Strain Rule' (in short ULSR) to emphasise its application over a finite time interval. Although root growth of the Arabidopsis seedling is in general not steady [5] in some conditions during development this is approximately the case [34]. Here we will investigate also how steady growth can be attained and maintained. Like embryonic roots we will therefore start the models from a minimal set of precursor cells ( Figure 1B). Subsequently the models will evolve to a situation where a larger population of cells occupies the meristem, elongation and mature part of the roots. In this process the growth of the root will accelerate initially, but then ideally converges to a state where a constant size of the meristem and elongation zone is obtained and the root growth rate becomes constant. For a uni-directionally growing root the following definition is proposed.
The growing root is in a steady state if all points at each position along the growth axis have a constant strain rate.
For practical reasons in the next sections our evaluation of stable and realistic growth will be largely based on three simplified criteria: 1. The overall root growth should reach a steady state.  (Table S1) with growth defined as a constant increment to a cell's target area. Cells in inner files are narrower and therefore longitudinal strain will be higher. This leads to growth irregularity with faster growth for central files creating an imbalance in the observed strain rates and distortion of the originally horizontal cell walls across the diverse layers. + Similar lengths of meristematic and mature cells, connected by a smooth transition [5]. + An approximately two-fold size range at any position [32].
2. Since only averaged data or progress curves have been reported of cell length distributions along the primary root axis of Arabidopsis, we have generated a new curve ( Figure 2) based on plants cultured in the same growth conditions as in [5]. 3. The underlying regulation of cell division and growth should not be in conflict with the ULSR.

Cell-autonomous regulation does not produce realistic root kinematics
We start off by considering strictly cell-autonomous regulation of cell growth and division as this appears to be the most elementary form of regulation. Cell-autonomous regulation implies here that essential developmental transitions are only governed by endogenous cellular programs independent of any external signals (excluding their direct physical connection in the symplast). Whereas we know this is unlikely the case, some plant growth studies point towards regulation based on cells acting as autonomous, pre-programmed units. Typically, cell behaviour is proposed to be a function of the number of events or the time passed with respect to some reference point (a defined cellular event). These types of regulation are called counters and timers, respectively. A counter keeping track of the number of cell divisions of a cell before exiting the proliferative phase has been repeatedly proposed [34][35][36]. Timers have been used in various growth models for instance on vertebrate segmentation [37] or the cell cycle [12,38,39]. It has indeed been found that in wild type roots under diverse treatments, cells spend about the same time traversing elongation zones of very different sizes (even during accelerated growth) suggesting that cells enter the EZ with a preprogrammed duration for expansion activity [5]. In his seminal study Green [26] includes a meristem doubling model in which the  (3) specifies the transition between division and elongation zone (DZ and EZ, respectively); column (4) specifies the transition to mature (differentiated) cells (mature zone or MZ) based on timing since the release from the QC or a spatial signal at a fixed distance from the root apex; column (5) specifies whether division rate is determined via a timer or sizer mechanism; and column (6) how cellular growth rates are defined. Developmental events can be determined to happen after a fixed duration ('Timer'), a fixed number of divisions ('Counter'), a fixed cell size ('Sizer'), and a fixed distance from the root apex ('Ruler'). For Models 10-12 more complicated regulatory mechanisms are specified. In Model 4 extra random noise was added to the timer ('+ noise'). For Model 5 cell division is dependent on a timer mechanism, with a different cell cycle time for inner and outer layers ('Timer (variant)'). Model 6 uses a uniform size (area) criterion for cell division, irrespectively of cell geometry ('Sizer (uniform)'). For Model 7 the size criterion is adapted to the difference in the width between inner and outer layers, thereby acting effectively as a length criterion. In the other (non-cell autonomous) cases the cell division sizer differs between layers. It is usually assumed that DZ and EZ have a characteristic elongation rate. More details of the models can be found in Table S1, Text S1 and in the corresponding figures. doi:10.1371/journal.pcbi.1003910.t001 Figure 2. Experimental cell length distribution. Cell length distribution along the principal growth axis (distance from the QC/ quiescent centre) determined for a typical 7-day-old Arabidopsis seedling root. Experimental set-up and growth conditions were in accordance with Beemster and Baskin [5]. The data points represent epidermal cell lengths. The 'polyloc method' was used for curve fitting (full line, cf. Methods). doi:10.1371/journal.pcbi.1003910.g002 meristem consists of a number of cells that undergo one cell cycle after which the proximal half of them enters the non-dividing state.
To test whether cell-autonomous regulation can produce realistic root apex structures we constructed a set of models in which cells in contact with the columella/QC maintain their capacity to divide. Upon release from the QC, these cells are allowed to divide another two or three times, based on cells keeping track of the number of divisions or the time passed since release from the columella/QC (Models 2-7, Tables 1 and S1, Figures 3 and 4). The loss of proliferative capacity is then followed by a fixed time interval of accelerated growth. In Model 2 (Table 1 and S1) the exit from the proliferative stage is determined by a timer, with cell division synchronous among cells at the same longitudinal position. Starting from a simple grid, after a transient phase, a roughly linear increase in length is achieved in accordance with the first of the proposed criteria ( Figure 3A). The synchronized cell division results in regular stepwise increases in cell numbers ( Figure 3B). All cells in one file derive from one cell (clones) occurring in groups (of increasing multiples of 2) with the same developmental state (i.e. number of divisions and time since QC release). Consequently, discrete sets of developmentally synchronous cells with similar size and growth rate arise (see Figures 3C and 3D) leading to sharp developmental transitions. The use of a counter instead of a timer to decide when cells exit proliferation as in Model 3 (Table 1 and S1), produces essentially the same outcome ( Figure S1). Similarly to what Green described with his meristem doubling model, the size of the proliferation zone periodically varies ( Figures 3C and S2), growing exponentially and losing half of its size as a complete group of clonal cells exits proliferation and enters accelerated growth. Importantly, in vivo such patterns of synchronised behaviour have not been observed. Contrary to the third (ULSR), the second criterion is breached since cell length distributions along the zones are not smooth (contrary to [5]). Indeed, some aspects of realistic root growth can be reproduced, yet, others not based on this type of model.
Simulating perfectly synchronized cell divisions produces highly regular cellular grids dissimilar to root micrographs (compare Figure S1 with anatomical sections in for instance [7]). After adding random uniform noise (within +/225% of reference) to each individual cell cycle time a more naturally looking outcome results (Model 4, Figure 4A), with a smoother cell number increase in time ( Figure S3A). Importantly, in a spatial context the cell length distributions remain discontinuous, yet smeared out more ( Figure S3C). Logically with the partial loss of synchronicity, cells at the same position along the growth axis will not necessarily grow at the same relative rate since their number of divisions or the time since release from the QC varies. Hence, the proposed ULSR is violated sensu stricto. We have found, however, through our simulations that, to some extent, the symplast acts as a stabilising framework in that local, random mechanical perturbations are dissipated through its structure ( Figure 4A, with colouring as a function of the ratio actual area/target area, with red colours representing an outgoing or cell-expanding potential: for details see Methods). However, if even more cell cycle time noise is added or if perturbations are systematic, then growth defects eventually arise. For instance in Model 5 and Model 6 (Tables 1 and S1) a systematically shorter cell cycle time of the outer cell layers is propagated over multiple generations. Through timer or counter based rules the cells in the outer layers undergo the fast acceleration phase much earlier, leading to tissue structure defects ( Figure 4B and 4C). Cell-autonomous regulation is prone to this since a local defect is easily clonally amplified, especially if occurring at an early stage in development, without external signals to constrain the effect within spatial bounds.
Likely, robust regulation also includes well-coordinated cell division and cell growth, leading to a stable cell size distribution over time. Contrary to a simple pre-programmed cell cycle timer, which does not respond to fluctuations in growth rate, a size-based cell division mechanism ('sizer') might provide this extra stability. Size-based control mechanisms are well-established in (plant) cell cycle modelling (e.g. [40,41]. If cell division is implemented to happen at an absolute cell size (or length rather as in Model 7) cell sizes (or lengths) are kept in strict bounds in the DZ, yet synchronicity of cell division is reduced ( Figure S3B).
We conclude that despite a quasi-steady growth, in the absence of spatial cues stable developmental zones and smooth gradients in cell lengths are not feasible based on the investigated types of strict cell-autonomous regulation. It makes sense that to have cells of the same generation behave according to a smooth positional gradient requires a spatial signal, even in the presence of noise.
Non-cell-autonomous regulation provides a solid basis for realistic root kinematics Model 8 (Tables 1 and S1) was constructed to test whether including spatial signals (non-cell-autonomous regulation) can improve resemblance with experimental data. Counter and timerbased rules for exit of proliferation, start of accelerated growth and maturation were replaced by a fully independent spatial signal that marks these transitions at fixed positions from the QC (like a 'ruler'). Simulations with this concept model lead to growth, which to a good approximation can be described as steady and linear ( Figure 5; Video S1). The spatial boundaries directly limit the total mass production occurring within those predefined zones. Importantly, the simulated data of (epidermal) cell lengths along the root apex are similar to what we obtained experimentally ( Figure 2), and the fitted cell length distributions are similar to what was reported in previous kinematics studies as well [5]. The strain rates are either low (in the DZ) or high (in the EZ) changing rather abruptly going from DZ to EZ for a single simulated root (Figures S4A). Using averaged strain rates from multiple simulations (with different random seeds for mechanical equilibration) shows that the strain rate curves become more smooth and bellshaped (although still skewed) as was observed in various studies [42,5,43] ( Figure S4B). Accordingly, Beemster and Baskin [5] described sharp changes in the slope of the velocity profiles for individual root samples versus averaged data (compare their Figure 2A (including inset) to Figures S4C and S4D). This rapid change was also found by van der Weele et al. [44]. Furthermore, since no conflicts with the ULSR appear to exist, non-cellautonomous regulation can effectively fulfil all proposed criteria for realistic root growth.
From a biochemical point of view the most likely spatial signal to act as a 'ruler' would be a phytohormone or combination of phytohormones undergoing long-distance transport through the root tissue resulting in a morphogenetic gradient. Auxin is arguably the most prominent root morphogen. It has been demonstrated that a stable auxin gradient can be formed which could act as the primary signal determining growth and cell division [12]. In the simulations reported with that model the growing root tip is regularly truncated providing an auxin source at a constant distance from the root tip. Via Model 9, which is a vertex-based variant of the latter model (Table 1 and S1) we investigated whether auxin patterns effectively reach a steady state with a source at a fixed position that becomes displaced further and further from the root tip due to growth. We specifically explored how spatially regulated growth and division (as in Figure 5) affect the auxin gradient (not conversely). It followed that an important factor is the definition of the auxin source: as a  Table S1). Upon release from the QC, cells can divide for a pre-programmed time duration (A-C) or number of times (D), followed by a fixed time duration of accelerated growth. (A) Total root length versus simulation time. It takes approximately 75 hours for the first cells to start accelerating growth (indicated by '*'). Subsequent sets of quasi synchronously dividing cells then grow exponentially for a fixed time, yielding a roughly linear area increase. (B) Total cell number versus simulation time. The cell number is built up exponentially till 12 (columns)62 3 = 96 'clones' are formed. These then undergo accelerated growth. One cell cycle later the next pool of 96 cells is ready to do the same and cell numbers increase by constant steps. (C) Detailed snapshots (at 99 hours, with diamond shaped markers, and 108 hours, with circular markers) of approximate cell length along the main growth axis. Subpopulations with narrow length distributions can be seen corresponding to dividing, accelerating and mature cells. The accelerating cell population occupies an increasing area until adding to the mature zone. The 'polyloc method' was used for curve fitting (cf. Methods). (D) Simulation output of Model 2 at 99 h (Table S1). The imposed growth and division rules have resulted in a highly regular grid with distinct zones of similar cell length (division zone (DZ) and elongation zone (EZ) are indicated). Areal strain rates ('AS' as defined in Methods) are mapped on the cellular grid, showing the elongation zone as a distinct region of relatively uniform accelerated growth. doi:10.1371/journal.pcbi.1003910.g003 (constant) net import from the upper plant parts (and) or as some form of local production. In the first variant of Model 9 (in the presence of first order auxin degradation) only the total auxin level of the root slowly converges to a steady state ( Figure S5A), whereas the concentration is steadily diluted (after an initial increase) by the growing root ( Figure S5B and S5C). If we define developmental transitions with a stable spatial signal, therefore this type of auxin source in principle does not support steady root growth. Rather, it might be suitable to produce temporary responses. A different behaviour emerges with variants of Model 9 that instead use local (root-based) auxin production: either with cellular production proportional to size or with a constant production rate per cell. In both cases the total auxin level ( Figure S5D and S5G) increases proportionally to the area increase ( Figure S5E and S5H), and the auxin concentration over the total root slowly (especially for the area-based production mode) converges to a steady state (Figure  Figure 3, yet with maximally 625% noise added to cell cycle time (except first division 610%, cf. Table S1 -Model 4; see also Figures S3A and S3C), which demonstrates that the symplast can act as a stabilizing framework that dissipates 'noisy' mechanical perturbations through its structure. Colouring is according to growth potential (defined in the Methods section) as a measure for 'turgor pressure'. (B) Timer-based developmental transitions can amplify tissue distortion. Simulation output of timer mechanism with layer-specific cell cycle time (Model 5 at simulation time 90 h). Upon release from the 'QC' cells can divide 3 times. For the outer 6 cell files CCT = 900 min, the inner 6 CCT = 1080 min. Afterwards, cells undergo accelerated growth until 4440 min after 'QC'-release. Therefore the inner cells lag behind if their release started more or less synchronous as is the case here (maximally 610% noise). Synchronicity is an important factor determining stable patterning in itself. Counters are sensitive as they multiply and therefore consistently amplify differences in CCT. Cells in the centre are lagging behind in terms of growth. As simulation time passes the irregularities eventually result in severely distorted patterns. Colouring is according to growth potential (defined in the Methods section) as a measure for 'turgor pressure'. (C) Simulation of Model 6 (Table S1) at simulation time 42 h, with a sizer-based cell cycle. One sizer, imposing division at a defined absolute cell size, is used despite differences in width of cells at similar positions along the main growth axis. Outer cell files have wider cells which reach the critical size before those of inner files. Therefore they undergo a much earlier exit from proliferation starting accelerated growth earlier, resulting in cell shape/tissue distortion. Cells in the centre are therefore lagging behind in terms of growth rate. Cell-autonomous regulatory systems appear inherently sensitive to this effect. Colouring is according to growth potential, GP or 'turgor pressure' (cf. Methods). doi:10.1371/journal.pcbi.1003910.g004 S5F and S5I). Production modes can be combined as well, yet this does not change the picture with typically one mode being (or becoming) dominant (results not shown). In accordance with previous work, the two-dimensional shape of the auxin gradient depends on whether a-polar (diffusive or AUX/LAX importermediated) or polarized (PIN-mediated) transport is dominant, with in the first case a gradual auxin concentration increase towards the lower (distal) root combined with a smooth auxin peak near the 'QC' cells ( Figure S6). If PIN-transport is dominant a much more pronounced local auxin maximum occurs at the central cell lines just proximal to the 'QC' cells. The dominance of diffusive versus a-polar transport is determined in this model by the relative values of the diffusion and PIN transport kinetic constants rather than their absolute values (see for instance the similarity of Figures S6A  and S6E, and Figures S6B and S6F). With the model variant based on constant auxin import that sharp peak shifts in time, connected to the growing root apex. However, it eventually disappears by growth dilution (Figure 6A). Local production modes again show the shifting of the sharp peak, however, in those cases the sharp peak slowly but gradually becomes stable ( Figure 6B and 6C). We conclude that typical auxin-based models can produce stable gradients if auxin production scales with the growing root. Importantly, unless auxin transport is completely a-polar, there is always an additional horizontal (transversal) concentration gradient besides the longitudinal gradient ( Figures 6D, S6, S7, S8, S9). This corresponds to what can be seen in various reporter studies [18,45]. How this relates to growth regulation is investigated in the next section.
Having obtained a steady and stable spatial auxin gradient, the question remains how it is interpreted by the cells as they move through it. In Model 10 (Tables 1 and S1) cell growth and division are directly determined by the local auxin concentration, e.g. according to specific step functions (division above a minimum size and above an auxin threshold, slow growth above the same auxin threshold, fast growth above a lower threshold). In this case no steady state was obtained with severely distorted root growth ( Figure 7A). Because of the radial gradient, cells at the same position along the growth axis are instructed to consistently produce different strain rates thereby violating the ULSR. With auxin concentration determining cell division, cell size distributions are not arranged in parallel along the longitudinal axis but follow the peaks and valleys in the auxin landscape (to some extent this has been observed, (e.g. [46]). When tuning the auxin transport parameters to a more dominant diffusion regime the transversal gradient remains present still precluding stable growth, although the distortions of the cellular grid are less pronounced. Some tissue locations (for instance the vascular layers near the QC) have a much lower tendency to grow (blue colours) compared to others (red coloured outer layers near the QC) and inhibit and distort growth ( Figure 7B). Alternative scenarios based on spatial or temporal derivatives of the auxin concentration pattern (as proposed for instance for Drosophila organ development [47]) do not provide more regularly spaced developmental cues (results not shown). We conclude that in itself a strictly auxin-based growth model is unlikely to produce stable and realistic root growth. Whereas a simple diffusion-based auxin transport mechanism (e.g. with a source near the 'QC') could resolve this issue, this spatial complexity resulting from polar auxin transport is important for other purposes, like auxin's role in regulating formative divisions around the 'QC' [17]. In the following we investigate two opposing theories proposed to circumvent this differential growth problem in response to a lateral (radial) auxin gradient.
It has been recently highlighted that hormones seem to regulate root growth in a layer (/tissue)-specific way [9]. Whereas short-range (paracrine) signals have been put forward as a hypothetical solution to transmit developmental signals transversally through the layers [9] we have tested whether it is sufficient for cell layers to 'communicate' through their mechanical connections. Model 11 (Tables 1 and S1) represents this concept of 'layer (/tissue)-driven growth' by selecting one arbitrary cell file (represented by 2 cell columns) in which the cells undergo turgor-driven growth. The other cell files are forced to passively follow (Figure 8, see legend and Dataset S1 for implementation). Simulations show that one cell file can effectively drive steady growth of the whole organ, yielding smooth cell length distributions (Figures 8B and 8C). Strictly, a conflict exists with the ULSR since one cell file is consistently growing faster. Nevertheless, if the cells in the other layers can respond rapidly (adjusting their target areas or 'turgor pressures') then the difference in strain rate apparently remains sufficiently low to cause much overall tissue distortion. This rapid response requires transversal communication between the layers. Biologically, signal transduction could potentially be mediated by the mechanical connections between the layers or otherwise a secondary chemical signal would have to be transmitted. Irrespectively, we conclude a layer-driven mechanism can solve the problem of layer-specific differences in the instructive hormone (here auxin) gradient.
A different solution to reconcile the central role of auxin in root growth with the ULSR may consist of introducing regulatory cross-talk into the models. In other words to complement or modulate the auxin signalling other spatially non-uniform signals, like transcription factors [17,48] or other hormones [49][50][51][52] are implicated. Based on proposed regulatory interactions [53,54], we have implemented a model (Model 12, Tables 1 and S1, Text S1) centred around the cross-talk between auxin as the cell proliferation promoting factor and cytokinin as the cell differentiation promoting factor ( Figure 9A). Auxin transport was described as both diffusive and PIN-mediated, with local production and an external source as in the previously discussed models. Cytokinin transport was described as strictly diffusion-based, also with local and external sources, the local production rate being repressed by auxin. This led to a shallow gradient, the concentration decreasing towards the tip (as in [55]). In this simplified model cytokinin upregulates the transcription factor SHY2 which is typically strongly expressed at the so called transition zone between DZ and EZ [54]. Increasing SHY2 leads to repression of PIN transporter levels and a resulting auxin signal drop that reduces the meristem size. Accordingly a local auxin response minimum has been observed before [18]. Auxin on the other hand has been proposed to counteract this effect by mediating SHY2 degradation, in part through gibberellin (GA) repressing ARR1 via the DELLA family protein RGA. In our simulation model we did not include any ARRs as variables. In accordance with others (see for instance [56]) auxin was assumed to directly inhibit cellular cytokinin production. GA was only included as a variable that determines where the EZ ends due to simple growth dilution (in accordance with [19]). The exit from proliferation and start of accelerated growth was determined by crossing a SHY2 concentration threshold (cf. Text S1). Interestingly, because of its inhibition of PIN expression, the steep increase in the SHY2 concentration coincides with a dip in the auxin concentration as well as the transition from a strongly polarized to a predominantly a-polar transport regime (Figures 9B and S10A) leading in general to much more linear spatial gradients in accordance with the ULSR (Figures 9C-E and S11). Since the first two proposed criteria for realistic growth were fulfilled too ( Figures S10B,C), this mechanism in principle allows for the required regular and stable zonation of the growing root apex, and can therefore potentially  (Table S1 -Model 9) with cells growing exponentially (in accordance with [32]) at the same low specific growth rate (except the 'cap' cells) up to a fixed distance of 240 mm from the tip and then undergoing 10-fold growth acceleration (up to a specific distance of 750 mm). Cell division was based on a sizer mechanism with noise added to have a more irregular, natural, tissue shape (up to a fixed distance of 240 mm from the tip). To not bias the analysis, growth and division did not depend on auxin. Four simulations based on differences in kinetic parameter values as well as some subtle differences in the differential equations were set up to represent different types of auxin sources in the model (see also Text S1, Figure S5 Table S1 -Model 9), with both external and local (area-dependent) auxin sources (and sinks). Kinetic parameters: s = 100 (mm 2 min) 21 ; F = 2.10 7 min 21 . Figures S6, S7, S8, S9 illustrate the dependence of the shape of the auxin gradient on the parameters used here (in a non-growing root). More information on the kinetic equations can be found in Text S1. doi:10.1371/journal.pcbi.1003910.g006  Table S1) in which cell division and (slow) growth are only possible above a fixed threshold of auxin concentration. Below this threshold (13.5 AU) and above a second lower threshold (8.8 AU) cells undergo accelerated growth. From an early stage of growth on strain rates are unbalanced leading to tissue distortion. The malformations accumulate and cell divisions are predominantly taking place in the central layers as determined by the auxin gradient. Colouring according to areal strain rates ('AS', cf. Methods) (B) Snapshot of a simulation with Model 10, but with a more dominant diffusion regime (parameters as in Figure 6D, with as threshold for accelerated growth ,60000 AU and for growth termination ,40000 AU) leading to a less pronounced lateral gradient (at 90 h). This produces less severe tissue distortion, but still severely inhibits growth and leads to unrealistic cell size distributions. Colouring is according to growth potential ('GP', as defined in the Methods section) as a measure for 'turgor pressure', showing a central region at the apex which opposes growth of the outer cell layers (indicated by blue versus red colours). doi:10.1371/journal.pcbi.1003910.g007 reconcile multiple roles for auxin in patterning along different directions. Furthermore, with Model 12 we could qualitatively predict how adding external auxin and cytokinin affects the meristem size (increasing and decreasing in size, respectively) as demonstrated experimentally by Beemster and Baskin ([34]; Figure 9C-E). We conclude that this computational model most effectively captures the basic growth characteristics of the Arabidopsis root and represents an ideal starting point to develop more advanced computational kinematic models which can predict root growth under more diverse conditions and perturbations.

Discussion
We have constructed and simulated different models that represent steady symplastic growth of the Arabidopsis root tip. We compared diverse regulatory mechanisms and found out which of them can adequately reproduce crucial properties of primary root growth, according to three well-defined criteria (steady-state, realistic cell length distribution, and ULSR) that allow rigorous comparison with experimental observations of in vivo growing roots of Arabidopsis thaliana. A necessary condition for stable unidirectional symplastic growth [33,26] was re-interpreted and re-formulated as a strain (rate) rule ('ULSR') for those mechanisms to conform to. It is mainly this third criterion that warrants the use of our sophisticated vertex-based simulations rather than the more simple approach of modelling a single cell file such as in [19], which could potentially produce kinematic output such as steady state length growth and cell length profiles. As soon as differences in strain rates or complex transport phenomena occur, the use of a two-dimensional tissue representation becomes necessary. Our simulations do suggest that symplastic structures are resilient to differences in strain rates as long as they are small and short-lived. How far or fast perturbations are precisely transmitted through  (Table  S1). Simulation time 109.5 h of model for which auxin concentration is 'interpreted' by the two layers of border cells (in analogy with endodermisspecific growth regulation by GA [78], a different tissue layer, for instance the epidermis, could be equally effective: result not shown) and translated into an increase in the target area of those cells (cf. Methods). The other cell layers are programmed to follow passively by re-setting their target areas to their actual areas after every simulation step in accordance with a small resisting force w.r.t. the layer that is controlling growth. Colouring is according to growth potential ('GP', as defined in the Methods section) as a measure for 'turgor pressure', showing border cells drive growth of neighbouring cells to the extent that their target areas are smaller than their actual areas (slight blue colour). (B) Plot of root length versus simulation time shows steady linear organ growth from 94 h on after a long preparatory phase to construct a realistic starting grid with a stable auxin gradient (code details in Dataset 1). (C) Plot depicting the cell length along the principal growth axis at step 103.5 h of the simulation with a model equivalent to Model 8 but with the growth driven by the 3th and 10th layer as in Model 11. Note that cell lengths vary smoothly from DZ to EZ similar to Figure 2 and 5C. doi:10.1371/journal.pcbi.1003910.g008 plant tissue and to which extend this affects organ growth are intriguing questions that go beyond the scope of this study. A more accurate representation of cell wall mechanical properties (such as [21]) and of tissue and organ structure (possibly in 3D) would then be desirable and would further enhance kinematic and microscopic resemblance with real root tissues.
Considering the complexity of molecular interactions implicated in growth regulation of the Arabidopsis root apex and considering the presence of many gaps in our understanding, we have opted rather to define rules for cell growth and division inspired by previous studies and enriched them with molecular interactions whenever useful. Whether regulation takes place at the supracellular ('organismal') level versus at the level of individual cells ('cellular') forms the subject of a long standing discussion [27]. This distinction is probably overly polarized as indicated by experimental observations [29][30][31]. In fact, signals that vary on the cellular, tissue-and organ-scale are known to affect pattern formation considerably. We have kept this conceptual distinction nevertheless and classified the proposed regulatory mechanisms as cell-autonomous (based on local pre-programmed rules) or noncell-autonomous (affected by external, spatial signals).
Cell-autonomous mechanisms like timers, counters and sizers are readily implemented with the logical expressions of a computer program. In a biological context such mechanisms typically require more complex designs [57][58][59]. Various sizer and timer based cell cycle models have been reported, making it reasonable representations of cell behaviour [38][39][40]. Given the conserved nature of the eukaryotic cell cycle machinery these types of mechanisms are likely operating in the plant cell cycle [41,60]. Timers and counters operating at different spatial and/or temporal scales are even more obscure. For instance, a developmental counter that determines the exit of proliferation would have to keep track of the number of cell divisions a cell has undergone. The epigenetic state of a plant cell can reflect this history [61]. In the context of plant development, telomere shortening could play such a role [62]. After cells exit the DZ, DNA duplication cycles are expected to continue through the process of endoreduplication. The number of DNA copies could therefore serve as a direct developmental marker to the plant cell.
We have shown that stable growth is possible based strictly on counter and timer mechanisms. However, the absence of spatial cues precludes realistic primary root growth. Indeed, such strictly cell-autonomous mechanisms will logically lead to groups of (nearly) synchronously growing cells. Cell packages derived from each division of the initial cells will behave in a similar way irrespective of their position along the growth axis and produce growth zones changing periodically in size (when such a package leaves the DZ and enters the EZ). Even taking into account the inherent noise in cell behaviour, fixed developmental zones and smooth transitions in cell lengths are not feasible based on this type of regulation. On the other hand it cannot be excluded that some indirect mechanisms exist by which cells can circumvent the need for a spatial signal by deriving spatial information in an autonomous way. An example of this hypothesis could be the  (Table S1). (A) Schematic view of regulatory interactions between model variables (in italics) and PIN exporters. Dotted lines illustrate potential cross-talk with gibberellin (GA) signalling (auxin stimulating GA and GA inhibiting cytokinin signalling) not included in the model. GA is represented in the model as an independent signal that undergoes growth-dilution, thereby determining the exit from elongation [19]. gravity-sensing columella cells that can extract spatial information through statholits in the root gravitropic response as modelled in [15].
Considering the inherent limitations of strictly cell-autonomous mechanisms and given the extensive evidence on the importance of phytohormones, it seems unlikely that those mechanisms are the main determinants of developmental zonation in the Arabidopsis root. We have demonstrated that spatial signals (for instance stemming from a biochemical gradient) can be a direct and (based on the visual and kinematic comparison) effective manner of instructing morphogenesis. This has been shown for various other life forms like Arthropoda and Vertebrata [63,47,37]. Our simulations regarding the effectiveness of auxin as the primary signal controlling root growth show that based on local auxin production a stable auxin pattern can be produced, however this potentially a slow process. A constant production per cell seems the most effective to this end. A faster breakdown [12] or efflux rate might aid in faster convergence. In any case, with an external auxin source such a pattern gradually fades out through growth dilution. Polar transport results in a lateral (radial) concentration gradient which conflicts with the 'ULSR'. In fact our model does not even capture the extra volumetric dilution of auxin from the inner to the outer layers in three dimensions. Similar patterns were obtained by other studies like Grieneisen et al. [12] and Santuari et al. ( [18]; their model did also not include the apoplastic compartment) and supported by various reporter studies (e.g. [64,45]). Although Grieneisen et al. [12] have simulated stable (2D) growth of the root tip, the lateral gradients were not recognized as problematic probably due to the fact that relative movement of cells (cell sliding) is possible with their (cellular Potts instead of vertex-based) framework. We have proposed various strategies to circumvent the ULSR conflict. It has been pointed out that several hormones are exerting their regulatory effect on the root in a cell-layer specific way [9,14]. This provides a way out of the 'ULSR' conundrum if accompanied by rapid transversal transmission to the other tissue layers. Candidate molecules to act as secondary transported signal are only just surfacing. On the other hand, even layer-driven growth by direct mechanical transduction was successful in producing a realistic root phenotype according to the three defined criteria. The role of auxin must by all means be understood in the complex context of various downstream response factors (with variations in levels, localisation, etc. (e.g. [48,18]) and also of other hormones that interfere via their respective signalling pathway components. We constructed a model based on the antagonistic role of auxin and cytokinin in root development, with the SHY2 transcription factor as a central regulator of meristem size [65][66][67]53,54] and gibberellin (GA) dilution determining cell maturation. Simulations with this model were in accordance with the ULSR and reproduced visual and kinematic observations as well as the expected increase and decrease of meristem size upon addition of auxin and cytokinin [34]. By down-regulating PIN-mediated transport through the transcription factor SHY2, cytokinin effectively flattens the lateral auxin gradient, at the basal end of the division zone ('transition zone' [6]), thereby signalling the exit of proliferation and start of differentiation without conflicting the ULSR. GA dilution has been proposed before as part of a more intricate mechanism (including cell compartments and DELLA proteins, [19]) determining the exit of the elongation phase in a single cell row of the Arabidopsis root. We have used growth dilution of GA in a simplified form with a GA production term proportional to (initial) cell width to account for layer-specific differences in GA concentration. Cell elongation is defined to stop as soon as the GA concentration drops below a specific GA minimum, which depends on reaching a certain cell size (and accordingly length). This could account for the similarity in cell lengths at the end of the EZ which has been observed under different conditions in the Arabidopsis root [34].
Various definitions exist for the term 'robustness' [68]. We define it here in a broad sense as a property that allows a system to maintain its functions in the context of internal and external perturbations [69]. 'Function' pertains to maintaining stable growth upon internal changes or changing growth conditions (external). Some of our simulations showed that, in the absence of feedback from spatial cues, and driven by cell proliferation, cellular defects (for instance leading to differences in cell cycle time) can propagate. In contrast to cell-autonomous regulation a more dictatorial or hierarchical form of regulation (in this case via organ-scale spatial signals), is more robust to local perturbations. A possible trade-off of such a dictatorial system might be its sensitivity to mutations in the controlling network itself [70] as supported by experimental studies on PIN mutants [71,72]. At the same time, such a high sensitivity to modulations in the control system may be beneficial for plants, being sedentary and therefore requiring a great deal of developmental plasticity. This trade-off also applies to the degree of dependence of cell growth and division. In theory, a completely independent regulation of these processes allows for a much wider morphological space to be accessible. Evidently, reaching stable balanced growth would become far more challenging. As the mechanistic coupling of growth and division is still poorly characterized [73], we opted here to couple those processes to various degrees in an artificial way, for instance as synchronized counters and timers, or concurring spatial transitions. Integrating an explicit model of the cell cycle in our models would undoubtedly further progress our understanding of how the coordination of those processes affects root morphogenesis.
Starting from a simple grid and by stepwise evaluation of regulatory mechanisms for cell division and expansion, we have obtained a model that produces realistic root growth according to well-defined criteria. In that respect spatial regulation appears to be essential, presumably with auxin as the central signal that is modulated by feedback interaction with cytokinin. This represents an important step towards predicting the Arabidopsis root phenotype under various conditions using vertex-based modelling. As indicated above, more advanced representations of cell wall mechanics, cell cycle regulation, etc. would benefit this model, and are currently being developed.

Methods
A vertex-based modelling approach based on the modelling software VirtualLeaf [25] was used to represent a growing root apex that consists of nodes (vertices) connected by edges. Cells are polygons and cell wall segments correspond to the edges. Mechanically the edges are acting as linear springs which are opposing the tendency of cells to increase their area by virtue of their target areas. This is mathematically expressed as terms in a generalized energy function (or 'Hamiltonian'), which was slightly modified to represent balanced linear growth (in the sense that small cells contribute equally to the energy function as larger cells, if they have identical relative growth rates).
Original Hamiltonian: Where indices i and j sum over all cells and polygon edges, respectively, l A is a parameter setting the cells' resistance to compression or expansion, and l M is a spring constant. A T is the cell's target area, L T the wall element target length.
In the new Hamiltonian the first term is replaced by: The cell mechanics algorithm in VirtualLeaf randomly displaces all nodes to steer the system towards a global energy minimum. During this Metropolis Monte Carlo simulation any node displacement that leads to an energy decrease (or a small energy increase corresponding to a high Boltzmann probability) is accepted [25]. Added example input xml files contain the mechanical and numerical integration settings for simulation and illustrate the data structure used to describe the state of the cellular grid (Dataset S1). Besides the mechanical framework cells and cell walls are endowed with biochemical properties represented by reaction and transport equations and logical rules that determine when and how cell division occurs and how the cellular surface energy terms change. Our approach is to alter the structure of those equations and rules rather than just change parameter values to compare the output of the model versions.
The cellular organisation of the root tip was inspired by an earlier study [12] and always starts from a regular grid of rectangular cells, sometimes with narrower inner cell layers roughly in accordance to root anatomical characteristics [7]. The upper cell row always consisted of fixed nodes to enforce downward growth. Across all simulations the three lowest cell rows are not driven to grow in size (via their fixed target areas) and are meant to represent the complex of columella cells together with the quiescent centre (QC). Although this part of the root apex (root cap) has its own detailed developmental characteristics, we strictly focus here on the proximal/basal side of the QC while excluding the distal meristem as part of the growth process. The formative asymmetric divisions of the stem cells surrounding the QC [74] were neither included.
The models in principle describe growth in two dimensions but it was assumed that transversal expansion is severely constrained, which is in correspondence with experimental data indicating little or no tangential expansion (,1% h 21 : [75]).
Cell walls behave as elastic springs with vertically oriented walls contributing less to the energy function (typically scaled down 5 or 10 times). This is sufficient to allow fluent elongation of cells along the root axis while keeping cell shape adequately constrained as required by turgor pressure acting on the cell walls. The output repertoire is limited both by our a priori model assumptions and by computational limits compared to what can be observed in vitro or in vivo. Briefly, lateral root formation and tropic responses such as gravitropism, coiling, waving, skewing are not produced. Model simulations do allow assessment of root growth rate and size, as well as detailed kinematics (the latter is in principle even more flexible than in a classical experimental setting where sampling has temporal and spatial limitations. It is also possible to assess deformation of overall organ shape to some degree (given the strong horizontal constraints) but especially defects in internal tissue organisation. At the cellular level the resolution of the simulations is sufficient to observe cell shape deformations.
The experimental data for comparison of simulated cell length distributions with the real Arabidopsis primary root ( Figure 2) were obtained by kinematic analysis of one typical seedling grown in the same conditions as described in Beemster and Baskin [5]. Curve fitting was in general obtained by smoothing and interpolation at an interval of 5 to 25 mm using the kernel smoothing function 'locpoly' of the KernSmooth package [76] for the R statistical package (R foundation for Statistical Computing). Table S1 presents a detailed overview of the models described in this paper. Except for Model 9 and Model 12 (cf. Figure 6 and 9), all models started off with 2 divisions of the cells in the row above the 'QC' (fourth row from below Figure 1B) to reduce the initial size of the cells and in some cases to introduce some asymmetry (noise) in that cell row. For cells to enter a new developmental phase requires the definition of a relative starting point in time. We have modelled this by associating a clock to each cell which is reset to zero as long as the daughter cells after cell division remain in contact with the nongrowing cap structure. This results in cells in contact to the QC (hypothetically sensed by paracrine signalling for instance) remaining part of the stem cell niche ( = capable of self-renewal, 'stemness' [77]) and dividing at a pre-set frequency. Upon loss of contact with the QC-cap structure, the cells assume the identity of a typical meristematic (DZ) cell until eventually losing the potential for division after for instance a fixed time or number of divisions.
Cell-autonomous timer, counter, and sizer mechanisms were implemented similarly to the exemplary fragments of VirtualLeaf (pseudo-)code in Text S1. Pseudo-code fragments are also presented for the final Model 12 in Text S1.

Visualisation
In several instances visualization of simulation output consisted of colouring the cells according to either: 'AS' -Areal strain (rate): To obtain a simple, approximate measure of the local expansion of the cells, the change of the area relative to the area in the previous simulation step was accounted for.
Red colouring typically indicates a positive strain rate and blue colouring a negative strain rate.
Or: 'GP' -Growth potential: The ratio CellArea actual CellArea t arg et was used to colour the cellular grids in specific instances as a measure of the cellular growth potential (broadly interpreted as the 'turgor pressure'). Red colouring indicates an energy driving areal growth, blue colouring an energy growth driving cell shrinking. Or: 'AU' -Arbitrary concentration Units, with red to yellow colouring indicating increasing concentration values.

Supporting Information
Dataset S1 Model code.  [1][2][3][4][5][6][7][8]10,11), Wortel (2).xml (for model 9), and Wortel (3).xml (for model 12). They specify the typical parameter values for the numerical integration and Monte Carlo methods, input/output behaviour, as well as some kinetic constants in the header of the file. This is followed by the list of nodes, cells and walls of the cellular grid including their geometrical and biochemical attributes. (ZIP) Figure S1 Similar output for counter-and timer-based models. (A) Simulation output of Model 3 (Table S1) with the exit from proliferation defined by a counter mechanism. The imposed growth and division rules have resulted in a highly regular grid with distinct zones of similar cell length. (B) Simulation output of the timer-based Model 2 (but here without any noise added to the starting divisions of the tissue). This yields a very similar grid as in (A) at 99 h simulation time (the small differences are due to a few nodes in close proximity that have not collapsed due to the stochastic character of the Monte Carlo mechanical framework).  Table S1, see also Figure 4A). Note the smoothened curve. The '*' indicates from where steady growth starts. (B) Output of Model 7 (Table S1). Upon release from the QC cells undergo 3 divisions based on reaching a cell layer-specific size ('sizer'). As for other strictly cell-autonomous mechanisms, cells belong to groups of similarly sized and synchronously growing cells. Cell division is less synchronized which leads to a smoothened increase in cell numbers. (C) Cell length along the growth axis at time step 91.5 h shows broader cell length distributions (blue dots) when noise is added (Model 4, Table S1) compared to the red dots produced with Model 2 (Table  S1, same data as in Figure 3C (D) Velocity profile as in (C), but obtained through binning and averaging data from ten randomized simulations. The shape of the curve is smoother than in (C). The data points were fitted with a progress curve using the 'polyloc' method (cf. Methods). (TIF) Figure S5 Temporal auxin patterns in a growing root. Evolution of the total auxin levels (A,D,G), the total simulated root area (B,E,H), and the auxin concentrations (C,F,I), corresponding to Figure 6A-C. A strictly external source of auxin combined with local sinks results in saturation of the total auxin in the root (A), which together with a steady increase in surface area (B) leads to the auxin concentration converging to zero (C). A strictly internal auxin source (production rates proportional to cell areas) results in the total auxin level (D) increasing proportionally to the area increase (E), and the auxin concentration slowly converging to a steady state (F). A strictly internal auxin source (production rates constant per cell) results in similar patterns as for (D-F) with a steady state in auxin concentration (I) as a function of the total auxin level (G) and area (H).  Figure S6) should lead to similar output. Fixed parameter values were as in Table S2. (TIF) Figure S12 Model 12 is robust to small parameter variations. The simulation output of Model 12 (yellow colouring; arbitrary units: AU) is shown here for a 10% increase of different parameter values related to hormone transport: (A) simulation based on the reference parameter set (Table S2); (B) D[0] perturbed; (C) k ex perturbed; (D) k ex perturbed; (E) s perturbed; (F) D [1] perturbed. The output is highly similar, which is also the case if these parameter values are decreased by 10% (results not shown), demonstrating local robustness/stability of the simulated output to changes of these parameters. (TIF) Table S1 Model overview. Overview of the models used in this study. Various categories w.r.t. developmental decisions are presented. Column (3) specifies the transition between division and elongation zone (DZ and EZ, respectively) with in parentheses the number of division or time since release form the QC; column (4) specifies the transition to mature (differentiated) cells based on timing since the release from the QC or a spatial signal at a fixed distance from the root apex; column (5) specifies whether division rate is determined via a timer or sizer mechanism; and column (6) how cellular growth rates are defined. Developmental events can be determined to happen after a fixed duration ('Timer'), a fixed number of divisions ('Counter'), a fixed cell size ('Sizer'), and a fixed distance from the root apex ('Positional', i.e. a ruler). For Models 10-12 more complicated regulatory mechanisms are specified. In Model 4 extra random noise was added to the timer (+/2 max. 25%). For Model 5 cell division is dependent on a timer mechanism, with a different cell cycle time for inner and outer layers. Model 6 uses a uniform size (area) criterion for cell division, irrespectively of geometry. For Model 7 the size criterion is adapted to the difference in the width between inner and outer layers, thereby acting effectively as a length criterion. In the other (non-cell autonomous) cases the cell division sizers differs between layers. It is usually assumed that DZ and EZ have a characteristic elongation rate (expressed as a constant addition to the target area or a fraction of the current cell area). Layer-specific differences are indicated (by column numbers) if necessary. All models are identical in that the columella/QC cells (lower 3 cell rows) cannot divide. For timer and counter based models the contact with the QC cells is important to keep them dividing. As soon as they lose contact the respective counter or timer starts to increase. For noncell-autonomous models a morphogen threshold determines if cells keep dividing. For each model cell division occurs through the insertion of a new horizontal wall dividing the cell in half. For models using hormone transport parameter values were used in accordance with known modelling studies such as [12,17], and Wabnik et al., 2010 [79]. In general reaction and transport equations are discussed in Text S1. A detailed parameter description can be found in the respective figures or in Table S2. (DOCX) For Model 12 the auxin reaction and transport parameters D[0], k_export, and k_import, were adjusted to produce an auxin gradient with a more pronounced and central maximum such as found in various reporter studies (see for instance Fig. 3c in [12] and Fig. 1A in [18]). Figures S11 and S12 illustrate parameter dependent behaviour and sensitivity. The model behaviour is robust to hormone transport parameters changes and even larger changes can in principle be accommodated based on the balance between a-polar and polar auxin transport and stable negative feedback regulation of auxin and cytokinin.

(DOCX)
Text S1 Background information on model implementation.

(DOCX)
Video S1 Realistic root growth through non-cell-autonomous regulation. Movie obtained through simulation of Model 8 (Table S1) for 84 hours shows that steady growth is reached after a transient time (of roughly 30 hours). Colouring is in accordance with areal strain rates ('AS', cf. Methods) and shows that the EZ, which is typically more intensely coloured, reaches a stable size, while moving down. (MP4)