Spatiotemporal modelling of hormonal crosstalk explains the level and patterning of hormones and gene expression in Arabidopsis thaliana wild-type and mutant roots

• Patterning in Arabidopsis root development is coordinated via a localized auxin concentration maximum in the root tip, requiring the regulated expression of specific genes. However, little is known about how hormone and gene expression patterning is generated. • Using a variety of experimental data, we develop a spatiotemporal hormonal crosstalk model that describes the integrated action of auxin, ethylene and cytokinin signalling, the POLARIS protein, and the functions of PIN and AUX1 auxin transporters. We also conduct novel experiments to confirm our modelling predictions. • The model reproduces auxin patterning and trends in wild-type and mutants; reveals that coordinated PIN and AUX1 activities are required to generate correct auxin patterning; correctly predicts shoot to root auxin flux, auxin patterning in the aux1 mutant, the amounts of cytokinin, ethylene and PIN protein, and PIN protein patterning in wild-type and mutant roots. Modelling analysis further reveals how PIN protein patterning is related to the POLARIS protein through ethylene signalling. Modelling prediction of the patterning of POLARIS expression is confirmed experimentally. • Our combined modelling and experimental analysis reveals that a hormonal crosstalk network regulates the emergence of patterns and levels of hormones and gene expression in wild-type and mutants.


Fig. S1
Trend in average root auxin concentration in wild type and mutants.            Methods S1 Using ImageJ to analyse experimental images.

Methods S2
Method for discretising the root and for implementing numerical simulations.
Notes S1 Comparison of modelled auxin concentration trend with experimental DII-VENUS data in the literature.

Fig. S2
Modelling results show that PIN and AUX1 auxin carrier proteins localise predominantly to the plasma membrane in the wildtype. All parameters are included in Table S1. (a, b) Localisation of PIN proteins to the plasma membrane at different root locations marked by red arrows. An example of PIN protein polarity is shown in (a). (c, d) Localisation of AUX1 proteins to the plasma membrane at different root locations marked by red arrows. profile (Werner et al., 2003; www.plantcell.org, Copyright American Society of Plant Biologists).
PIN1 proteins localise mainly to the vascular cell files but with a weak signal in the epidermal and cortical tissues, consistent with . PIN1 concentration data were collected for the vascular tissues, plotted and compared to profiles from the vascular and pericycle cell files from the model. The model profile is a plot of the average PIN1 concentration in each cell tier cross-section of the root rather than each grid point since the large difference between PIN concentrations in the plasma membrane and cytosol make grid point plots difficult to read (Fig. 7 in the main text).
The same process was followed for the PIN2 images. The data were captured from the 2 external regions and then combined. However, given the location of the PIN2 proteins in the lateral root cap, epidermis and cortical cells and the tapering root shape it was more difficult to capture data than for the PIN1 proteins in the central vascular region. Fig. S9 shows a) quantification of PIN2 protein by analysing the experimental images ( Fig. 1, ; and b) modelling predictions for the patterning of PIN2 protein.
There was a reasonable match in concentration trends in the wild type and PLSox but not for the other mutants. This may be due to the overly simplistic structure of the root model where the region of PIN2 expression in the model is an even rectangle of cells consisting of the 2 exterior cell files of equal width extending for the full length of the root. In reality, the region of PIN2 expression in the root is not a simple rectangle; it consists of the lateral root cap (which only extends to the EZ zone and is not included in the modelled root structure), and the epidermal and cortical cell files starting about 5 cell tiers above the QC region (Muller et al., 1998).

Fig. S10
Modelling results for PLSm transcription patterning in wild type. All parameters are the same as in Table S1.

