Parsimonious Model of Vascular Patterning Links Transverse Hormone Fluxes to Lateral Root Initiation: Auxin Leads the Way, while Cytokinin Levels Out

An auxin maximum is positioned along the xylem axis of the Arabidopsis root tip. The pattern depends on mutual feedback between auxin and cytokinins mediated by the PIN class of auxin efflux transporters and AHP6, an inhibitor of cytokinin signalling. This interaction has been proposed to regulate the size and the position of the hormones’ respective signalling domains and specify distinct boundaries between them. To understand the dynamics of this regulatory network, we implemented a parsimonious computational model of auxin transport that considers hormonal regulation of the auxin transporters within a spatial context, explicitly taking into account cell shape and polarity and the presence of cell walls. Our analysis reveals that an informative spatial pattern in cytokinin levels generated by diffusion is a theoretically unlikely scenario. Furthermore, our model shows that such a pattern is not required for correct and robust auxin patterning. Instead, auxin-dependent modifications of cytokinin response, rather than variations in cytokinin levels, allow for the necessary feedbacks, which can amplify and stabilise the auxin maximum. Our simulations demonstrate the importance of hormonal regulation of auxin efflux for pattern robustness. While involvement of the PIN proteins in vascular patterning is well established, we predict and experimentally verify a role of AUX1 and LAX1/2 auxin influx transporters in this process. Furthermore, we show that polar localisation of PIN1 generates an auxin flux circuit that not only stabilises the accumulation of auxin within the xylem axis, but also provides a mechanism for auxin to accumulate specifically in the xylem-pole pericycle cells, an important early step in lateral root initiation. The model also revealed that pericycle cells on opposite xylem poles compete for auxin accumulation, consistent with the observation that lateral roots are not initiated opposite to each other.


Introduction
The development of vascular tissue was a critical innovation in the evolutionary history of plants, providing a mechanism for long-distance transport and structural support that was crucial for their colonisation of land. In order to fulfil these functions, vascular tissues must be appropriately positioned to ensure the continuity of the vascular strands, forming a transport and communication network between distant tissues and organs. The vascular tissues also comprise wood, which is a critical global resource. Improving our understanding of vascular development in plants is therefore an important challenge with yields for both fundamental and applied biology.
The precise and consistent layout of the vascular cylinder (stele) of Arabidopsis thaliana makes it an ideal system in which to study vascular patterning. The Arabidopsis stele is comprised of a central xylem axis flanked by intervening procambial cells; these separate the xylem from the phloem poles, which are located perpendicular to the xylem axis ( Fig 1A). The xylem axis consists of two cell types, protoxylem and metaxylem; protoxylem cells are positioned at the margins of the axis and mature earlier than the metaxylem cells in the centre [1]. Specification of the protoxylem cells is one of the earliest events in post-embryonic vascular patterning [2] and so has been the focus of attempts to unravel the genetic and hormonal interactions involved in vascular patterning. It has recently been suggested that the correct positioning and specification of protoxylem in the Arabidopsis stele depends on a mutually inhibitory interaction between two major classes of plant hormones, auxin and cytokinins [3]. Although recent modelling efforts are consistent with this proposition [4], the relevance and plausibility of an informative spatial pattern in cytokinin levels in this context remains unclear.
Auxin and cytokinins have long been known to interact antagonistically in a wide range of processes in plants [6][7][8], such as meristem maintenance in the shoot [9] and root [10], control of axillary branching [11,12], and lateral root initiation [13]. Recently, researchers have begun to elucidate the molecular and mechanistic bases of this antagonism. Auxin and cytokinins are known to influence each other's synthesis [14][15][16], metabolism [16,17] and signalling and response machinery [3,9,10,[18][19][20]. Moreover, cytokinin has emerged as a key regulator of auxin transport [3,10,21], and this regulation has been shown to play a central role in protoxylem specification [3]. Auxin transport depends on the activity of molecular transporters; while auxin can either passively permeate into cells or enter them via AUX1/LAX importers, its efflux must be mediated by exporters such as the PIN class of efflux carriers [22,23]. Auxin may also move symplastically via plasmodesmata [24], and transporters of other families, such as the PGPs [25][26][27][28] and NRTs [29,30], may interact with the PINs or otherwise modulate auxin movement. The PIN proteins, however, are the only auxin exporters known to be polarly localised within cells, and their localisation along specific regions of the plasma membrane has been shown to be an important factor determining the overall auxin distribution and flux pattern [31][32][33]. In the current work, we investigate how cytokinin signalling The network is active in all cells but shown in only two, representative of the procambium/pericycle and the xylem axis; the faint elements are proposed to be downregulated in the corresponding domain. (I) Simulation of a longitudinal section of a root shows uneven auxin accumulation throughout the stele, with outer cells (pericycle) accumulating higher levels of auxin (parameter settings and PIN localisation as described in [5]). Arrowheads mark the protoxylem position in cross sections; scale bars represent 10μm. In the model schematic, yellow arrows represent PIN transporters and pink arrows represent passive auxin influx. modifies the auxin distribution and flux via the PIN transporters in the Arabidopsis root tip and explore the possible functional consequences.
Cytokinin signalling is essential for correct vascular patterning. The stele of wol plants, which have a mutation in the cytokinin receptor CRE1 leading to a severely attenuated cytokinin response, is smaller than in wild type and is composed entirely of protoxylem cells [2,34]. Cytokinin inhibits the specification of protoxylem; repression of cytokinin signalling at the protoxylem position by AHP6 (Fig 1B) is required for normal protoxylem specification, and treatment with exogenous cytokinin or a mutation in AHP6 results in the spread of cytokinin signalling to the protoxylem position and a concomitant loss of protoxylem [35]. A domain of high cytokinin signalling is observed in the procambial cells in cross sections of wild-type roots (Fig 1C, [35]) with a complementary domain of high auxin signalling in the xylem axis ( Fig  1D, [3]). Expression of PIN1 and PIN7 overlaps with the domain of cytokinin signalling ( Fig  1E and 1G), which has been shown to regulate the expression of PIN7 and the localisation of PIN1 on lateral plasma membranes in the procambium [3]. Coupled with the specific expression pattern of PIN3 (Fig 1F), this is proposed to result in transport of auxin away from the regions with high cytokinin signalling and towards the xylem axis, where it accumulates and promotes xylem identity. Auxin accumulation also activates AHP6 expression at the protoxylem position [3]. Given that AHP6 in turn inhibits cytokinin signalling [35], the data suggest that a feedback loop involving hormone transport dynamics and effectively mutually inhibitory interactions provides a mechanism for positioning the xylem axis in Arabidopsis (Fig 1H, [3]). While AHP6 expression along the xylem axis is regulated by the HD-ZIP III transcription factors [36], the mechanism by which cytokinin signalling is excluded from the metaxylem remains unclear, though the weak phenotype of the ahp6 mutant suggests that a redundant partner of AHP6 may be active in the metaxylem.
It is clear that cytokinin signalling plays an important role in vascular development. The patterns in cytokinin signalling, such as those observed via pARR5::GUS, are a combined consequence of both the spatial distribution in the levels of cytokinin and the strength with which individual cells respond to those levels of cytokinin. Since it is known that hormonal cross-talk can regulate the cellular cytokinin signalling and response machinery, it remains unclear to what extent patterns in cytokinin signalling reflect actual underlying distributions in cytokinin levels rather than reflecting patterns in the response machinery and hence indirectly the auxin pattern in the tissue, since this response depends on auxin. In the most extreme scenario, cytokinin could even be spatially homogeneous, with variations in signalling resulting solely from the expression patterns of the receptors and auxin-mediated modifications of the downstream effects of cytokinin. To shed light on the auxin dynamics guiding vascular patterning and the relative importance of cytokinin levels versus the cytokinin response machinery in this process, we here investigate auxin-cytokinin interactions in a realistic spatial context. Computational modelling provides tools with which to address complex feedbacks such as these and has been applied successfully in similar contexts [31,37,38]. Our model identified additional components of the genetic-hormonal network coordinating vascular development and has uncovered a mechanism by which regulated transverse auxin transport increases the robustness of the xylem-axis maximum and, unexpectedly, impinges on lateral root initiation.

Computational Methods
Simulations were conducted using a grid-based spatial model of auxin and cytokinin movement in which cells occupy an experimentally reasonable area and the cell wall is explicitly included; inputs to the model were cell and tissue geometry, as well as cell polarity. While the model does not attempt to explain the polarity patterns themselves (it assumes them as an initial condition), it does allow us to explore the concentration and flux patterns resulting from the patterns in polarity over the tissue. Our focus is to study how transporter regulation can alter these patterns given the feedbacks through auxin and cytokinin. The model can be considered an extension of a previously published model of longitudinal auxin transport in the Arabidopsis root [31]. While the earlier model only considered auxin and the role of transporters in an unregulated fashion, the current work explicitly considers the regulation of transporter activity under the control of auxin and cytokinin. These hormones, in turn, dynamically move within the tissue due to diffusion and membrane transport properties. We strive for a parsimonious model that allows us to gain general insights into the patterning mechanisms while maintaining sufficient proximity to the core players involved to allow our model to be compared and contrasted to mutants and hormone treatments.
For the purposes of our model, we explicitly capture three exporters and one importer under independent control of the regulatory dynamics described in [3] and [39]. The PIN proteins are treated as identical in terms of factors such as affinity or transport efficiency, differing only in how they are regulated and localised. However, an additional weaker PIN is used to represent weak expression of PIN7 in the phloem position, differing only in having a lower permeability. We do not explicitly model PIN transcription and translation, nor do we include posttranslational processes that might impact PIN activity. Instead, we use a simplified approach wherein PIN 'expression' in our simulations represents the PIN activity level resulting from all of these processes. Furthermore, an additional, unregulated efflux transporter is included with apolar expression in the outer tissue layers (endodermis, cortex, epidermis), while the AUX1/ LAX importers are represented by an influx transporter with apolar localisation in all stele cells. This importer is auxin-upregulated only in the context of lateral root initiation. Finally, our model makes the simplifying assumption of uniform longitudinal flux through the cross section, which we implement by a uniform constant influx as well as a uniform efflux of auxin and cytokinin into and out of the plane of study. Note that within this parsimonious twodimensional (2D) modelling setting, influx into the cross-sectional plane cannot be distinguished from biosynthesis, nor efflux out of the cross-sectional plane from degradation. All the processes that are included in the model are described hereinafter.
Synthesis and degradation of auxin and cytokinin. Biosynthesis and degradation, or alternatively, influx into and efflux out of the 2D plane of study, of auxin and cytokinin occurs in all cells according to the following basic equation: where C represents the concentration, b the biosynthesis rate (or influx into the plane), and δ the degradation rate (or efflux out of the plane). When we simulated cytokinin gradients, b was confined to specific cell types in order to create a localized source of influx/biosynthesis, as an extreme attempt to generate an informative gradient, while δ was set to provide the desired gradient steepness; the value of b and K A cyt (see below) were altered in those simulations to retain consistent average cytokinin signalling levels. Parameter values for the different simulation conditions are provided in Tables 1-3. Transport and diffusion of auxin and cytokinin. Under default conditions, cytokinin was allowed to diffuse freely within cells, as well as across the plasma membrane and throughout the apoplast, with no change in the diffusion coefficient. As we have not found a published diffusion coefficient for cytokinin, we used the diffusion coefficient of auxin (from [31]), which is a molecule of approximately the same size, composition, and shape. Auxin was allowed to diffuse freely within cells and within the apoplast at a lower rate. Diffusion of both auxin and cytokinin is in accordance with Fick's law: Although auxin can diffuse freely within cells and in the cell walls, the plasma membrane represents a barrier to auxin diffusive flux. There is a small amount of "leakage" across the membrane, taken into account by a background influx and efflux permeability term. However, auxin flux across the plasma membrane occurs primarily via the action of transporters. In the absence of transporters, efflux is severely limited (represented by a low level of leakage), while influx occurs more readily but is still strongly increased by the presence of a transporter. Auxin flux across the plasma membrane is represented in our simulations by the following equations: J ¼ ÀðE pin P e PINn ÞC in þ ðE aux P i AUXn ÞC out þ ðP i bgn ÞC out for membranes with at least one PIN expressed; ÀðP e bgn ÞC in þ ðE aux P i AUXn ÞC out þ ðP i bgn ÞC out otherwise; : where n is the inward directed unit vector perpendicular to the plasma membrane, and C in and C out represent the auxin concentration immediately adjacent to the membrane at the cytosolic and cell wall side, respectively. E pin and E aux represent the relative activity level of the transporters, and P i and P e represent the permeability of the membrane to influx and efflux, respectively. The influx transporter is uniformly localised on all plasma membranes of all stele cells; influx permeability is separated into a smaller passive component, P i bg , and a larger active component due to the transporter, P i AUX . In the absence of efflux transporters, the passive efflux permeability P e bg is used, while membrane elements on which PIN transporters are localised have a higher efflux permeability, P e PIN . When multiple PIN transporters are localised at the same membrane element, their activities are summed.
Regulation of transporter expression levels. The expression level of the individual transporters was determined according to a hormonal regulation model based on experimental data. Auxin upregulates PIN3 [10] and, during lateral root initiation, AUX1 [40]. In contrast, when calculating the expression levels of PIN1 and PIN7, the concentrations of both auxin and cytokinin have to be taken into account [3]. PIN1 and PIN7 are upregulated by cytokinin, but simultaneously down-regulated by auxin, because auxin activates AHP6, thereby inhibiting the cytokinin response which would upregulate the PINs. The auxin-AHP6-cytokinin-PIN pathway has been simplified in this model as a direct downregulation by auxin of PIN1 and PIN7. These regulatory interactions are represented in our parsimonious model by the following equations: : where E up and E down represent the hormonal up-and down-regulation of the auxin transporters in our parsimonious model, captured by simple sigmoidal Hill functions of transporter activity (PINs and AUX1) as a function of the average hormone concentration within the cell (x): The parameter K A i represents the hormone level (x) that gives rise to a half maximum transport activity. It adopts three possible values, depending on the transporter/hormone combination: K A AUX for auxin-dependent AUX1 regulation; K A PIN for auxin-dependent PIN regulation; and K A cyt for cytokinin-dependent PIN regulation. Please see Table 2 for the values used. Please note that the values of K A cyt were chosen to give rise to an average PIN activity level of % 70% in the stele (% 20-50% in simulations with a steep cytokinin gradient). Cytokinin patterning. Under the default conditions, cytokinin movement in the model occurs solely via diffusion. A cytokinin gradient, such as shown in S7A Fig, was then generated by defining source cells (either the proto-and metaxylem or the phloem, as indicated in the figure captions) and altering the degradation rate to tune the slope of the gradient, as discussed above. In those simulations, a degradation rate of δ = 0.24s −1 was used, producing a very steep gradient with a characteristic length λ = 11μm and a root mean displacement of ffiffiffiffi ffi x 2 p ¼ 15mm. We also explored the potential impact of cytokinin confinement by the plasma membrane. While in most of the simulations cytokinin was considered to be able to diffuse freely throughout the cross section, in a subset of simulations the plasma membrane presented a barrier to cytokinin movement, implemented through a bidirectional membrane permeability for cytokinin, P cyt , which we varied. Cytokinin flux across the plasma membrane was then given by: Furthermore, these simulations included a separate and lower diffusion coefficient for cytokinin in the apoplast, D cyt apo , as was done for auxin in this study and in [39].
Auxin-driven cytokinin biosynthesis. To explore the possibility that auxin-driven cytokinin biosynthesis could be capable of generating an informative pattern in cytokinin levels, as recently suggested [41], we also modelled the feedback of the auxin levels on cytokinin biosynthesis. This was done by replacing the cytokinin synthesis parameter b (in Eq (1)) by a sigmoidal function of the average auxin concentration within the cell (parameters given in Table 4): wol mutant and cytokinin treatment. We implement the wol mutant by directly modifying both PIN localization and intensity to closely match our own experimental observations (see . In addition, the anatomy of the realistic wol section is based on microscopy of a wol plant in order to account for changes in size and cell number. Likewise, to implement cytokinin treatment, we again directly modified PIN localization and intensity following the experimental observations (S4D- S4F Fig). In both cases, cytokinin response and levels are not changed and the PIN modifications are not an indirect consequence of changes in signalling interactions.
Numerically solving the equations. Auxin and cytokinin diffusion, membrane transport, production and decay were calculated numerically on the discretised grid points. This was done by numerically solving in an integrative fashion both the reaction part and the spatial coupling of the above equations using an alternating-implicit direction method [42]. Diffusion takes place when two grid points belong to the same cell or when both are part of the cell wall. Permeability is involved whenever two grid points are separated by a plasma membrane. Production and breakdown take place at each grid point belonging to a cell. As usual, appropriate tests were performed to ensure numerical stability and absence of grid effects. Heat maps. Concentrations of auxin and cytokinin are visually shown within the tissue context in two different manners: either through a linear rainbow heatmap ranging from minimum to maximum concentration values (except for Fig 1I, in which we use a logarithmic colour map), or as a "DR5" colour-scheme inspired by auxin response. The latter also uses a linear map scaled between minimum and maximum values, but is designed to use colour intensity to emphasize the regions with high concentrations. Transporter activity is shown using this same colour scheme, now representing the transporter activity level on a per cell basis. The colour bars are shown within the figures or referred to within the captions.
The laser intensity and gain settings were adjusted to not saturate the signal in control roots in order to see a reduction or change in pattern in the triple mutants and an increase in the TCS::GFP signal in the cytokinin movement experiment. Imaging was done on Leica SP5 II HCS A, GFP settings, with a 488 argon laser for both experiments.
The TCS::GFP images were analysed with imageJ. Intensity values were extracted from a rectangle 48 μm in height bounded by the epidermis on the upper left and the cell wall between the QC and vascular initials at the lower edge. These were then averaged in Excel and plotted as a scatter graph with error bars indicating the 95% confidence interval.
Quantitative RT-PCR. Col-0 seeds were plated on 0.5X MS growth media and grown in a Sanyo growth cabinet for 5 days. Approximately 50 seedlings were transferred to liquid 0.5X MS growth media for either mock (DMSO) or IAA (1 μM) treatments and incubated with shaking. After 0.5 h and 2 h time points 2 mm root tip segments were harvested into liquid N. 100 ng of total RNA isolated with Thermo Scientific GeneJET Plant RNA Purification Kit was DNase-treated and oligodT-primed cDNA synthesis was performed with Revert Aid Premium Reverse Transcriptase (Thermo Scientific). Diluted cDNA was used as a template in qPCR with HOT FIREPol EvaGreen qPCR Mix Plus (no ROX) (Solis BioDyne) and the primers listed in Table 5. The Bio-Rad CFX384 was used with one cycle 95°C for 15 minutes, 50 cycles each consisting of 95°C for 10 s, 60°C for 10 s and 72°C for 30 s, one cycle 95°C for 10 s followed by melt curve analysis. Data was normalised to ACT2 expression and relative expression calculated with the comparative CT method. Statistical significance was analysed using R packages nlme and gmodels by fitting a linear mixed model with a fixed effect for sample and a random effect for the biological repeat. The model contrasts were calculated and resulting p-values were adjusted for multiple comparisons. The experiment was performed in quadruplicate.

Results
In order to investigate the regulatory network discussed above, we implemented its dynamics in a spatially explicit computational model of auxin transport within a cross section of the root, taking into account cell shape, polarity of transporters and the presence of cell walls. In addition to cytokinin regulation of PIN1 and PIN7 [3], we also included auxin upregulation of PIN3 [10] and, during lateral root initiation, of auxin importers such as AUX1 [40]. While we present some of our results using a rainbow-coloured linear heatmap to display the auxin concentration, it is often useful to highlight the auxin pattern in a manner which can be more readily compared with experimental data. We therefore also use a blue heatmap ranging from the minimum to maximum auxin concentration with an intensity curve designed to be similar to the more familiar output of the auxin reporter DR5. (See Materials and Methods for further details.)

Lateral patterning of the PINs is necessary
The auxin distribution in the root can be predicted through simulations of auxin transport dynamics which take into account both transporter localisation and tissue layout [31]. Simulations based on realistic longitudinal sections of a root using a typical structured polarity pattern (for details, see [5]) predict both a longitudinal and a radial gradient of auxin within the vascular cylinder, with the highest levels found in the marginal cell file, which represents the pericycle ( Fig 1I). A naïve rotation of this pattern around the central axis of the root would yield a stele with higher auxin signalling in all pericycle cells. This is clearly inconsistent with the observed pattern of auxin signalling, which is radially asymmetric (Fig 1D, [3]). Models which analyse auxin dynamics within a longitudinal cross-section ignore transverse PIN localisation, which includes tangential PIN localisation (i.e., clockwise and counterclockwise PIN localisation within a transverse cross-section) and radial asymmetries in inward and outward PIN Transverse Hormone Fluxes in Arabidopsis Vascular Patterning localisation. However, the clear pattern observed in transverse cross-sections suggests that these elements may be important for the positioning of auxin along the transverse axis. To evaluate the role of the fluxes generated by transverse and radial PIN localisations, we implemented a computational model of auxin transport within transverse root sections. The flow of auxin and cytokinin into the 2D plane of observation (i.e. the root's cross section) was modelled by assuming a constant and homogeneously distributed influx of the hormones into the plane, while the hormones' efflux out of the plane was emulated through an effective homogeneous decay rate. Given that we focus here on studying the contribution to patterning by the radial and transversal fluxes mediated by transporters and their regulation, we did not assume any prepatterns in hormone entry or exit along the longitudinal axis, unless specified otherwise.
(See Material and Methods for equations and simulation details.) We conducted simulations with two different layouts, a 'realistic' cross-section that was generated through segmentation of an experimentally obtained confocal cross-section of the root, and a simplified, algorithmically generated 'geometric' cross-section (Fig 2A) to allow for a more general analysis of the dynamics. Our model takes into account the spatial structure of cells and the presence of cell walls as well as their organisation within a tissue. This is important because diffusion of the hormones occurs both within the cells and within the cell walls. In the case of auxin, diffusion is not allowed to occur over the plasma membranes; instead, known and regulated influx and efflux permeability rates over the plasma membrane are used. By contrast, cytokinin is allowed to diffuse unrestrictedly through the tissue, unless otherwise specified. An important feature of the model is that the PIN levels within a cell dictate the permeability rate of the auxin efflux through its plasma membrane, assuming a linear relationship between plasma-membrane localised functional PIN and PIN expression. Likewise, the presence of influx carriers, and their expression level, when regulated, determines the influx. Increasing transporter activity leads to a linear, non-saturating permeability increase in the appropriate direction (see Material and Methods).
Our spatial model explores the concentration and flux patterns resulting from underlying polarity patterns in the transporters, which we take as an initial condition to the model. For the purposes of our parsimonious model, PIN proteins are treated as identical in terms of factors such as affinity and transport efficiency, differing only in how they are regulated and localised. However, an additional weaker PIN was used to represent weak expression of PIN7 in the phloem position, differing only in having a lower permeability rate. Furthermore, an additional efflux transporter with apolar localisation is included in the endodermis, cortex, and epidermis, to represent the reported expression and localisation of PIN2 [51], PGP1 and PGP19 [26] in the outer tissue layers. Localisation of the AUX1/LAX importers is always apolar.
This template, and its underlying elements and formulations, lends itself readily to study how transporter regulation, due to the auxin-cytokinin feedback interactions explained above, can alter hormonal patterns by modifying the permeability contributions of the transporters. Specifically, the current work captures three exporters and one importer under independent control, taking into account the regulatory loop described in [3]. The model therefore includes hormonal regulation of the levels of the auxin transporters PIN1, PIN3 and PIN7, which are positioned according to experimental observations (Fig 2D-2F, [3]), and, in certain simulations related to lateral root initiation, of an auxin importer positioned in all stele cells (Fig 2B, see discussion below). Transporter efflux (and, during lateral root initiation, influx) contributions are altered by the hormones on a cellular level, using the average concentration of auxin and cytokinin in each cell (Fig 2I). Simple sigmoidal response curves (Hill-functions) are used, either in an activating or in an inhibitory fashion (Fig 2G and 2H). (See Material and Methods for equations and parameters.) Auxin upregulates PIN3 [10] and, during lateral root initiation, AUX1 [40]. By contrast, when calculating the expression levels of PIN1 and PIN7, the concentration of both auxin and cytokinin have to be taken into account [3] in a multiplicative fashion. Note that cytokinin is not transported polarly. We begin with simulations in which cytokinin is considered to influx into and efflux out of all cells in the cross-sectional plane equally. PIN1 and PIN7 are upregulated by cytokinin, a process which is simultaneously downregulated by auxin, due to AHP6 activation, which inhibits PIN upregulation by cytokinin. Our model incorporates this as a direct inhibition by auxin of the PIN1 and PIN7 upregulation by cytokinin (see Eq (5) and further). Note that in order to achieve a parsimonious model, the exporter in the epidermis, cortex and endodermis is not regulated (Fig 2C).
Our simulations include a non-specified, generalised auxin importer to represent the activity of the AUX1/LAX importers. AUX1, LAX1, and LAX2 are expressed in the stele near the QC, while LAX3 is not (Fig 3, [47,48]); taken together, AUX1, LAX1, and LAX2 expression includes the entire stele. Based on these observations, we positioned the auxin importer in our simulations throughout the stele (Fig 2B) to represent the sum of AUX1, LAX1 and LAX2 activity. While AUX1, the best characterised auxin importer in Arabidopsis, was reported to be upregulated by auxin in roots based on a microarray experiment [40], later observations found no change in the expression of an AUX1 marker line following auxin treatment [48], raising doubts regarding the generality of auxin-regulation of the auxin importers. We therefore investigated the regulation of the importers in this developmental context by conducting a qRT-PCR on 2 mm long sections of Arabidopsis root tips treated with 1 μM indole-3-acetic acid (IAA). Expression of AUX1, LAX2, and LAX3 showed a slight but significant increase in response to auxin (S1 Fig). However, the increase was relatively mild compared with genes that are known to be auxin-responsive, such as IAA2, WOX5 (S1 Fig) or AHP6 [3]. We therefore chose to not include regulation of the auxin importer level in our model in general, even though the qRT-PCR results do not rule out the possibility of indirect auxin regulation of the importers or cell-specific changes in response to auxin. We have also conducted simulations in which the importer is regulated by auxin and found no significant differences in the results (compare S2  Fig and Fig 4F-4H). The observation that AUX1 is induced by auxin during lateral root formation [39,40] suggests that this is at least the case in more mature tissues. We therefore do include this regulation when specifically analysing the process of lateral root initiation.

Regulation of PIN expression is necessary to generate observed auxin patterns
To assess the importance of hormonal regulation in this context, we first considered the impact of a transverse pattern in transporter positioning without including the regulatory network. To do so, we conducted simulations in which the PINs were not hormonally regulated and were therefore constantly expressed at full strength (Fig 4A), which we refer to hereafter as 'static' simulations. In both the geometric (Fig 4)  This suggests that the known transporter localisation is sufficient to generate appropriate wild-type auxin patterns within the different layouts. However, the simulated and experimentally observed auxin patterns did not match when cytokinin-deficient (wol mutants) or cytokinin-treated roots were considered, even though in the model the PINs were again positioned according to their experimentally observed subcellular localisation Note that while localisation of the PINs was altered (according to observations reported in [3]) in the in silico cytokinin treatment and wol mutant, cytokinin levels and response were not changed. The auxin pattern in the static simulations of wol (Fig 4C, S3B Fig) was broader than the observed pattern ( Fig 4J, S3E Fig), with auxin accumulating in the endodermis as well as in the stele. Likewise, static simulations of cytokinintreated roots (Fig 4D, S3C Fig) showed auxin accumulation in a broad domain that corresponds poorly to the reported pattern ( Fig 4K, S3F Fig).
We therefore tested what patterns are predicted when we take into account the observations that PIN levels are dynamically regulated by auxin and cytokinins [3,52]. (For details, see Materials and Methods.) These 'dynamic' simulations, in which we included hormonal regulation of the expression level of the PINs (Fig 4E), generated auxin patterns which were more similar to experimental observations in all three conditions, particularly wol (compare the dynamic results in Fig 4F-4H with Fig 4I-4K, in contrast to the static results in Fig 4B-4D).
We next tried to understand why the dynamic regulation of PIN levels is important in this patterning process. An inspection of the in silico auxin concentration profiles in wol with and without regulation suggests that introducing the regulation may promote the stele to effectively act as a sink, thereby generating the experimentally observed auxin pattern. To verify this, we compared the flux pattern of auxin in the static and full dynamics simulations of wol roots. We decomposed the flux at each position into a radial (inwards/outwards) and an angular (clockwise/counter-clockwise) component (Fig 5A). The stele as a whole acts as a sink, with a negative (i.e., inwards) net flux in both the static and dynamic simulations (Fig 5B). Unexpectedly, the magnitude of the inward radial flux in the stele is slightly higher in the static simulations. Thus, the difference in the resultant patterns with and without regulation is not due to regulation causing or amplifying the stele to act as a sink.
We therefore explored whether the spatial pattern of the auxin flux might account for the difference between the static and dynamic simulations of wol. In general, the flux pattern is quite similar in static and dynamic simulations (S5 Fig). There is no significant difference in the direction of the flux under the two conditions, but the magnitude of the flux is uniformly greater in static simulations, since hormonal regulation reduces the activity of the stelar PINs in the dynamic simulation. As a result, the apoplast accumulates more auxin in the static than in the dynamic simulations, as is seen by a comparison of auxin levels under the two conditions ( Fig 5C and 5D, zoomed for easier comparison in Fig 5E and 5F). We propose that the increased auxin in the apoplast in static simulations is available for uptake by the endodermis, cortex, and epidermis, resulting in increased auxin levels in these layers. Diffusion in the apoplast thus provides an avenue for auxin to 'leak' out of stele in the static simulations.

Active auxin import is required for protoxylem specification
It has been shown that the pin1 mutant presents defects in protoxylem formation while pin7 does not [3]. We therefore simulated the pin1 and pin7 mutants in silico by completely removing either PIN1 or PIN7 from the model. Consistent with earlier experimental findings, auxin is severely reduced in the xylem axis in simulations of pin1 (S6A and S6D Fig) but not pin7, although the auxin maximum is diffuse in pin7 simulations (S6B and S6E Fig).
While the pin1 and pin7 mutants have been examined in this context [3], experimental data on the role of the importers is lacking. We therefore investigated the role of the auxin importers in xylem specification. Initial examination of the simulated DR5 pattern in the auxin importer mutant did not reveal major differences from wild type (S6C and S6F Fig). However, an important difference became evident when we plotted the concentration of auxin accumulating in each cell of the xylem axis in simulated roots (Fig 6A and 6B). Despite the spread of the auxin maximum in simulated pin7 roots and its retention in the axis in the auxin importer mutant, auxin levels in the xylem axis are higher in pin7 than in the importer mutant (Fig 6B). In our simulations, the xylem axis of the importer mutant has only one-third as much auxin as in the wild-type xylem axis. Since auxin promotes protoxylem specification, we hypothesised that the lower auxin levels in the importer mutant may result in protoxylem phenotypes. We therefore examined the protoxylem of Arabidopsis plants with mutations in the auxin importers to test our prediction that these mutants would show defects in protoxylem formation. Although the aux1-21, lax1 and lax2 mutants all exhibited defects in protoxylem formation, the defects were not significantly more frequent than in wild-type plants (Fig 6C). However, roots of the triple mutant aux1lax1lax2 had a significantly higher frequency of protoxylem defects (multiple pairwise Fisher's exact test, Hommel-corrected, p < 0.01). Consistent with this, expression of the auxin-regulated DR5rev::GFP marker was reduced in the aux1lax1lax2 mutant in 11/14 roots (Fig 6F and 6G). We next conducted a more sensitive assay by germinating the mutants at a cytokinin concentration which normally causes protoxylem defects in about 50% of wild-type plants [35] in order to determine whether a reduction in active auxin import would make these plants more sensitive to cytokinin. Under these conditions, lax1 plants showed no difference from wild type, but specification of protoxylem in lax2 was, surprisingly, resistant to cytokinin (Fig 6D; multiple pairwise Fisher's exact test, Hommel-corrected, p = 0.0118). Both aux1-21 and aux1-22 also showed significantly more protoxylem defects than wild-type plants (Fig 6D and 6E; multiple pairwise Fisher's exact test, Hommelcorrected, p = 0.0479 and 0.0431, respectively), demonstrating that the robust formation of an auxin signalling maximum in the xylem axis is dependent on active auxin import by AUX1, as predicted by the model.

Vascular patterning does not require a cytokinin gradient
So far, the model operated without generating a gradient of cytokinin in the stele, leading to the surprising implication that no spatial pattern in cytokinin levels is required for correct and robust auxin patterning in the stele. Nevertheless, two different mechanisms have been proposed by which cytokinin levels could be distributed in a graded fashion. First, rootwards transport of cytokinin via the phloem [52] raises the possibility that a gradient of cytokinin centred around the phloem position might play a role in transverse patterning. While earlier simulations suggest that this gradient is unlikely to act as a source of positional information [4], they were based solely on the connectivity network between cells, rather than the spatial structure of the cells and their connections involving the apoplast. Secondly, a recent modelling effort concluded that local cytokinin biosynthesis in the xylem is crucial for vascular patterning in Arabidopsis [41]. In this model, cytokinin synthesis in the xylem axis generates a gradient of cytokinin levels and response which organises the tissue by controlling cell divisions.
In considering the potential importance of certain cells (either phloem or xylem) as sources of cytokinin, we first asked what properties of movement would be required to generate a gradient in the cytokinin levels. The gradient produced by any freely diffusing substance depends on its diffusion coefficient, D, and rate of degradation, δ; the ratio of the two determines the characteristic length of the gradient, l ¼ ffiffi ffi D d p ; which is the distance at which the concentration has dropped to 1 e = , or about 37% of the concentration at the source. The characteristic length is a good indicator of the spatial range over which a gradient can be informative; its scale should more-or-less correspond to the size of the tissue under consideration. To form a cytokinin gradient with a characteristic length on the order of the Arabidopsis stele (50μm) while assuming a diffusion coefficient similar to auxin, a degradation rate of around δ = 0.24 s −1 is needed. Equivalently, when assuming a degradation rate similar to auxin, a diffusion coefficient of around D = 2.5 × 10 −3 μm 2 /s is needed. Both values, however, seem to be very unrealistic and are inconsistent with the known spread of cytokinin. It therefore seems unlikely that an informative gradient of cytokinin can be generated solely via diffusion.
Even though an informative gradient seems unlikely, we nevertheless examined whether such a gradient could play a role in determining the auxin accumulation pattern in our simulations. This was evaluated by allowing cytokinin influx or production to occur solely at either the phloem pole cells or the xylem axis cells. Even with an extremely high cytokinin degradation rate to guarantee a steep gradient centred around the position of the phloem poles (S7A Fig [4], these findings indicate that no spatial pattern in cytokinin distribution is required to generate an auxin maximum in the xylem axis, a spatial pattern in cytokinin levels that is anyway unlikely. However, while no variations in cytokinin levels are required, the simulations do present variations in cytokinin response and signalling (due to the emerging auxin patterning, as each cell independently responds to auxin, which suppresses the cytokinin response machinery).
While our simulations indicate that a pattern in cytokinin levels is not required, it remains possible that such a pattern is relevant for vascular development. We therefore examined the plausibility of a cytokinin gradient more thoroughly. Given the extreme parameter constraints required to form a significant gradient in cytokinin levels, we considered whether other biophysical properties of the tissue might be able to relax those constraints. For example, a higher diffusion coefficient might be compatible with an informative cytokinin gradient if cytokinin movement is obstructed by plasma membranes. To test this, we ran simulations in which cytokinin movement was affected by plasma membranes, rather than being able to freely pass through them, by introducing a bidirectional cytokinin membrane permeability. To further release the constraints on cytokinin movement required to form a gradient, we considered a lower diffusion coefficient for cytokinin within the apoplast, as is the case for auxin [53]. We varied both parameters to find a regime in which localised cytokinin synthesis could form an informative gradient within the stele via diffusion, as suggested by a recent study [41]. We tested localised cytokinin production in either the xylem axis or the phloem. Remarkably, both apoplastic diffusion and cytokinin membrane permeability must be set to unrealistically low values to generate a gradient. When cytokinin membrane permeability was set equal to the (low) value for auxin efflux when there are no PIN transporters along the plasma membrane, diffusion was still sufficient to generate a very even, almost flat cytokinin distribution (Fig 7  and S8 Fig, leftmost column). Generating a weakly uneven distribution required a membrane permeability of 10% of the auxin value (Fig 7 and S8 Fig, second column), while a permeability of at most 1% was necessary to generate a strongly uneven distribution (Fig 7 and S8 Fig, third and fourth columns). Even under these extreme conditions, the uneven distribution could not be called a 'morphogen gradient', because cytokinin was high in the source cells and evenly low elsewhere, without gradation in between (Fig 7 and S8 Fig, top row); the pattern is thus completely cell-autonomous. To generate a weak morphogen gradient, it was necessary to also reduce the diffusion coefficient in the apoplast to 1% of the symplastic value (Fig 7 and S8 Fig,  third and fourth column of middle row); a strong gradient required an apoplastic diffusion coefficient only 0.1% of that in the symplast (Fig 7 and S8 Fig, third and fourth column of bottom row). The reason is that the apoplast forms a contiguous space, causing cytokinin levels to even out over the tissue whenever transmembrane transport is very low. The key difference with auxin is that auxin transport can be polarly organised, allowing the slope of gradients to increase many orders of magnitude [54]. So far, there is no indication for such polar transport of cytokinin.
Until now, we have assumed a prepattern in cytokinin biosynthesis (or influx) in our efforts to explore the feasibility of forming an informative pattern of cytokinin levels. Alternatively, cytokinin biosynthesis could be regulated by auxin, generating a feedback loop that might enable the emergence of a formative gradient in cytokinin levels, as suggested by [55]. To explore this possibility, we implemented auxin-regulated cytokinin biosynthesis. Similar to the results of simulations with prepatterned cytokinin synthesis, auxin-driven cytokinin synthesis Severely reduced cytokinin movement is required to establish a cytokinin gradient across the Arabidopsis root via diffusion from a local biosynthesis source in the xylem axis. To establish a cytokinin gradient without altering the diffusion coefficient of cytokinin, we introduced bidirectional membrane permeability. Decreasing the membrane permeability (moving right within rows) concentrates cytokinin in the source cells (here, the xylem axis), but no morphogen gradient is established, because the diffusion within the apoplast evens out the distribution within the rest of the root. Decreasing the diffusion coefficient (downwards in columns) in the apoplast transforms the step-wise cytokinin distribution pattern into a true gradient. Parameter values are given in Table 6.  only manifested cytokinin patterning under extremely unlikely parameter settings. As soon as the apoplastic cytokinin diffusion coefficient or the cytokinin permeability rate was increased beyond an unreasonably low level, no informative cytokinin gradient was formed (S9 Fig).
The fact that none of our simulation attempts are capable of generating cytokinin patterning within reasonable parameters settings reflects a fundamental property regarding the spatial and temporal scales of morphogen gradients generated through localized production, undirected movement and decay. In such a system, the root mean square displacement, the distance which molecules are typically able to move before breakdown, is calculated as ffiffiffiffi ffi l. Thus, the distance molecules are typically able to move is directly related to the characteristic length, which means that within the context of source-decay mechanism fixing the characteristic length implies fixing the length of communication. Consequently, cytokinin molecules diffusing in a gradient with a characteristic length of 50μm would move, on average, around 71μm before breaking down, independent of the specific choice of effective diffusion or decay [54]. Thus, any choice of parameters to form a cytokinin gradient on this scale implies that cytokinins are able to move less than 100 microns (less than the typical length of a mature root cell), which would preclude a role for it as a signal in longitudinal root-shoot communication, contrary to observations [52,56].
Although severe restrictions on cytokinin movement seem unlikely based on our theoretical treatment and biological considerations of long-distance cytokinin communication, we also experimentally evaluated the validity of this reasoning. To do so, we treated Arabidopsis roots with cytokinin and measured the intensity of the fluorescent TCS::GFP cytokinin response marker [57] throughout the root cross section. TCS intensity increased in response to cytokinin within six hours, with a more marked increase seen following overnight treatment (Fig 8). Importantly, the increase was not confined to an initial response in the outer layers, but was observed at the same timescale throughout the root, implying that cytokinin was able to effectively penetrate into the stele within a timescale of hours, if not minutes. This is inconsistent with the very slow and limited cytokinin movement required in our model, reinforcing the biological implausibility of an informative cytokinin gradient forming on these spatial scales via an effective diffusion-driven process (note that an effective diffusion does not exclude that part of the observed spread might stem from cytokinin-triggered cytokinin biosynthesis).
Given the unlikelihood of an informative cytokinin gradient on the scale of radial patterning within the stele, along with experimental and computational data demonstrating that positional information from xylem-or phloem-derived cytokinin is not essential [4,52], we conclude that cytokinin most likely evens out, and we therefore continue to operate in the remaining simulations under realistic cytokinin biosynthesis and transport parameter settings, which consequently give rise to an even distribution of cytokinin. Under these settings, the results of simulations with local variations in biosynthesis (either due to a prepattern or due to feedbacks) or with homogeneously distributed biosynthesis become undistinguishable.

Polar PIN1 localisation generates an auxin circuit
While the observed subcellular localisation of PIN7 is clearly apolar, immunolocalisations of PIN1::GFP suggest polar localisation [3]. PIN1 is therefore polarly localised within the procambium cells in our simulations, positioned to transport auxin away from the phloem poles and towards the xylem axis ( Fig 2D). However, since the experimental data cannot rule out apolar PIN1 localisation, we also investigated the consequences of this possibility in our model ( Fig  9A). Apolar PIN1 localisation causes only a small change in the quantity of auxin accumulating in the axis (Fig 9B and 9C), but further examination revealed dramatic differences in the auxin flux pattern throughout the stele.  Polar localisation of PIN1 towards the xylem axis results in increased co-ordination of the fluxes through the procambium and pericycle (Fig 10A and 10B), which can be seen by a detailed examination of the flux patterns (S10 Fig). In simulations with apolar PIN1 localisation, the radial flux alternated between inwards and outwards along the circumference of the pericycle; by contrast, roots with polar localisation of PIN1 presented consistent inwards radial flux in the pericycle (despite spikes of outwards flux) except at the protoxylem poles (Fig 10A,  S10C Fig). Likewise, the direction of the angular flux within the pericycle presented strong variations in the simulations with apolarly localised PIN1, but was consistently directed towards the xylem axis in the simulations with polar PIN1 (Fig 10A, S10D Fig). In addition, in polar simulations angular flux was consistently towards the phloem poles in the endodermis ( Fig  10A, S10F Fig), but this was not the case in apolar simulations (S10F Fig). In the procambium, radial flux is inwards in simulations of roots with polar PIN1 (again, the protoxylem position is an exception), despite spikes of outwards flux in the cell walls, while roots with apolar PIN1 have a domain of outward flux in the procambium spanning over a quarter of the circumference (Fig 10A and 10B, S10A Fig). Unlike in the pericycle, the angular flux in the procambium does not show strong variations; however, simulations with polar PIN1 show four quadrants of angular flux, with each protoxylem pole receiving auxin via flux from two of the quadrants ( Fig  10A, S10B Fig). By contrast, simulations of roots with apolar PIN1 localisation only have two domains of angular flux, both oriented towards the same protoxylem pole (Fig 10B, S10B Fig).
The co-ordinated fluxes resulting from polar subcellular localisation of PIN1 generate an auxin circuit. Auxin flows towards the phloem poles and then inwards; in the procambium, it flows towards the xylem axis and then out through the protoxylem poles (Fig 10A). This circuit is absent in simulations with apolar localisation of PIN1 in the procambium (Fig 10B), which show outward flux at the protoxylem poles but disordered flux in the pericycle and procambium. The auxin circuit leads to higher levels of auxin accumulation within the axis (Fig 9D), despite having greater levels of flux within the stele. Despite these significant differences, both polar and apolar PIN1 localisations generate an auxin maximum in the xylem axis, prompting us to consider whether the auxin circuit may confer other advantages, such as increased robustness.

The auxin circuit confers robustness and flexibility
In order to explore whether the auxin flux circuit generated by polar PIN1 localisation confers additional robustness or other benefits to the developing root, we investigated the resilience of the auxin maximum to transient AUX1 activation (and thus auxin accumulation) in pericycle cells in simulations with either apolar or polar localisation of PIN1 and PIN7. Upregulation of AUX1 in xylem-pole pericycle cells is one of the earliest steps in lateral root initiation [39]. We replicated this in simulations by activating AUX1 in a focal pericycle cell for a period of 120 seconds in a transverse root section located shootwards of previous simulations, where lateral root primordia are primed [58] (Fig 11A, 11D, 11G and 11J). Since the PIN3 expression pattern changes quickly shootwards of the QC [3] and its regulation is not as well understood, we only included PIN1 and PIN7 in these simulations; furthermore, we introduced auxin-regulated AUX1 expression, but limited to only the pericycle cells. We also evaluated the effect of PIN polarity in the pericycle independently of the procambium.
Simulations in which PIN1 and PIN7 were apolarly localised in the pericycle consistently retained auxin accumulation in the xylem axis when an AUX1 pulse was provided to a pericycle cell, regardless of the polarity of the PINs in the procambium cells (S11 Fig). However, this was not the case when PIN1 and PIN7 were polar in the pericycle. When polar localisation of PIN1 and PIN7 in the pericycle was combined with apolar localisation in the procambium, a pulse of AUX1 activation was sufficient to stabilise AUX1 upregulation and auxin accumulation in a xylem-pole pericycle cell. In geometric cross sections, activation within any pericycle cell could trigger the response, while in realistic root layouts it was limited to a subset of xylem-pole pericycle cells. However, auxin always accumulated in the same pericycle cell in geometric layouts, regardless of which cell received the AUX1 pulse (Fig 11B, 11E, 11H and 11K). In realistic layouts, a subset of the xylem-pole pericycle cells accumulated auxin following an AUX1 pulse (S12E and S12H Fig), while the other xylem-pole pericycle cells could not compete with the xylem axis for auxin accumulation (S12B and S12K Fig). By contrast, introducing polar PIN1 and PIN7 localisation in the procambium resulted in two important differences: (i) any xylempole pericycle cell (but no other pericycle cell) can present auxin accumulation and stable AUX1 expression; and (ii) this is only triggered specifically in the cell that received the AUX1 pulse (Fig 11C, 11F, 11I and 11L). The same result was found for both geometric and realistic root layouts (S12 Fig). Thus, the coordinated flux pattern generated by polar PIN1 localisation provides a mechanism for auxin to accumulate specifically in the xylem-pole pericycle cells, whenever another mechanism is able to trigger an initial (transient) increase in AUX1 expression in those cells. Possible mechanisms that could trigger such an initial increase have been discussed elsewhere (see [59] for a review).
In Arabidopsis, lateral roots are initiated by the division of xylem-pole pericycle cells [1]; the initiation is triggered by auxin accumulation in the pericycle cells [60] and tends to alternate between opposite xylem poles [61]. The alternating left-right pattern of lateral root positioning is dependent on AUX1 [58]. To evaluate whether this behaviour adheres in our computational model, we conducted simulations of roots with polarly localised PIN1 and PIN7 in which AUX1 was simultaneously activated in opposite pericycle cells for 120 seconds in both geometric (Fig 12A, 12D, 12G and 12J) and realistic (Fig 13A, 13D, 13G and 13J) cross sections. In no case did AUX1 remain active in both pericycle cells following a pulse; in simulations with polar PINs, one of the two focal cells invariably retained AUX1 and thus accumulated auxin, while the other subsided (Fig 12C, 12F, 12I and 12L; Fig 13C, 13F, 13I and 13L). When PIN1 and PIN7 were apolarly localised in the procambium, simultaneous AUX1 pulses resulted in an auxin pattern similar to a single pulse; geometric cross sections consistently accumulated auxin in the same cell, while realistic cross sections only accumulated auxin in a subset of the xylempole pericycle cells (Fig 12B, 12E, 12H and 12K; Fig 13B, 13E, 13H and 13K).
To further illustrate the potential of a "triggered" xylem-pole pericycle cell to suppress other xylem-pole pericycle cells to get triggered, we introduced a delay between the application of the AUX1 pulse in the two cells. For a sufficiently long delay, the cell which received the pulse earlier consistently accumulated auxin at the expense of the other cell. For realistic cross sections, a delay of 5s was sufficient (S13 Fig); for the geometric cross section, a 100s delay was sufficient, while a 20s delay was too brief (S14 Fig). This indicates a bistability induced by the transverse auxin flows; the system ensures that either one or the other pole can accumulate auxin. A situation in which both poles accumulate auxin is unstable, with tiny initial biases related to AUX1 and auxin levels and tissue layout determining which cell will eventually present persistent accumulation.

Discussion
To date, efforts to understand auxin transport have focused primarily on the longitudinal dimension [31,[62][63][64][65]; recently, however, experiments have indicated a role for transverse auxin transport in vascular patterning [3], and computational research has begun to address this question [4,41]. In this work we have demonstrated that the experimentally derived model is sufficient to maintain an auxin maximum along the xylem axis. We have also identified a novel role for active auxin influx in this network, and discovered a connection between subcellular polarity and the emergence of robustness and flexibility at the level of the whole tissue. These insights into the dynamics of transverse auxin transport will play an important role in future attempts to construct three-dimensional models of auxin transport and root growth in Arabidopsis.
Our simulations indicate that the mutually inhibitory auxin-cytokinin loop [3] is capable of maintaining an auxin maximum in the xylem axis. In reality, however, the auxin maximum becomes gradually confined to the protoxylem, while our simulations consistently show sustained auxin accumulation throughout the xylem axis, akin to the pattern observed near the quiescent centre in Arabidopsis roots [3]. Since the regulatory dynamics captured in our model do not recapitulate this change in the auxin pattern, we suggest that other factors are required to drive it; for example, PIN3 expression is known to expand into the metaxylem cells in the same domain where auxin signalling disappears from them [3] and is a plausible candidate.
In addition, regulation of the expression level of PIN transporters by this network is required in order to produce the auxin pattern observed in several conditions, indicating that not only the localisation of the PINs but also the correct regulation of the strength of their expression plays an important role in auxin patterning. Our analysis of the auxin fluxes in simulations of wol-like roots with and without dynamic regulation of the PINs demonstrates that, despite identical localisation of the PINs, dynamic changes in their expression level lead to local differences in the auxin fluxes, causing differences in auxin distribution and flux at the tissue scale. In addition to its relevance to future computational models of auxin transport, this insight has implications for experimental biologists. If differences in expression level alone (i.e., without changes in PIN localisation) can lead to different patterns, this highlights the need to record not only the localisation of the PINs in any particular context, but also how strongly they are expressed and how this is regulated, which may involve changes in PIN expression or abundance (e.g., [66]).
While the addition of PIN regulation reconstitutes the observed auxin distribution in wol, dynamic simulations of cytokinin-treated roots do not show a good match with the experimental data, though they correspond more closely than static simulations. The simulations show unexpectedly high levels of auxin in the procambium, phloem, and pericycle cells. Several possible explanations for this discrepancy exist. Our simulations conflate AUX1, LAX1, and LAX2 into a single importer, ignoring differences in their expression pattern and auxin uptake ability. In fact, LAX2, the importer present throughout the stele (Fig 3K), takes up auxin at a lower rate than AUX1 and LAX1 [48], which have a more specific expression pattern (Fig 3I and 3J). As a result, our simulations likely overestimate the amount of auxin taken up by these tissues. Furthermore, while the simulations include cytokinin-induced changes in the expression of the PIN transporters (S4D-S4F Fig), the possible effect of cytokinin on the importers is not included. Cytokinin-induced regulation of the AUX1/LAX genes or changes in their expression pattern may be able to explain the persistent mismatch. Indeed, it has recently been shown that expression of AUX1 and LAX2 is reduced by cytokinin treatment [67].
Our simulations indicate that active auxin import plays an unexpected key role in accumulating a sufficient level of auxin in the xylem axis to result in protoxylem specification. We confirmed this prediction by observing a high proportion of roots with defective protoxylem in the triple auxin importer mutant, aux1lax1lax2; furthermore, protoxylem formation was found to be hypersensitive to cytokinin in the aux1 mutants. These findings are consistent with earlier studies, in which aberrant phyllotaxis and defective secondary xylem differentiation were found in double, triple, and quadruple importer mutants, but not in the single mutants [44,49], although more recent work has also reported discontinuities in the vascular strands of the leaves of lax2 mutants [48], suggesting that subtle phenotypes may be found in single importer mutants as well. In our analysis of the root tip, we unexpectedly found that protoxylem formation is resistant to cytokinin treatment in the lax2 mutant. This may be due to the downregulation of LAX2 by cytokinin [67]; lower LAX2 expression in procambium cells adjacent to the xylem axis may reduce their ability to compete with the xylem axis for auxin, allowing auxin to accumulate in the axis despite other changes induced by cytokinin. However, contrary phenotypes are observed in the root and leaf, suggesting that the picture may be more complex. It is also interesting to note that AUX1 and LAX3, but not LAX1 or LAX2, have been reported to play a role in lateral root initiation [48], while AUX1 and LAX1 or LAX2 (but not LAX3) have been reported to regulate phyllotaxis [44]. Here, we have identified a role for AUX1 and LAX1 or LAX2 (but not LAX3) in vascular specification in the primary root, suggesting the possibility that specific LAX genes may act in concert with AUX1 in different developmental contexts, amplifying and buffering the auxin pattern generated by the PIN transporters.
We found that an even distribution of cytokinin throughout the root is sufficient to activate the PINs in a way that results in efficient auxin transport to the xylem axis and patterning; a gradient centred around local sites of cytokinin biosynthesis, whether in the xylem or the phloem, is not required and is unlikely. Our results demonstrate the difficulty in generating a cytokinin gradient on the cross-sectional scale of the Arabidopsis root via an effectively diffusion-driven process. Forming a reliably informative gradient within the stele requires a membrane permeability of 0.01μm/s and an apoplastic diffusion coefficient of 0.6μm 2 /s. Altering these parameters results in either a steepening of the gradient to form a step-like cytokinin pattern (in which the cytokinin is effectively trapped in the source tissues) or a flat distribution of cytokinin throughout the root (when cytokinin can escape from the source cells).
To generate a graded distribution spanning the stele, it was necessary to significantly reduce the diffusion coefficient of cytokinin within the apoplast. In retrospect, the reason for this is clear. The source cells can retain high cytokinin levels when membrane permeability is low, but rapid cytokinin diffusion within the apoplast allows cytokinin to readily move to all of the remaining cells in the stele, which will take it up based on the membrane permeability. While lowering membrane permeability can trigger a difference in cytokinin levels, it is still the diffusion coefficient which determines whether the difference takes the form of a steep drop-off or an informative smooth gradient. In addition to demonstrating the implausibility of an informative cytokinin gradient forming on this scale via diffusive processes (whether or not combined with apolar transport), this also highlights the importance of explicitly including the apoplast in computational models of plants so as to avoid such artefacts. This finding is reinforced by spatial-temporal considerations of hormone diffusion within the context of this system: in order to be informative on the transverse scale of the Arabidopsis stele (around 50μm), a gradient would need to have a characteristic length of roughly the same order (Fig 14). Generating such a gradient via diffusion would require the degradation rate or diffusion coefficient of cytokinin to differ by five orders of magnitude from those used for auxin here (see Materials and Methods) and elsewhere [31,39]; such values have in fact been used to produce cytokinin gradients in other computational models [38]. However, such a large difference seems unlikely, given that auxin and cytokinin are molecules of approximately the same size, composition and shape. Furthermore, a consequence of kinetic parameter values that would give rise to a potentially informative transverse gradients would be that the root mean square displacement of cytokinin (the distance which molecules are typically able to move before breakdown) would be limited to approximately 75μm, which is inconsistent with the observation that cytokinin travels over long distances in planta [52].
In addition, the timescale on which cytokinin patterning takes place would also be four orders of magnitude slower than for auxin. We have shown that application of cytokinin at the surface of the root leads to activation of a cytokinin response marker within the stele, at a distance of approximately 50μm, within hours of treatment. We cannot rule out the possibility of a systemic response to exogenous cytokinin application; however, such a mechanism would seem to preclude the possibility of deriving positional information from cytokinin. While our data are not sufficient to provide a quantitative estimate of the cytokinin movement rate, they do establish a lower limit on the order of tens of microns per hour. This is on the scale of the upper limit of the parameters required to form an informative gradient in our simulations. Spatial patterning of cytokinin on a scale appropriate to affect vascular patterning in the root therefore seems unlikely, despite evident patterns in cytokinin signalling and response within these tissues. While cytokinin remains important in such a scenario, the patterning process is ultimately driven by auxin.
We therefore suggest that cytokinin is unlikely to form informative gradients in contexts at the scale of transverse vascular patterning in developing Arabidopsis roots. We propose that the observed spatial patterns in cytokinin signalling (Fig 1C) may instead result from patterning of the cytokinin response machinery, perhaps due to modifications by auxin signalling. Alternatively, these patterns may arise from local, directed cytokinin transport, potentially combined with localised production and/or degradation, a possibility which we cannot exclude. Regardless of the preferred scenario, our analysis raises questions and challenges regarding the mechanics and dynamics of cytokinin transport which must be further addressed by experiments and more modelling alike. Furthermore, given that the parameters used here and elsewhere [38] to generate a locally informative cytokinin gradient are incompatible with longdistance movement of cytokinin, our findings highlight how local patterning mechanisms may have implications at the scale of the whole organism. It is therefore important when constructing verbal or conceptual models to consider constraints which are introduced by conflicting requirements at different spatial scales. Quantitative considerations can lead to different qualitative mechanisms for how patterning occurs. Investigating the plausibility of a cytokinin gradient in quantitative terms has led us to a qualitatively different way of thinking about the role of cytokinin, where the hormone is spatially homogeneous and patterning emerges from the response machinery.
This contrasts with the perspective reflected in a recent computational study on the same network [41], which argued that local cytokinin biosynthesis is crucial for vascular patterning. Although the model in [41] distinguishes between cytokinin levels and response, the cytokinin response only serves to regulate cell division; the auxin/PIN dynamics are driven by cytokinin levels, not cytokinin response. A gradient in cytokinin itself is thus crucial for the operation of their model; indeed, they establish a cytokinin gradient by using a membrane permeability equivalent to the lowest value we used in Fig 7. However, their results rely on the biologically implausible requirement of extremely little to no cytokinin diffusion within the apoplast. Our findings therefore fundamentally contrast with those reported in [41], in which severely restricted cytokinin movement establishes a gradient in cytokinin levels which drives PIN localisation and thus auxin dynamics; the auxin pattern then modulates the cytokinin response to control cell division. Thus, a gradient in cytokinin levels is a driving element in their model. By contrast, while we concur that local biosynthesis may be important to ensure that enough cytokinin is available to the developing root, we argue that a cytokinin gradient is not required for the patterning, and it is, in fact, challenging to establish such a pattern on the cross-sectional scale of the Arabidopsis root. Although directed cytokinin transport would allow such a pattern to be established straightforwardly, there is currently no indication of such a process.
We note that similar constraints apply in the embryo, where spatial scales even smaller than in the root make the presence of an informative cytokinin gradient even less likely. Although the embryo is too small for the formation of an informative cytokinin gradient solely via diffusion, a cytokinin signalling pattern has been proposed as a source of spatial information in this context [41], raising the question of how this signalling pattern is established. A cytokinin gradient has also been proposed to play a role in patterning processes in other plant organs (e.g., [38]), and these proposals will likewise have to be re-evaluated in light of our findings. While cytokinin transport or patterned cytokinin perception may address these difficulties, this only pushes the patterning question back a level: how, then, are the patterns in perception or transport established?
Another outstanding question is how low cytokinin signalling levels are maintained in the metaxylem position. In the current work, auxin accumulates throughout the xylem axis, including the metaxylem, and inhibits cytokinin signalling. However, this inhibition is mediated via AHP6, which is expressed in the protoxylem but not in the metaxylem [35]. Although auxin signalling is seen throughout the xylem axis [3], consistent with our model, this alone is insufficient to account for the low levels of cytokinin signalling in these cells. The same assumption was used in a model of the auxin-cytokinin loop in the embryonic root [41], where auxin was considered to inhibit the cytokinin response in all cells, even those without AHP6 expression, such as the metaxylem. Implicit in these approaches is the assumption that an auxin-responsive gene acts in the place of AHP6 in the metaxylem to inhibit cytokinin signalling. This assumption was made explicit in another modelling effort, which introduced an additional, hypothetical negative regulator of cytokinin, suggested to be a CKX or type A ARR, to simulate a cytokinin response pattern which matched experimental observations [4]. However, in this case the putative cytokinin repressor was not regulated by auxin, but was directly positioned in the metaxylem cells. By explicitly including the activity and regulation of AHP6 and its unknown partner, future modelling efforts could test the effect of different regulatory mechanisms (e.g, differences in auxin sensitivity between AHP6 and its redundant partner), hormone treatments or mutant combinations; the predictions from these simulations would help guide experimental biologists to identify the missing components in this network.
In our simulations, the auxin distribution pattern within the maximum in the xylem axis results from an auxin circuit generated by the polar subcellular localisation of PIN1 and PIN7 in the procambium and pericycle. To date, the challenges in detecting PIN polarity in crosssections have impaired a definitive experimental determination of the subcellular localisation of the PINs in this context [3]. Our simulations have shown that polar efflux within the stele organises the auxin flux into a circuit directed towards the xylem axis. In addition to providing robustness to the auxin maximum along the xylem axis, this auxin circuit provides a mechanism favouring the xylem-pole pericycle cells, the very cells which will later initiate lateral root primordia. In the absence of the auxin circuit (i.e., if PIN1 and PIN7 are apolarly localised), geometric factors determine which pericycle cells (if any) are sensitive to increased auxin import leading them to accumulate auxin. For example, in our simulations of a geometric cross section with apolar PIN1 and PIN7, the auxin flux in the procambium is oriented towards a single protoxylem pole (Fig 10B) because of the geometry of the cross section, which contains a series of consecutive cell walls that are perfectly aligned from the central stele cell up to this cell. This bias is not seen in simulations with polar PIN localisation, in which the auxin circuit stabilises the xylem axis against such perturbations; increased auxin import persists in the xylem-pole pericycle cells but is transient in other pericycle cells. Thus, the auxin circuit not only buffers the xylem axis against fluctuations in auxin import in other cells, but also provides the flexibility needed for auxin to accumulate in the xylem-pole pericycle, an event which plays a critical role in lateral root initiation [39,60]. Furthermore, our results have shown that pericycle cells on opposite xylem poles compete for auxin accumulation, consistent with the observation that lateral roots are initiated in an alternating left-right pattern, never appearing concurrently at both sides [61]. Nevertheless, we cannot describe the factors which lead to one pole being favoured over the other in any given instance, although we have shown that the relative timing of activation is important. A full understanding of this phenomenon is likely to require incorporating existing longitudinal and transverse auxin transport models into a threedimensional model of the root.
In addition to their specific implications for development in Arabidopsis thaliana, our simulations have also uncovered a novel mechanism by which processes can be made sensitive to the geometry of a tissue. Such a sensitivity, which may be desirable in particular contexts, is dependent on the polar subcellular localisation of a transporter and can be modulated by the effect of other transporters and network components. In this context, we have shown how such a mechanism can alter the transport dynamics to generate a coordinated flux pattern at the scale of a tissue or organ. This circuit thus serves to link patterning at several levels, connecting subcellular polarity to tissue-level organisation. To establish a cytokinin gradient without altering the diffusion coefficient of cytokinin, we introduced bidirectional membrane permeability. Decreasing the membrane permeability (moving right within rows) concentrates cytokinin in the source cells (here, the phloem poles), but a gradient is not established because the diffusion within the apoplast evens out the distribution within the rest of the root. Decreasing the diffusion coefficient (downwards in columns) in the apoplast transforms the step-wise cytokinin distribution pattern into a true gradient. Parameter values are given in Table 6. (TIF) S9 Fig. Auxin-driven cytokinin biosynthesis also requires severely reduced cytokinin movement to establish a cytokinin gradient across the Arabidopsis root. All cells in the cross section can synthesize cytokinin in an auxin-dependent manner. As in the simulations with fixed cytokinin biosynthesis (Fig 7, S8 Fig), decreasing the membrane permeability (moving right within rows) does not establish a gradient because the diffusion within the apoplast evens out the distribution within the rest of the root. Decreasing the diffusion coefficient within the apoplast (moving downwards within the columns) transforms the step-wise cytokinin distribution pattern into a true gradient. Parameter values are given in Table 6 and Table 4.