Methods S1 Using ImageJ to analyse experimental images
In this work, experimental images were analysed using ImageJ (http://imagej.nih.gov/ij). With ImageJ, it is possible to define regions on an image (such as the vascular cell cylinder in Fig. A  below) and determine an estimate for the relative hormone response or protein concentration profiles by measuring the relative signal intensity along the selected region. Regions were defined by a series of consecutive rectangles which progressively diminished in size towards the distal end of the root tip. Relative signal intensity data were collected from each rectangle and concatenated to give a relative hormone response or protein concentration profile for the selected cell files derived from an original experimental image. For example, Fig. A below shows a PIN2 fluorescent image . Using ImageJ, the epidermal and vascular cell files were defined by rectangles, and the relative wt PIN2 concentration profiles for each cell file were extracted. The data acquired were plotted using MATLAB (Fig. A below).

Root structure
We defined a root structure following the Grieneisen et al. (2007) model (Fig. A below). The root is 10 cells wide with 4 epidermal cell files, 2 border/pericycle cell files and 4 vascular cell files, with 3 distal tiers of columella cells (Fig. A below). The root is 35 cells in length, including 3 tiers of columella cells, 12 tiers in the meristematic zone (MZ) and 20 tiers in the elongation zone (EZ). Cell wall thickness is 2 µm, and the width of all cells is 20 µm. Cell length varies along the longitudinal root axis. Cells in the columella and MZ regions are 28 µm long, including the cell walls. Since cortical cell lengths in the transition zone were shown to increase at a relatively constant rate from 28 µm to 64 µm over a root distance of approximately 280 µm in 10 d-old plants (Beemster & Baskin, 1998), the Grieneisen et al. (2007) root structure was modified by increasing cell lengths at the same rate from 28 µm at the proximal end of the MZ, by 15% per cell tier, to a length of 64 µm over 7 cell tiers into the EZ. The length of these 7 cell tiers was 28 µm, 32 µm, 36 µm, 42 µm, 48 µm, 56 µm and 64 µm, to give an overall root tip length of 1594 µm.
To facilitate the calculation of average concentrations in different parts of the root, 3 cell types were established. Cell type 1 was defined as the 4 epidermal cell files along the entire root length (i.e. including some columella cells); type 2 cells are the 2 pericycle files (including some columella cells); and type 3 cells are the 4 vascular cell files (including the quiescent centre and some columella cells).

Grid point representation of the root
In the model the root is divided up into 2 µm by 2 µm areas, each of which is represented by a grid point (GP , Fig. B below). Each GP can have different properties, depending on its location in the root, and is assigned a numerical property code. The root is therefore represented by a root map which is a 100 by 797 matrix of GP values. For modelling purposes, the plasma membrane

Shoot-root boundary
is included in cell wall GPs and the code for each cell wall GP can vary depending on the property of its associated plasma membrane. The properties assigned to each GP allow rules to be established governing such processes as species biosynthesis, decay, activation or inactivation, or species flux between adjacent GPs. Equations are set up to govern the flux between each GP and its 4 nearest neighbour (NN) GPs located to the N, S, E or W (Fig. B  below). For example all species can diffuse within the cytosol but only the hormones can cross the plasma membrane (PM) into the cell wall. Cytokinin (CK) and ethylene (ET) can diffuse across the PM while auxin crosses via PIN and AUX1 carrier proteins. Once within the cell walls, the hormones can diffuse between NN GPs. PIN and AUX1 proteins can cycle between the cytosol and the NN cell wall GP but cannot diffuse between adjacent cell wall GPs. Individual kinetic equations and parameter values for these processes are described in Table S1.

Fig. B
Grid point representation of a cell in the MZand example of flux between nearest neighbour (NN) grid points. Green, cell walls; yellow, cytosol; red line, plasma membrane; red dots, grid points in the cytosol; black dots, grid points in the cell wall; black arrows, diffusion between NN; red arrow, auxin cell influx mediated by AUX1 carrier protein; blue arrow, auxin cell efflux mediated by PIN carrier protein.

Root map
The root map consists of a digital matrix. Each point in the root map has a code which describes the property of that point within the root. An example of this (Fig. C below) shows the coding for a single tier of cells within the meristematic zone. The cytosol is represented by code 0 and the cell wall GPs by codes 1, 2 and 3. The plasma membrane is not represented as a separate entity within the root map but is modelled as part of the cell wall, with different plasma membrane properties (such as polar PIN placement) being included in cell wall properties. We adapted the approach of prescribed PIN protein placement used by Grieneisen et al. (2007). PIN protein localisation to the cell wall was modelled by specifying different rate constants (low, medium and high) governing PIN flux from the nearest neighbour cytosolic GP to the cell wall GP, depending on the individual cell wall codes 1, 2 and 3 respectively (kinetic equations and parameter values are included in Table S1). There are two additional cell wall codes defining PIN activity at the root-shoot cell wall boundary. For the vascular and pericycle boundary cell walls the rate constant for PIN transport is set to zero. Therefore, there is no root to shoot auxin efflux in these cell files. For the epidermal boundary cell walls, the root to shoot auxin efflux is directed by PIN proteins. PIN protein is cycled from the nearest neighbour cytosolic grid point to the cell wall grid points using different rate constants: 1 is low, 2 is medium, 3 is high (kinetic equations and parameter values for localising PIN proteins to cell wall are given in Table S1).

PIN and AUX1 protein localisation to the plasma membrane
The plasma membrane (PM) is not represented as a separate entity within the root map but is modelled as part of the cell wall, with individual plasma membrane properties being included in cell wall properties. Therefore for modelling purposes PIN and AUX1 proteins can be located in a PM/cell wall unified structure. The concentration of PIN protein at a PM/cell wall GP is determined by the cycling and recycling of PIN from the cytosol to and from the PM/cell wall and degradation of PIN within the PM/cell wall. There is no exchange of PIN proteins between PM/cell wall GPs. Cycling of PIN to a PM/cell wall GP is determined by the concentration of PIN at the nearest neighbour cytosolic GP (affected by cytosolic biosynthesis, degradation and flux) and the rate constant for localising PIN proteins to the PM/cell wall GP. The rate constant can be either low, medium or high (kinetic equations and parameter values for localising PIN proteins to cell wall are given in Table S1) and is specified by the cell wall code (1, 2 or 3 respectively) in the root map. The rate of recycling back from the cell wall to the cytosolic GP is determined by PIN concentration at the PM/cell wall, and is inhibited by the cytosolic auxin concentration (Paciorek et al., 2005), and kinetic equations and parameter values for localising PIN proteins to the PM/cell wall are given in Table S1.
AUX1 cycling from the cytosol to and from the PM/cell wall is similar to PIN cycling. It is determined by AUX1 concentrations at the cytosolic and cell wall grid points and two rate constants, one for flux in each direction. These rate constants do not vary between different cell walls. Again, AUX1 degrades within the PM/cell walls as well as in the cytosol and there is no exchange of AUX1 between cell wall GPs.

Flux at the shoot-root boundary
The shoot-root boundary cell wall is located at the proximal end of the root structure. Ethylene and cytokinin can diffuse between the cytosolic GPs of the most proximal cell tier and the boundary cell wall GPs. Auxin shoot to root flux, inhibited by downstream ethylene signalling, is assumed to occur in the pericycle and vascular cell files from the boundary cell wall GPs into the cytosolic GPs, but auxin efflux in the opposite direction is not allowed. Experimental evidence (Suttle, 1988;Chilley et al., 2006) indicates that a relatively high ethylene signalling response inhibits the transport of auxin from the shoot to the root tip. However, the molecular basis of this inhibition is unclear. We assume that a molecule or molecules, X, located downstream of ethylene signalling, inhibit the transport of auxin from shoot to root (Liu et al., 2010. Thus, at the shoot-root boundary cell wall, we assume the auxin flux from shoot to root is are constants, they can be mathematically grouped into one parameter (we have kept them separate in order to maintain the biological clarity). Therefore, the actual auxin concentration at the cell wall GP of the shoot-root boundary, , is not important on its own for calculating auxin shoot to root flux. In addition, in the model, shoot to root auxin flux is regulated by the concentration of X at the root cytosolic GP that is adjacent to the shoot-root boundary cell wall. Biologically, this implies that the ethylene downstream response at the shoot-root boundary is equal to or similar to the ethylene downstream response at the nearest neighbouring cytosolic GP in the root. This is a reasonable assumption as the two GPs are nearest neighbours.
For the most proximal epidermal cells, PIN-mediated auxin efflux from the cytosolic GPs into the boundary cell wall GPs is allowed but influx is not. The concentrations of ethylene and cytokinin in the boundary cell wall GPs are set such that there is a relatively smooth transition of concentration levels from the root to the shoot. Biologically, this implies that the concentration of ethylene and cytokinin at the shoot-root boundary is equal to or similar to their respective concentrations at the nearest neighbouring cytosolic GP in the root. This is a reasonable assumption as the two GPs are nearest neighbours.

Notes S1 Comparison of modelled auxin concentration trend with experimental DII-VENUS data in the literature
We have compared modelled auxin concentration levels with both IAA2::GUS data and DII-VENUS expression data that we have experimentally measured (see main text for detail). Here we further compare our modelling results with DII-VENUS data in the literature (Brunoud et al., 2012;Band et al., 2014). DII-VENUS is a fluorescent reporter expressed under the 35S promoter, which is localised to the nucleus and quickly degrades in the presence of auxin. It is considered that DII-VENUS fluorescence more accurately represents relative auxin levels (the inverse of DII-VENUS levels) than other reporters (Brunoud et al., 2012;Band et al., 2014). The DII-VENUS data in the literature (Brunoud et al., 2012;Band et al., 2014) were used to compare relative experimental auxin response levels with modelled auxin concentrations for different cell types in the root tip.
The root structure in our model includes epidermal, pericycle, vascular, columella and quiescent centre (QC) cells. Since the root architecture in the literature (Band et al., 2014) is more complex than the relatively simple structure in our model, we approximately classify the cell types in our model to match the cell types in the literature (Band et al., 2014), as summarised in the following To better compare the relative auxin levels for different cell types, we set the auxin level at the QC to be 1, for both experimental data (Band et al., 2014) and modelling results. Fig. A below shows the comparison between experimental data (Band et al., 2014) and modelling results. The relative auxin levels are plotted in descending order, so that experimental and modelling trends can be easily compared. The trend of the modelled auxin levels for 5 cell types (i.e., quiescent centre, stele, endodermis, epidermis meristem, and cortex meristem) is similar to the trend observed experimentally (Band et al., 2014). However, the epidermis and cortex in the elongation zone (marked by *) and the columella (marked by **) are markedly different. These discrepancies between modelling results and experimental observations could be explained as follows.
DII-VENUS expression level measurements rely on the rates of production and decay of DII-VENUS. DII-VENUS expression is under the control of the 35S promoter, while the rate of decay depends both on the levels of auxin and the auxin co-receptors TIR1 and AFB1-5 (Brunoud et al., 2012). Homogenous expression of 35S and the auxin co-receptors is therefore necessary to allow representative comparison of relative auxin levels using the DII-VENUS reporter. It was reported that DII-VENUS fluorescence was ubiquitous in tir1 afb1 afb2 afb3 quadruple mutant roots, which were also significantly less sensitive to auxin (Brunoud et al., 2012).The co-receptors TIR1, AFB1 and AFB3 were shown to have very low relative expression (as fusion proteins) in the columella and lateral root cap cells (Fig. S5, Brunoud et al., 2012). This could result in the underestimation of relative experimental auxin levels in these cell types and account for the difference between the Band and model results for the columella. Using a non-degradable reporter, mDII-VENUS, the 35S promoter was shown to have significantly increased expression in the transition and elongation zones of the epidermis and cortex (Fig S6A,B in Brunoud et al., 2012), which again could result in underestimation of relative experimental auxin levels in these areas. The DII-VENUS data in the literature (Band et al., 2014) suggest that auxin responses in the epidermis and cortex elongation zone are higher than all other regions of the root except the QC. Moreover, by taking into account the effect of increased 35S expression, auxin responses in the epidermis and cortex elongation zone could possibly exceed those in the QC. This result is not apparent in experimental imaging using the auxin reporter IAA2::GUS (Fig. 3). This discrepancy leads to a number of possibilities. First, the non-degradable reporter, mDII-VENUS, may not fully reflect 35S expression levels in the elongation and transition zones. Second, other unknown factors may suppress DII-VENUS in these zones. Third, there may be additional suppressors of the IAA2::GUS reporter in the transition and elongation zones resulting in variable reporter sensitivity to auxin in different regions of the root. As noted in Band et al. (2014), the higher auxin response derived using DII-VENUS data in the elongation zone brings into question the hypothesis that a gradual decrease in auxin levels from the QC maximum determines root developmental zones Grieneisen et al., 2007).
After taking the above factors into account, our modelling results are in reasonably good agreement with experimental results derived using DII-VENUS data (Band et al., 2014).
In addition, we further compare the modelled auxin concentration level with experimental measurement using DII-VENUS for three cell types (Brunoud et al., 2012) in the meristem proximal to the QC. Since the root architecture in the literature (Brunoud et al., 2012) is more complex than the relatively simple structure in our model, we approximate the classifications of the cell types in our model in terms of the cell types in the literature (Brunoud et al., 2012), as summarised in the following table.
Cell type in the literature (Brunoud et al., 2012) Corresponding cell type in our model stele vascular, meristem endodermis and pericycle pericycle, meristem cortex and epidermis epidermis, meristem Our modelling result (Fig. B below) shows that vascular, pericycle and epidermal cells have a high, medium and low auxin level, respectively. This trend is in agreement with experimental observations (Fig. 2B in Brunoud et al., 2012). In summary, we have compared our modelling results with experimental DII-VENUS data from the literature (Brunoud et al., 2012;Band et al., 2014). Our modelled auxin concentration trends are in reasonable agreement with experimental observations.
While the modelling equations must be formulated in specific forms for describing the kinetics of all processes as described in Fig. 2, many parameter sets can be fitted against experimental data as the number of parameters is much more than that of experimental observations. By examining parameters randomly, we find that, when a parameter changes, if we allow at least one or more other parameters to change, we can find a new set of parameters that meet the criteria for model fitting (see section 'Model fitting reveals that both PIN and AUX1 activities must be restricted to certain ranges in order to generate correct auxin patterning'). The model using the new set of parameters also makes correct predictions, leading to the conclusions drawn in this work. For example, we have systematically examined all results in this paper for three randomly-generated sets of parameters that meet the criteria for defining a wildtype, and find that all modelling results for the three sets of parameters are qualitatively similar. The set of parameters we have used to generate all results in this work is included in Table S1.
We have subsequently evaluated model sensitivity to changes in the parameter values defining the wildtype root (Table S1). The parameters used for sensitivity analysis include diffusion rates, auxin decay rate, cytokinin biosynthesis rate and AUX1 and PIN transcription rates. The valuation of model sensitivity has shown that the model is robust to variations in the parameter values, as detailed below.
Changes to the auxin diffusion constant have little effect on average auxin concentrations and auxin patterning in the root. If the auxin diffusion constant is very small (< 20 µm 2 s -1 ), the auxin transport rate is reduced and auxin patterning no longer emerges. Changes to the cytokinin and ethylene diffusion constants have an even smaller impact on auxin concentrations and patterning than changes to the auxin diffusion constant. The PIN protein diffusion constant affects the ratio of PIN proteins in the plasma membrane to that in the cytosol. At very low PIN diffusion constant values (<20 µm 2 s -1 ), PIN proteins predominantly localise to the cytosol, limiting auxin transport such that no auxin maximum emerges. As the PIN diffusion constant increases, PIN at the plasma membrane increases, which in turn increases the rate of auxin efflux and an auxin maximum emerges. Auxin patterning can emerge for a wide range of PIN diffusion constants (>=20 µm 2 s -1 and <= 1000 µm 2 s -1 ). Changing the AUX1 diffusion constant has similar effects and the formation of the auxin maximum is very robust to changes in the AUX1 diffusion rate (>=5 µm 2 s -1 and <=1500 µm 2 s -1 ).
The coordinated rates of expression of PIN and AUX1 proteins control auxin transport through the root. The relative expression level of the two transporter proteins determines the proportion of auxin in the cytosol to the cell walls, with high relative PIN expression driving auxin into the cell walls and high relative AUX1 expression concentrating auxin in the cytosol. The importance of the coordinated permeability of the two transporter proteins is also demonstrated by simultaneously increasing PIN and AUX1 protein expression levels.
Due to the discrepancy between modelling results and experimental observations for cytokinin patterning, it is important to evaluate the model sensitivity to the parameter changes relating to cytokinin. Auxin patterning proved to be very robust to changes in all parameters relating to cytokinin diffusion, biosynthesis and degradation.
Model sensitivity analysis also shows that auxin patterning always emerges if the auxin decay rate constant is reduced to 50% or increased to 200% of the reference wildtype value.
In summary, model sensitivity analysis shows that the modelling results are robust to variation in parameter values and that there is a wide range of values that meet the criteria for generating correct auxin patterning.