Mathematical Modelling of Auxin Transport in Plant Tissues: Flux Meets Signalling and Growth.

Plant hormone auxin has critical roles in plant growth, dependent on its heterogeneous distribution in plant tissues. Exactly how auxin transport and developmental processes such as growth coordinate to achieve the precise patterns of auxin observed experimentally is not well understood. Here we use mathematical modelling to examine the interplay between auxin dynamics and growth and their contribution to formation of patterns in auxin distribution in plant tissues. Mathematical models describing the auxin-related signalling pathway, PIN and AUX1 dynamics, auxin transport, and cell growth in plant tissues are derived. A key assumption of our models is the regulation of PIN proteins by the auxin-responsive ARF-Aux/IAA signalling pathway, with upregulation of PIN biosynthesis by ARFs. Models are analysed and solved numerically to examine the long-time behaviour and auxin distribution. Changes in auxin-related signalling processes are shown to be able to trigger transition between passage- and spot-type patterns in auxin distribution. The model was also shown to be able to generate isolated cells with oscillatory dynamics in levels of components of the auxin signalling pathway which could explain oscillations in levels of ARF targets that have been observed experimentally. Cell growth was shown to have influence on PIN polarisation and determination of auxin distribution patterns. Numerical simulation results indicate that auxin-related signalling processes can explain the different patterns in auxin distributions observed in plant tissues, whereas the interplay between auxin transport and growth can explain the 'reverse-fountain' pattern in auxin distribution observed at plant root tips.

Introduction

Plant growth and development is tightly controlled by the spatial distribution of the plant hormone auxin. Auxin distribution patterns are organ specific (Petersson et al. 2009) and may be classified into two general types: spot and passage patterns. Spot patterns are characterised by local maxima of auxin concentrations and are observed in primordium initiation of leaves and flowers, as well as formation of lateral roots (Benková et al. 2003; Dubrovsky et al. 2008). High concentration of auxin at these points promotes cell growth and division, leading to organ development. Passage patterns are characterised by files or networks of (neighbouring) cells which have higher auxin concentrations than those surrounding them and are observed principally in the leaves and roots. In developing leaves, auxin distribution becomes arranged in a passage pattern forming networks, and the leaf veins are formed along these networks (Biedroń and Banasiak 2018).

Auxin transport and distribution in a plant tissue are controlled by the family of auxin-efflux carrier protein PIN-FORMED (PIN) (Leyer 2005) and auxin-influx carrier AUXIN RESISTANT1 (AUX1) (Yang et al. 2006). PIN is necessary for the formation of heterogeneous auxin distributions observed in plants (Okada et al. 1991). PIN proteins are localised to the plasma membrane of cells where they are then responsible for active transport of auxin out of the symplast. While it is clear that some form of feedback mechanism exists that links auxin to the polarisation of PIN (Chen et al. 2014, 2012; Robert et al. 2010), the exact nature of this feedback remains unclear (Feng and Kim 2015; Gao et al. 2015). One key hypothesis for the mode of the feedback mechanism is chemically via a so-called canalisation effect, where auxin flux through a cell membrane has a positive effect on PIN localisation to that membrane; however, there is also evidence for a strain-based mechanism (Homann 1998), where PIN is localised to the membranes with higher mechanical strain. Differential expression of AUX1 is also required for auxin pattern formation in some tissues (Swarup et al. 2001); however, most cells have symmetric distributions of AUX1 (Kleine-Vehn and Friml 2008).

The dynamics and transport of auxin in a plant tissue are also regulated by an auxin-related cellular signalling pathway. Auxin influences gene expression via the so-called ARF-Aux/IAA signalling pathway (Lavy and Estelle 2016). The signalling pathway describes a mechanism where auxin influences the levels of the family of gene transcription factors AUXIN RESPONSE FACTOR (ARF), via an interaction with Aux/IAA transcriptional repressors. Thus, auxin modulates gene response by controlling the levels of ARFs, through which it plays a role in primary root growth (Wan et al. 2018), root hair formation (Zhang et al. 2018), fruit growth and flowering. It has also been shown that the auxin-related signalling pathway has an influence on PIN dynamics by having roles in governing its biosynthesis (Paciorek et al. 2005; Vieten et al. 2005), degradation (Abas et al. 2006), and polarisation (Sauer et al. 2006). Despite the clear importance of the ARF-Aux/IAA signalling pathway however, it is likely that this mechanism alone is not enough to explain all auxin responses and details of other auxin-related signalling processes are emerging (Leyser 2018).

Although the interactions between auxin, PIN, and the auxin-related signalling pathway are essential for the transport and heterogeneous distribution of auxin in a plant tissue, which are necessary for growth and development of plants, the exact mechanism of nonlinear coupling between these processes is not yet completely understood. Thus, the use of mathematical models to investigate the validity of possible interaction mechanisms is important to better understand the dynamics and pattern formation in auxin distribution in plant tissues.

There are several results on mathematical modelling of auxin transport through plant tissues, assuming that auxin influences the polarisation of PIN proteins in cell membranes. The flux-based transport enhancement approach (canalisation), where flux of auxin out of the cell through the membrane has a positive feedback on the localisation of PIN to this membrane, has been used to generate realistic branching patterns observed in leaf vein formation (Feugier et al. 2005; Fujita and Mochizuki 2006) and has also been analysed in Feller et al. (2015), Stoma et al. (2008). When considering auxin transport through both apoplast and symplast, the auxin-dependent PIN distribution has been modelled by assuming that PIN proteins preferentially localise towards neighbouring cells with high auxin concentration. This approach was employed to generate spot-type patterns in auxin distribution observed in phyllotaxis (Heisler and Jönsson 2006; Jönsson et al. 2006) and auxin channels (Merks et al. 2007). Further models considering influence of external auxin sensors on PIN distribution have also had success in capturing passage patterns in solutions of mathematical models including the apoplast (Wabnik et al. 2010), although the biological relevance of this mechanism has been questioned (Feng and Kim 2015; Gao et al. 2015). The problem of generating different types of patterns in auxin distribution via unified mechanisms was addressed in Cieslak et al. (2015) by considering the notion of ‘unidirectional fluxes’ with a model based on Petri nets and in Hayakawa et al. (2015) where the influence of non-flux-based feedback of auxin on PIN polarisation was described by auxin-dependent PIN degradation. Both of these models demonstrated that a change in a single parameter could lead to switching between passage and spot patterns in auxin distribution in a plant tissue. Mathematical models have also been used to show how the distribution of auxin in the plant root tip is maintained (Band et al. 2014; Mironova et al. 2010). An excellent summary of various mathematical models of polar auxin transport may be found in van Berkel et al. (2013).

In this work, we derive and analyse novel mathematical models for nonlinear interactions between auxin-related signalling processes, PIN and AUX1 dynamics, intercellular auxin flux, and growth of a plant tissue. For our modelling, we primarily assume a flux-based mechanism of PIN localisation of a similar form as in Hayakawa et al. (2015), coupled with a detailed model of the auxin signalling pathway. We show that including the interplay between auxin-related signalling pathway and dynamics of PIN proteins in the mathematical model for auxin transport allows us to obtain both spot- and passage-type patterns in auxin distribution, depending on the values of the model parameter representing the rate of binding of PIN to auxin-TIR1. Using linear stability analysis, we determine the range of model parameters for which homogeneous patterns are stable. This analysis identifies possible mechanisms for the formation of heterogeneous auxin distributions in plants and possible interaction points between auxin and PIN responsible for homogeneous, spot and passage patterns, respectively. By considering model parameters that would generate oscillatory dynamics in auxin concentration in the model for auxin-related signalling pathway in a single cell, we show that the coupling between PIN dynamics, auxin transport, and cellular signalling processes can explain the formation of oscillatory auxin responsiveness observed in the basal meristem of plant roots (De Smet et al. 2007). Numerical simulations of the mathematical model for auxin transport, coupled with PIN dynamics, signalling processes, and auxin-dependent growth, suggest that cell growth can be one of the mechanisms underlying the formation of the ‘reverse fountain’ of auxin flow in plant root tip (Grieneisen et al. 2007). Modelling and simulations of interactions between auxin-related signalling pathway and apoplastic auxin transport demonstrate dependence of pattern formation on assumptions on the mechanisms of auxin flux between symplast and apoplast and PIN localisation to the cell membrane. The incorporation of auxin-related signalling processes, tissue growth, and strain-dependent PIN polarisation into mathematical models for auxin transport, analysis of oscillatory dynamics in auxin and PIN levels in plant tissues and of formation of ‘reverse fountain’ in growing tissues, and comparison between different mechanisms for auxin flux and PIN localisation are the main novel contributions of the modelling and analysis presented here.

Materials and Methods

It is observed experimentally that cellular auxin mediates the dynamics of PIN via its signalling pathway, whereas PIN regulates the heterogeneous distribution of auxin in tissues by controlling auxin flux between cells (Abas et al. 2006; Paciorek et al. 2005; Sauer et al. 2006; Vieten et al. 2005). It is further known that auxin influences the plant growth on the cellular and organ levels (Fendrych et al. 2018; Reinhardt et al. 2000).

In this work, we derive and analyse new mathematical models for nonlinear interactions between auxin flux, auxin-related signalling pathway, PIN and AUX1 dynamics, and plant cell growth. Incorporating the signalling and growth processes into mathematical models for auxin transport allows us to investigate the influence of cellular processes on the distribution of auxin in plant tissues.

Geometric Setting

In our models for auxin dynamics, a plant tissue is represented by a regular lattice of N cells of square shape, and equal size and dimensions, as shown in Fig. 1. In modelling auxin transport through a plant tissue, we shall consider two cases: (i) assuming direct interactions between neighbouring cells as in Fig. 1a and (ii) distinguishing between auxin dynamics in symplast and apoplast. In the second case, we split the apoplast (middle lamella and plant cell walls) so that each cell has an equal portion of apoplast surrounding it. Then, on a regular lattice the geometry of a plant tissue will be given by squares representing the cell inside, surrounded by four equal, regular trapeziums representing the apoplast, as shown in Fig. 1b. Similar geometric representations have been used in previous models (Wabnik et al. 2010).

Fig. 1

Schematics of the tissue geometry used for numerical simulations. a Simple geometry considering only intracellular space and cell membrane, with auxin flux considered to occur directly between cells. Here represents the volume of cell i, and represents the size of the portion of the membrane of cell i between cells i and j. b Schematics of a plant tissue where the domains representing the apoplast are equally divided between neighbouring cells, and passive auxin flux also occurs in the apoplast. Here represents the volume of apoplast compartment bordering cell i between cells i and j, and represents the size of the border between apoplast compartments (i, j) and (i, k)

Mathematical Model for Auxin-Related Signalling Pathway

In plant cells, auxin is perceived by the TRANSPORT INHIBITOR RESPONSE 1 (TIR1) receptor protein (Dharmasiri et al. 2005; Kepinski and Leyser 2005). Upon perception by and direct binding to TIR1, auxin enhances the interactions between TIR1 and Aux/IAA by acting as a ‘molecular glue’ (Tan et al. 2007), and the enhanced interaction between TIR1s and Aux/IAAs leads to the degradation of Aux/IAAs (Salehin et al. 2015; Wang and Estelle 2014), see Fig. 2. When auxin concentrations are low, Aux/IAAs repress activity of AUXIN RESPONSE FACTOR (ARF) by directly binding to ARFs and inhibiting their transcriptional ability. Correspondingly, a rise in auxin levels leads to degradation of Aux/IAAs, releasing repression of ARFs by Aux/IAA. ARFs enhance the transcription of auxin-responsive mRNAs, including Aux/IAA, and may interact with the binding site as single monomers, dimers, and two monomers simultaneously. ARF binding sites may also interact with ARF–Aux/IAA complexes and with ARF and Aux/IAA as single molecules. Hence in the mathematical model, in each cell of a plant tissue we consider production and degradation of auxin, its binding to TIR1, and dissociation from TIR1. We also consider production of Aux/IAA from mRNA, Aux/IAA binding to and dissociation from auxin-TIR1, degradation of Aux/IAA from the Aux/IAA–auxin–TIR1 complex, and binding to and dissociation from ARF. We further consider dimerisation of ARF monomers and splitting of ARF dimers into ARF monomers. Finally, we consider mRNA transcription to be enhanced by ARFs as monomers, dimers, and double monomers, and inhibited by Aux/IAA as a single molecule interfering with an ARF monomer, and as the ARF–Aux/IAA complex. The total concentrations of TIR1 and ARF are considered to remain constant.

Fig. 2

Schematic of the auxin-related signalling pathway. We assume PIN interacts with the auxin signalling pathway similar to Aux/IAA-type auxin-response proteins, and hence, PIN is degraded due to activation of the auxin-related signalling pathway (Color figure online)

In Abas et al. (2006), it was shown that auxin influences the degradation of PIN proteins via a mechanism similar to that by which auxin influences the degradation of Aux/IAAs. Auxin has also been shown to enhance PIN biosynthesis by controlling its gene expression through the ARF-Aux/IAA pathway (Paciorek et al. 2005; Vieten et al. 2005). From this limited evidence, we will assume in our model that auxin feedback on PIN biosynthesis operates via the ARF-Aux/IAA signalling pathway, specifically that ARF upregulates mRNA encoding PIN. We will assume that the mRNAs of Aux/IAA and PIN are identically regulated and that PIN binds to the auxin–TIR1 complex whereupon it may be marked for degradation. Thus, low levels of auxin lead to auxin transport being inhibited due to repression of PIN biosynthesis, medium levels of auxin lead to increase in auxin transport due to enhanced biosynthesis of PINs, and high levels of auxin lead to its transport being inhibited due to enhanced degradation of PIN proteins, as shown in Fig. 2. Hence in the mathematical model, in each cell of a plant tissue we consider production of PIN from mRNA, association of PIN to auxin-TIR1 and its dissociation from auxin-TIR1, and PIN degradation from the PIN–auxin–TIR1 complex.

Assuming spatial homogeneity of signalling processes in each cell, the dynamics of auxin signalling pathway on the cell level can be described by a system of ordinary differential equationscompleted with initial conditions given by initial concentrations of signalling molecules, specified in the analysis and numerical simulations of the mathematical model below, where the subscript i denotes to which cell the variable belongs, , and N is the total number of cells. Here, mRNAs are denoted by , cytosolic PIN is denoted by , auxin is denoted by , Aux/IAA is denoted by , TIR1 is denoted by , auxin–TIR1 complex is denoted by , PIN-auxin-TIR1 is denoted by , Aux/IAA-auxin-TIR1 is denoted by , ARF monomers are denoted by , ARF–Aux/IAA complexes are denoted by , and ARF2 dimers are denoted by . A list of all variables considered in our models can be found in Appendix Table 1. Model (1) is similar to the model for auxin signalling pathway derived in Middleton et al. (2010), with the inclusion of PIN as a secondary auxin-response protein.

Table 1

Listing of variable in Eqs. (1)–(8)

Variable	Meaning
m_i	Cellular mRNA concentration
r_i	Cellular Aux/IAA concentration
s_i	Cellular TIR1 concentration
c_i	Cellular auxin-TIR1 concentration
v_i	Cellular Aux/IAA-auxin-TIR1 concentration
e_i	Cellular PIN-auxin-TIR1 concentration
f_i	Cellular ARF concentration
g_i	Cellular ARF-Aux/IAA concentration
w_i	Cellular ARF2 concentration
a_i	Cellular auxin concentration
p_i	Cellular PIN concentration
u_i	Cellular AUX1 concentration
A_ij	Apoplastic auxin concentration
P_ij	Membrane-bound PIN concentration
U_ij	Membrane-bound AUX1 concentration

Parameter is the rate of mRNA production, is the rate of mRNA degradation, is the ratio of ARF-dependent mRNA production to ARF2- and double ARF-dependent mRNA production, and , , , , and are the binding thresholds to the relevant binding site of ARF monomers, ARF dimers, ARF–Aux/IAA complexes, molecules of ARF and Aux/IAA, and two molecules of ARF. The rate of Aux/IAA translation is , whereas and are the binding and dissociation rates of Aux/IAA and auxin-TIR1, and are the binding and dissociation rates of Aux/IAA and ARF, and is the degradation rate of Aux/IAA from Aux/IAA-auxin-TIR1. By and , the binding and dissociation rates of auxin and TIR1 are denoted, whereas and are the binding and dissociation rates of two ARF proteins, and are the binding and dissociation rates of PIN and auxin-TIR1, and is the rate of degradation of PIN from the PIN–auxin–TIR1 complex.

Auxin Transport in Plant Tissues

In the mathematical model for auxin transport in a plant tissue, we consider the dynamics of cellular auxin , the PIN-mediated flux of auxin between neighbouring cells, and the dynamics of cellular and membrane-bound PIN. The index ij denotes the membrane of cell i between two neighbouring cells i and j, e.g. denotes the size of the portion of the membrane of cell i between cells i and j.

Auxin is produced inside the cells with rate , degraded with rate , and transported between cells by membrane-bound PIN . Cellular PIN is translated from mRNAs with rate and its localisation to the cell membrane depends on the auxin flux through the membrane: stronger auxin flux through a specific membrane portion enhances localisation and leads to higher concentration of membrane-bound PIN in that part of the cell membrane (Fig. 3).

Fig. 3

Schematic of PIN-mediated auxin transport between two cells. Auxin (red circles) is transported from cell i to cell j by the efflux protein PIN (blue rectangles). In mathematical models, the concentration of auxin in cell i is denoted by , and the concentration of PIN localised to the portion of cell i’s membrane which neighbours cell j is denoted . The flux of auxin from cell i to cell j is denoted by and is assumed to positively feedback on the localisation of PIN to membrane portion ij between cells i and j (Color figure online)

Considering homogeneous distribution of membrane-bound PIN on each part of a cell membrane, see Fig. 1, the interplay between auxin flux and PIN dynamics is modelled by a system of strongly coupled nonlinear ODEswhere is short notation for and denotes the volume of the cell i. The flux of auxin between neighbouring cells i and j and the localisation of cytosolic PIN from cell i to membrane portion ij facing cell j together with dissociation of membrane-bound PIN back to the cell are given byHere H is a function describing the feedback of auxin flux on PIN localisation and is defined such that it is bounded between 0 and 1, increasing in , and . Parameter denotes the rate of PIN-mediated auxin transport, is the maximal rate of PIN localisation to the membrane, denotes the rate of PIN dissociation from the membrane, h is the flux-response coefficient, and is the flux threshold for positive feedback. In the response term H, flux is scaled by . The individual response terms were chosen such that for fluxes smaller (greater) than a threshold value , i.e. , the individual response would be approximately zero (one), ensuring strong positive feedback for large auxin fluxes.

Auxin-Dependent Tissue Growth

It has been observed that auxin can enhance cell growth in shoots (Reinhardt et al. 2000); however, auxin can inhibit primary root cell elongation (Overvoorde et al. 2010). Thus, in our model we consider the growth rate of a cell to be dependent on the concentration of auxin within the cell where denotes the length of either the horizontal or vertical wall of cell i, is the maximum growth rate, and is the threshold for half-maximal auxin-dependent growth rate. Here denotes that if is the horizontal (vertical) length of cell i then is the horizontal (vertical) length if cell j. The first Eq. (4a) corresponds to auxin-enhanced growth, whereas Eq. (4b) describes auxin-inhibited growth.

The growth of a cell is constrained by the cell wall and adhesion between cells leading to ‘tissue tension’, where slow-growing neighbouring cells will constrain growth of the cell, and fast-growing neighbouring cells will accelerate its growth. Hence, in our model we include a simple term for tissue tension such that growth rate of a cell is scaled by the ratio of the neighbouring cell length to the current cell length.

When considering signalling processes and auxin and PIN dynamics in a growing tissue, Eqs. (1)–(3) are modified by including the dilution effect due to growth:where () denotes the concentration of a chemical in cell i (membrane ij), and () denotes the reaction terms in the corresponding equations. Since in our model the cell shapes are simplified to be rectangular, is taken to be the cell length along the appropriate axis, and is the product of the length and width of the cell, where and are determined by (4) for the corresponding sides of the cell i.

Strain-Dependent PIN Localisation

There is evidence that plasma membranes undergoing higher strains have increased PIN localisation to them (Homann 1998). We model this mechanism by considering PIN localisation depending on the strain rate of the corresponding cell membranein addition to the auxin flux-related PIN localisation (compare with (3)), where is the strain-dependent rate of PIN localisation to the cell membrane. Here in the response term H flux is scaled by .

Symplast–Apoplast Model for Auxin Transport in Plant Tissues

Mathematical models considering direct flux of auxin between cells (see, e.g. Hayakawa et al. 2015; Stoma et al. 2008) provide a good framework to analyse the auxin transport through a plant tissue. However, along with active transport of auxin in/out of the cell it is important to consider the effect of passive flux of auxin through the apoplast. As described above, auxin is transported out of the cell symplast by membrane-bound PIN proteins. Due to the pH gradient between the apoplast and cytoplasm and weakly acidic nature of auxin, auxin passively diffuses from the cell wall into cell interiors; however, auxin is transported into cell symplast by membrane-bound influx proteins AUX1 at a much higher rate (Rubery and Sheldrake 1974; Yang et al. 2006). Auxin-influx protein AUX1 is synthesised within cells and then is trafficked to the cell membrane. Biosynthesis of AUX1 is known to be enhanced by auxin (Heisler and Jönsson 2006). Contrasting PIN, AUX1 is symmetrically localised in membranes for most plant cells (Kleine-Vehn and Friml 2008).

Thus, when considering both symplast (cell inside) and apoplast (plant cell walls and middle lamella), the mathematical model for auxin transport through a plant tissue, coupled with cellular signalling processes and dynamics of PIN and AUX1, in addition to Eqs. (1) and new equations for , , and , includes the dynamics of auxin in apoplast, cellular AUX1 , and membrane-bound AUX1 :where is the flux of auxin between cell i and ij-part of the apoplast, and denotes the localisation of AUX1 from cell i to membrane portion ij together with dissociation of back to the cell. Parameter is the translation rate of AUX1 from mRNA, and is the degradation rate of AUX1. A list of all variables considered in our models can be found in Appendix Table 1.

The transport of auxin across the ij-part of plasma membrane (part of membrane between cell i and ij-part of apoplast) combines active transport by PIN and AUX1 and a small contribution from passive diffusion. In apoplast, we consider passive diffusion of auxin between neighbouring apoplast compartments, denoted by and for different parts of apoplast, where denotes the size of the interface between apoplast compartments ij and ji, denotes the size of the interface between apoplast compartments ij and ik, and denotes the size of apoplast compartment ij. The passive fluxes of auxin through the apoplast and of AUX1 localisation to the membrane are given by:Here is the rate of passive flux of auxin through the apoplast, is the rate of AUX1-membrane localisation, and denotes the rate of AUX1-membrane dissociation. To analyse the emergence of patterns in auxin distribution and PIN polarisation in plant tissues, we shall compare two different types of auxin transport and PIN localisation: where features saturating auxin transport (Heisler and Jönsson 2006; Jönsson et al. 2006), is an extension of the flux considered in (3) by including the presence of apoplast, is a mechanism for PIN localisation, proposed in, e.g. (Heisler and Jönsson 2006; Jönsson et al. 2006), which, along with spontaneous localisation, specifies that higher auxin concentrations in neighbouring cells will cause PIN localisation to the membranes of the neighbouring cells, and is the mechanism for PIN localisation considered in (3), where J denotes the mechanism for auxin flux given by either (8a) or (8b), depending on which flux is considered in numerical simulations of the model.

Here is the rate of passive flux of auxin through the cell membrane, is the rate of PIN-dependent saturating auxin flux, is the rate of AUX1-dependent saturating auxin flux, is the rate of PIN-dependent non-saturating auxin flux, and is the rate of AUX1-dependent non-saturating auxin flux. Parameters , , denote the passive, PIN-dependent, and AUX1-dependent efflux of auxin, respectively, and , , denote the passive, PIN-dependent, and AUX1-dependent influx of auxin, respectively. Parameters , denote the concentration of auxin for half-maximal transport by PIN and AUX1, respectively. In localisation processes, is the rate of PIN membrane localisation, and is the rate of PIN membrane dissociation. Parameter denotes the proportion of PIN localisation that is auxin-dependent and is the half-maximal concentration of auxin for auxin-dependent PIN localisation. In , the dependence of PIN localisation on concentrations of auxin in neighbouring cells may be related to the fact that auxin-enhanced cell expansion places strain on the neighbouring membrane and thus enhances the PIN localisation (Homann 1998).

Numerical Methods and Implementation of Model Equations

Numerical codes for simulations of model Eqs. (1)–(3) or (1), (6)–(8) are implemented in Python, taking advantage of the Scipy module (Jones et al. 2001). Solutions were obtained using the scipy.integrate.odeint package which solves systems of ODEs using lsoda from the FORTRAN library odepack which can automatically select to use Adams (stiff) or BDF (non-stiff) methods, dynamically monitoring data to decide which method should be used (Hindmarsh 1983; Petzold 1983).

For numerical simulations, we consider two types of initial conditions: (i) small perturbations of the homogeneous steady state or (ii) zero concentrations for most molecules with the exception of TIR1 () and ARF () since the total amounts of TIR1 and ARF are conserved, and the conserved quantities were chosen as initial conditions. To calculate small perturbations around the homogeneous steady state, the homogeneous steady state was first calculated numerically, and then, in each cell the concentration of each component of the steady-state solution was multiplied by a random number between 0.9 and 1.1.

For certain simulations, we consider some cells to be either source or sinks. Compared to standard cells in the domain, in source cells the rate of auxin production is doubled and in sink cells the rate of auxin degradation is doubled. We solve symplast model (1)–(3) with initial condition as a perturbation of homogeneous steady state and periodic boundary condition to examine the emergence of spot and passage patterns in auxin distribution in plant tissue, as shown in Fig. 5. Zero-flux boundary conditions were considered to analyse the effect of boundary conditions on pattern formation. We solve model (1)–(4a) with the initial condition as a perturbation of the homogeneous steady state and periodic boundary condition to examine the effect of tissue growth on the emergence of spot and passage patterns, as shown in Figs. 8, 9, and model (1)–(4a),(5) to examine how varying the weighting between flux-induced and strain-induced PIN localisation affects auxin distribution in spot and passage patterns, as shown in Fig. 10. We solve symplast model (1)–(3) with the initial condition as a perturbation of homogeneous steady-state and zero-flux boundary condition to examine oscillatory dynamics in the auxin signalling pathway, as shown in Figs. 6, 7. We solve numerically model (1)–(4b) with zero initial condition and zero-flux boundary condition to examine the emergence of the reverse-fountain pattern of auxin distribution at the root tip, as shown in Fig. 11. For model (1)–(4b), (5), we consider zero initial condition and zero-flux boundary condition to examine how varying the weighting between flux-induced and strain-induced PIN localisation affects the emergence of the reverse-fountain pattern in auxin distribution at the root tip, as shown in Fig. 12. We compare solutions of model (1), (6), (7), (8a), (8c), model (1), (6), (7), (8a), (8d), model (1), (6), (7), (8b), (8c), and model (1), (6), (7), (8b), (8d), completed with periodic boundary conditions and initial condition as a perturbation of the homogeneous steady state, to examine the influence of different mechanisms for auxin transport and PIN localisation on patterns in auxin distribution, as shown in Fig. 13.

Fig. 5

Mathematical model including auxin signalling can generate both passage- and spot-type patterns of auxin distribution. Numerical solution of model (1)–(3) on a regular lattice of cells with . Green colour represents the concentration of cellular auxin, darker shades correspond to higher concentrations. Magenta colour represents the concentration of membrane-bound PIN, darker shades correspond to higher concentrations. a for low values of the rate of PIN binding to auxin-TIR1 , here , passage patterns of auxin distribution are formed. b for higher values of , here , spot patterns of auxin distribution can be formed, here blue borders indicate spots. All other parameter values are described in Appendix Tables 2 and 3. Periodic boundary conditions were used for both simulations (Color figure online)

Fig. 8

PIN polarisation aligns with oriented cell growth. Numerical solution of model (1)–(4a) on a regular lattice of cells, starting from the same initial conditions as in Fig. 5, but with cells undergoing auxin-dependent, horizontal growth. a In the passage parameter regime oriented growth has no effect on the placement of cells within the passage and only shifts the PIN alignment from vertical to horizontal for five cells. b In the spot parameter regime oriented growth disturbs the formation of spots, halting the emergence of one and enlarging another, and shifts the PIN alignment from vertical to horizontal for 36 cells. All parameters are described in Appendix Tables 2 and 3. Cells are represented in the reference configuration (Color figure online)

Fig. 9

PIN polarisation aligns with oriented cell growth. Numerical solution of model (1)–(4a) on a regular lattice of cells, starting from the same initial conditions as in Fig. 5, but with cells undergoing auxin-dependent, vertical growth. a In the passage parameter regime oriented cell growth has no effect on the placement of cells within the passage and only shifts the PIN alignment from horizontal to vertical for four cells. b In the spot parameter regime oriented growth disturbs the formation of spots, leading to the emergence of two new spots and modifying another into a small passage, and shifts the PIN alignment from horizontal to vertical for 35 cells. All parameters are described in Appendix Tables 2 and 3. Cells are represented in the reference configuration (Color figure online)

Fig. 10

Relative weighting of chemical and mechanical feedback on PIN localisation. Numerical solution of model (1)–(4a),(5) on a regular lattice of cells, starting from the same initial conditions as in Fig. 5a, but with cells undergoing auxin-dependent, horizontal growth, and with strain-dependent PIN localisation. a When chemical and mechanical PIN localisation are weighted equally, , four passage cells undergo small shifts and a total of 10 cells have altered PIN alignment, with one cell shifting from horizontal to vertical, and nine cells shifting from vertical to horizontal. b When chemical PIN localisation is dominated by mechanical PIN localisation, and , eight passage cells undergo small shifts and a total of 21 cells have altered PIN alignment, with three cells shifting from horizontal to vertical, and eighteen cells shifting from vertical to horizontal. All parameters are described in Appendix Tables 2 and 3. Cells are represented in the reference configuration (Color figure online)

Fig. 6

Oscillatory dynamics in the components of the auxin signalling pathway for appropriate rates of PIN membrane localisation. Numerical solution of model (1)–(3) on a regular lattice of cells, with zero-flux boundary conditions, and a parameter set listed in Appendix Table 2 that generates oscillatory dynamics in a single-cell model for auxin-related signalling pathway, a for PIN concentration is reduced to within physically realistic ranges and single spot cells with oscillating concentrations of auxin and cytoplasmic PIN are generated. b For the concentration of PIN on cell membranes rises above physically realistic ranges of M and numerical solutions display no oscillatory dynamics (Color figure online)

Fig. 7

Oscillatory dynamics are tissue-dependent and robust to growth. a Numerical solution of model (1)–(4a) on a regular lattice of growing cells with oscillatory parameters, see Appendix Tables 2 and 3, and zero-flux boundary conditions. Oscillatory single-cell spots are now joined by four-cell spots which also have oscillatory dynamics. b Numerical solution of model (1)–(3) on a regular lattice of cells, with zero-flux boundary conditions. Oscillatory parameters have been set for cells (5,4) and (5,7); all other cells have standard parameter values, see Appendix Tables 2 and 3. For modified oscillatory cells, oscillatory dynamics are not preserved (Color figure online)

Fig. 11

Influence of tissue growth on PIN polarisation contributes to the formation of ‘reverse-fountain’ auxin distribution patterns at the root tip. Numerical solution of model (1)–(4b) on a modified domain and different combinations of source and sink cells. a The central four cells in the top row are source cells, auxin flows from these cells to the base of the tissue with no flow from the base cells back up the tissue. b The outer four cells in the top row are sink cells, auxin flows from cells in the top five rows into these sinks, in the bottom three rows of cells auxin flows to the base of the tissue. c The central four cells in the top row are source cells and the outer four cells in the top row are sink cells, auxin flows from the source cells down the tissue, the central columns flow to the base of the tissue and the outer columns divert outwards to flow back up to the sink cells, resembling the reverse-fountain pattern observed at the root tip. Model parameters are described in Appendix Tables 2 and 3, with zero-flux boundary conditions. Cells are represented in the reference configuration (Color figure online)

Fig. 12

Strain-induced PIN localisation does not significantly affect the formation of reverse flows when weighted below chemically induced localisation. Model (1)–(4b),(5) solved on a modified lattice of growing cells, with source and sink cells as in Fig. 11. a Strain-induced PIN localisation is weighted equally with flux-induced PIN localisation, . Auxin flows from the source cells to halfway down the tissue where it branches out and then flows back up the outer cells. Auxin produced at the root tip also flows up the outer layer cells. b Strain-induced PIN localisation is weighted above flux-induced PIN localisation, , . Auxin flows directly from source cells to sink cells. Auxin produced in the central file of cells below the source cells does not flow to the base of the tissue, instead immediately flowing outwards to the outer cells where it then flows upwards to the sink cells. All parameters are described in Appendix Tables 2 and 3, with zero-flux boundary conditions. Cells are represented in the reference configuration (Color figure online)

Fig. 13

For the symplast–apoplast model, the mechanisms of auxin transport and PIN localisation determine the steady-state pattern in auxin distribution. a For model (1),(6)–(8)b),d) considering non-saturating auxin flux and flux-induced PIN localisation passage and spot patterns emerge similar to model (1)–(3). b For model (1),(6)–(8)a),d) considering saturating auxin flux and flux-induced PIN localisation similar patterns to the case with non-saturating auxin flux and flux-induced PIN localisation. c For model (1),(6)–(8)a),c) considering saturating auxin flux and concentration-induced PIN localisation a pattern of single-cell spots with high auxin concentration emerges. For model (1),(6)–(8)b),c) considering non-saturating auxin flux and concentration-induced PIN localisation the steady-state distribution is homogeneous (not shown). All parameters are described in Appendix Tables 2 and 4, with periodic boundary conditions (Color figure online)

Python code used to solve our models numerically is available online (Allen 2019).

Pattern Recognition

For two (neighbouring) cells A and B, auxin is defined as flowing from cell A to cell B if . We define two cells A and B as being connected if, potentially via some chain of other cells, auxin flows from A into B and from B into A, through different membranes. In the absence of any source or sink cells, we define a passage as a set of pairwise connected cells that splits the rest of the domain into two distinct sub-domains, for the three-cell case a passage occupies the entire domain, i.e. each cell is connected to every other cell. A spot is defined as either a single cell which experiences only auxin influx, or a set of pairwise connected cells that does not split the domain into sub-domains and experiences inflow from at least one neighbouring cell.

Linear Stability Analysis

Linear stability analysis of model (1)–(3) was performed to determine how changes in parameter values affect the auxin distribution in a plant tissue. The simplified domain considered in the stability analysis is a ring of three cells with each cell having two neighbours and communicating with every other cell in the domain. The outer boundary, i.e. cell borders which do not border other cells, has a zero-flux boundary condition. This is the smallest domain such that each cell communicates uniquely with every other cell. A similar approach was considered in Hayakawa et al. (2015).

We used standard numerical continuation techniques via the MatCont package in MATLAB (Dhooge et al. 2008), to examine the stability of steady-state solutions of model (1)–(3). This approach was used to determine the stability of the homogeneous steady-state solution of model (1)–(3) in order to locate the regions of the parameter space for which heterogeneous pattern formation is possible. We further used these methods to investigate the possibility of oscillatory solutions of model (1)–(3) upon variation of and .

To investigate the likelihood of occurrence of different pattern types when pattern formation can occur, model (1)–(3) are solved for various values of for different initial conditions defined as random perturbations of homogeneous steady state.

Results

Unless otherwise specified, parameter values are taken as the default values listed in Appendix Tables 2 and 3. To analyse the effect of auxin-related signalling processes on auxin transport and heterogeneous distribution in a plant tissue we considered model (1)–(3) in the three-cell ring domain for a wide range of parameters corresponding to the coupling between transport and signalling processes, i.e. the rate of PIN binding to auxin-TIR1 , and the sensitivity of the flux-feedback function h. The stability analysis shows that the increase in sensitivity of auxin-induced PIN degradation to auxin ( increases) leads to transition of stable heterogeneous solutions to homogeneous steady states, as shown in Fig. 4a. For sufficiently high values of h and appropriate values of , both spot- and passage-type patterns of auxin distribution are possible. Although both types of patterns were obtainable in the parameter region indicated in Fig. 4a, the value of has a great influence on the probability of each pattern emerging. To test the influence of parameters on probability of the emergence of specific pattern types, we used a lattice of 33 cells with periodic boundary condition. Specifically, we found that as increased the ratio of occurrences of passage patterns to spot patterns reduced dramatically, from as high as for to for , with fixed value . The type of boundary conditions also has an effect on the pattern formation for model (1)–(3). For zero-flux boundary conditions, the probability of the emergence of passage patterns was for and for . Since the value of h determines the range of values of for which heterogeneous patterns can emerge the probability of certain pattern types emerging for specific values of will also vary with h. Although and have similar roles in determining the degradation rate of PIN since they influence the rate of PIN-auxin-TIR1 binding, we found that upon varying from 0.5 to 50 for fixed values of between 1 and 250 the probability of emergence of specific pattern types was unchanged, with fixed value . Numerical simulation results for the cases of both passage and spot patterns in auxin distribution are included in Fig. 5 for . We tested a small set of values of , and for higher values of the concentration of membrane-bound PIN is increased leading to much higher concentrations of auxin in spot and passage cells compared to other cells and thus very well defined patterns, whereas for membrane-bound PIN levels resemble those in Fig. 5 (data not shown). When relating model (1)–(3) to symplast–apoplast model (1), (6)–(8) and considering in (3) concentration-induced PIN localisation as in (8c), spots formed of two cells that have strong PIN alignment between them emerge using parameter values that would lead to passages for (1)–(3) (data not shown).

Table 2

Parameter values for the model of auxin-related signalling pathway (1), calculated from non-dimensionalised values in Middleton et al. (2010)

Parameter	Description	Figs. 4a,b	Figs. 5a, 8a, 9a, 10a, 11, 12 (default)	Figs. 5b, 8b, 9b, 10b,13	Figs. 6, 7
α_m	Maximum mRNA transcription rate, μMmin^-1	0.5, 0–75	0.5	0.5	10
ϕ_m	Ratio of ARF-dependent to ARF2- and double	0.1	0.1	0.1	0.1
	ARF-dependent mRNA transcription rates
θ_f	ARF-DNA binding threshold, μM	1	1	1	1
θ_w	ARF2 binding threshold, μM	10	10	10	10
ψ_f	Double ARF-DNA binding threshold, μM2	0.1	0.1	0.1	0.1
θ_g	ARF + Aux/IAA-DNA binding threshold, μM	1	1	1	1
ψ_g	ARF-Aux/IAA-DNA binding threshold, μM2	0.1	0.1	0.1	0.1
μ_m	mRNA decay rate, min^-1	0.05	0.05	0.05	0.05
α_r	Aux/IAA translation rate, min^-1	5	5	5	5
β_r	Aux/IAA-auxin-TIR1 binding rate, μM^-1min^-1	5	5	5	5
γ_r	Aux/IAA-auxin-TIR1 dissociation rate, min^-1	5	5	5	5
μ_r	Aux/IAA decay rate, min^-1	5	5	5	5
α_p	PIN translation rate, min^-1	5	5	5	5
β_p	PIN-auxin-TIR1 binding rate, μM^-1min^-1	1–250	5	100	5
γ_p	PIN-auxin-TIR1 dissociation rate, min^-1	5	5	5	5
μ_p	PIN decay rate, min^-1	5	5	5	5
β_g	ARF-Aux/IAA binding rate, μM^-1min^-1	0.5	0.5	0.5	0.5
γ_g	ARF-Aux/IAA dissociation rate, min^-1	5	5	5	5
β_a	auxin-TIR1 binding rate, μM^-1min^-1	0.5	0.5	0.5	0.5
γ_a	auxin-TIR1 dissociation rate, min^-1	5	5	5	5
β_f	ARF dimerisation rate, μM^-1min^-1	0.5, 0.005	0.5	0.5	0.005
γ_f	ARF2 splitting rate, min^-1	5, 0.05	5	5	0.05
α_a	Auxin biosynthesis rate, μMmin^-1	0.5	0.5	0.5	0.5
μ_a	Auxin degradation rate, min^-1	0.5	0.5	0.5	0.5
S_tot	Total TIR1 present in cell, μM	10	10	10	10
F_tot	Total ARF present in cell, μM	10	10	10	10

Table 3

Parameter values for flux-based processes in models Eqs. (2)–(5)

Parameter	Description	Fig 4a	Figs. 4b, 6a, 7b	Figs. 5, 6b (default)	Fig. 7a	Figs. 8, 9, 10a,b	Figs. 11, 12a,b
ϕ_A	PIN-dependent auxin transport rate, μM^-1min^-1	0.005	0.005	0.005	0.005	0.005	0.005
λ	PIN membrane localisation rate, μmmin^-1	0.5	0.05	0.5	0.05	0.5, 0.5, 0.25, 0.1	0.5, 0.25, 0.1
ν	Strain-dependent PIN localisation rate, μmmin^-1	0	0	0	0	0, 0, 0.25, 0.4	0, 0.25, 0.4
δ_p	PIN membrane dissociation rate, min^-1	0.05	0.05	0.05	0.05	0.05	0.05
h	Sensitivity of PIN localisation to auxin flux	0-25	50	50	50	50	50
θ	Auxin flux threshold	2	2	2	2	2	2
χ	Maximal cell growth rate, μmmin^-1	0	0	0	1	1	1
θ_x	Auxin threshold for half-maximal growth rate, μM	0	0	0	0.8	0.8	5

Fig. 4

Analysis of parameter inference on solution types of model (1)–(3). a The boundary between pattern formation (shaded area) and homogeneous distribution (non-shaded area) is defined by an approximately linear relationship between h and . Minimum value of h presented here is 0.081. b As and are varied, model (1)–(3) undergoes a Hopf bifurcation and is able to have oscillatory solutions. For a single-cell model, i.e. , oscillatory solutions occur within the area bounded by the black solid line. For the three-cell geometry used for analysis model (1)–(3) are able to generate oscillatory solutions in the smaller area bounded by the dashed blue line (Color figure online)

Oscillations in the concentrations of targets of auxin-responsive ARFs have been observed experimentally in the protoxylem cells in the root meristem (De Smet et al. 2007). Considering the mathematical model for the auxin-related signalling pathway in a single plant cell, it has been shown in Middleton et al. (2010) that for certain parameter values solutions of the mathematical model can exhibit oscillatory dynamics. Here we demonstrated that our modified auxin signalling model (1)–(3) can have oscillatory solutions in the case when considering the dynamics in a single cell (with zero fluxes between cells) and in the case of PIN-mediated auxin transport in the three-cell domain, as shown in Fig. 4b. To analyse the effect of auxin transport on the oscillatory behaviour of component of the auxin signalling pathway inside cells in a tissue, we consider model (1)–(3) with the set of parameters for which oscillations in auxin concentration in the single-cell model would occur. We found that when considering oscillatory set of parameters for all cells in a tissue and reducing the rate of PIN localisation to the membrane by a factor of 10, i.e. , we obtain spot-type patterns in auxin distribution with some spot cells constituting oscillations in the levels of components of the auxin signalling pathway, which was not possible for the previously considered value of , as shown in Fig. 6. We tested a small set of values of , and observed oscillations in numerical simulations for but did not observe oscillations for . Interestingly, single-cell-type spots demonstrate oscillatory dynamics in auxin concentration, whereas four-cell-size spots do not present oscillations in the auxin concentration. When considering oscillatory parameters alongside horizontal growth however, (1)–(4a), not only are single spot cells with auxin concentration oscillating around a high value generated, but oscillatory dynamics are present in four-cell-size spots with smaller period, as shown in Fig. 7a. Assuming that different cells have different properties of the auxin signalling pathway, in two cells we considered the set of model parameters that would lead to oscillatory dynamics, with all other cells having standard parameter values, see Appendix Table 2. In this case, no stable oscillations in the auxin dynamics were observed and as the steady-state distribution of auxin we have low auxin concentrations in the two modified cells and high auxin concentration in the neighbouring cells which experience strong auxin flux from the modified cells, and all other cells have low concentrations of membrane-bound PIN (M), as shown in Fig. 7b. We further considered model parameters that would lead to oscillatory dynamics, i.e. as in Fig. 6, in a single-cell model in all cells within a specified radius from the central cell and standard parameter values, i.e. as in Fig. 5, in all other cells. In each case, apart from when every cell had parameters that would lead to oscillatory dynamics in a single-cell model, the oscillatory dynamics in auxin concentration were not persistent (data not shown). This suggests that in order to generate oscillations in the levels of targets of ARF observed in protoxylem cells, changes in the signalling pathways in all cells of a plant tissue are required.

The precise role of growth and its influence on the distribution of auxin in plant tissues is still open Korver et al. (2018). Incorporating the tissue growth in the auxin transport model (1)–(4a), where membrane-bound PIN is diluted on the growing membranes, we found that oriented, either horizontal or vertical, growth influences the overall PIN polarisation across the growing tissue, i.e. PIN preferentially polarises along the axis of growth; however, this effect does not alter the type of patterns in the auxin distribution in a tissue (i.e. passage or spot) but can alter its distribution, as shown in Figs. 8, 9. Interestingly, growth appeared to exert a stronger influence on the polarisation of PIN when considering parameter values that lead to the emergence of spot patterns in the absence of growth. We found similar results when considering horizontal and vertical growth simultaneously for identical parameters as in Figs. 8, 9 (data not shown). It seems that maximal growth rate has limited influence on the overall pattern formation; when was varied between 0.1 and 10 there were small changes in the exact concentrations, M for auxin and M for PIN; however, PIN polarisation patterns were unchanged. Incorporating strain-dependent localisation of PIN, see Eqs. (1)–(4a),(5), for equal or stronger weighting of flux-induced compared to strain-induced PIN localisation () we obtained similar patterns in the auxin distribution, as shown in Figs. 8a and 10a. However, when strain-induced localisation strongly dominates flux-induced localisation (), a significant reduction in concentrations of both auxin and PIN, compared to the cases where , is observed, as shown in Fig. 10. To ensure that maximum amount of PIN localised to a cell membrane is consistent with previous simulations, we considered a range of values of and such that . For , heterogeneous patterns do not form.

To analyse the effect of auxin-related signalling pathway and growth on auxin flux in the plant root tip, we consider model (1)–(4b) on a modified lattice of cells resembling a root tip. We assume that growth (cell elongation) occurs only along one axis, i.e. down the root, and is inhibited by high auxin concentrations and constrained by tissue tension. For consistency with auxin availability in the root, the bulk flow of auxin from shoot to root through the vascular bundle was simulated by including a source term in the central four of the top row of cells. In some simulations, we also included sinks in the epidermis cells, outer two cells on each side on the top row, since it is assumed that some auxin is evacuated from the root tip along the epidermis (Swarup et al. 2005). The steady-state solutions of model Eqs. (1)–(4b), considering zero initial conditions and no strain-dependent PIN localisation (), are presented in Fig. 11. When there are source cells only, auxin flows from the apex of the tissue to the base where it settles, as shown in Fig. 11a. When there are sink cells only auxin flows from cells in the top five rows of cells to the sinks, and auxin in the bottom three rows of cells pools at the base of the tissue, as shown in Fig. 11b. When both source and sink cells are included a reverse-fountain pattern similar to those observed at the root tip emerges where auxin flows down the root and both pools at the tip and branches out to flow back up the outer layers of cells, as shown in Fig 11c.

To analyse the effect of strain-induced PIN localisation to cell membranes on auxin flux in a plant root tip and its reverse flow, we considered Eq. (5) for a range of values of parameters (rate of chemical localisation) and (rate of mechanical localisation) such that so that maximum amount of PIN localised to a cell membrane is consistent with previous simulations. For , reverse flow patterns may still be generated, when auxin still flows from source cells to the base of the tissue as well as branching out and back up the outer cells, when auxin flows only partway down the tissue from the source cells before branching out to flow back up the outer cells and does not flow to the base of the tissue, as shown in Fig. 12a. For , formation of a reverse flow pattern was completely inhibited, with auxin flowing directly from source cells to sink cells, and auxin flowing from the base of the tissue to the sinks, as shown in Fig. 12b.

Table 4

Parameter values for the auxin flux and PIN and AUX1 dynamics considered in model Eqs. (6)–(8), parameters based upon those in Heisler and Jönsson (2006)

Parameter	Description	Fig. 13
α_u	AUX1 biosynthesis rate, μMmin^-1	5
μ_u	AUX1 degradation rate, min^-1	5
θ_u	Saturation of auxin-induced AUX1 biosynthesis, μM	1
ϕ_a	Auxin membrane permeability, μmmin^-1	0.55
κaef	Fraction of protonated auxin in cell	0.004
κain	Fraction of protonated auxin in wall	0.24
ϕ_p	Saturating PIN-induced auxin membrane permeability, μmmin^-1	0.27
ϕ~_p	PIN-induced auxin membrane permeability, μmμM^-1min^-1	0.27
κpef	Effective PIN-induced auxin efflux	4.67
κpin	Effective PIN-induced auxin influx	0.034
θap	Saturation of PIN-induced auxin transport, μM	1
ϕ_u	Saturating AUX1-induced auxin membrane permeability, μmmin^-1	0.55
ϕ~_u	AUX1-induced auxin membrane permeability, μmμM^-1min^-1	0.55
κuef	Effective AUX1-induced auxin efflux	0.045
κuin	Effective AUX1-induced auxin influx	3.56
θau	Saturation of AUX1-induced auxin transport, μM	1
ϕ_A	Rate of auxin diffusion in apoplast, μmmin^-1	67
ω_u	Rate of AUX1 localisation to membrane, μmmin^-1	0.5
δ_u	Rate of AUX1 dissociation from membrane, min^-1	0.05
ω_p	Maximum rate of PIN localisation to membrane, μmmin^-1	0.5
κ_p	Fraction of PIN localisation due to auxin feedback	1
θpa	Threshold for half-maximal auxin-dependent PIN localisation, μM	1
δ_p	Rate of PIN dissociation from membrane, min^-1	0.05

To analyse the effects of the mechanisms for auxin transport through apoplast and PIN localisation on the formation of auxin distribution patterns in a plant tissue, we solved model (1), (6)–(8) numerically on a regular lattice of cells with parameters as in Appendix Tables 2 and 4. When considering non-saturating auxin flux and flux-induced PIN localisation, i.e. mechanisms of the same form as in model without apoplast, (1),(6)–(8b),(8d), then behaviour is similar to the case without apoplast, with both passage and spot patterns able to emerge, as shown in Fig. 13a. When considering non-saturating auxin flux and concentration-induced PIN localisation, (1),(6)–(8b),(8c), then the steady-state auxin distribution is homogeneous (not shown). When considering saturating auxin flux and flux-induced PIN localisation, (1),(6)–(8a),(8d), then similar patterns emerge as in the case with non-saturating auxin flux and flux-induced PIN localisation, as shown in Fig. 13b. When considering saturating auxin flux and concentration-induced PIN localisation, i.e. mechanisms of the same form as considered in Heisler and Jönsson (2006), (1),(6)–(8a),(8c), then single-cell spots in auxin distribution emerge, as shown in Fig. 13c. When varying , the proportion of auxin-induced PIN localisation in (8c), between 0 and 1 and numerically solving (1),(6)–(8a),(8c), we found that heterogeneous spot patterns were only generated for values of , and homogeneous distributions were generated for (data not shown). For the simulations in Fig. 13, we used parameter values in the signalling pathway consistent with those used to generate spot patterns in model (1)–(3), since when using a smaller value of consistent with that used to generate passage patterns in model (1)–(3) the higher concentrations of PIN on cell membranes led to pooling of auxin in the apoplast compartments adjacent to high PIN-expressing membranes (data not shown). We performed numerical simulations for a set of values of and found that auxin pooling in apoplast compartments occurred for , whereas auxin concentrations were more realistic for . Disruption of heterogeneous auxin distribution and reduced uptake and accumulation of auxin is observed for the symplast–apoplast model for auxin transport (1), (6)–(8) when AUX1 is not included, , or PIN is overexpressed compared to AUX1, , (data not shown), agreeing with experimental observations (Okada et al. 1991; Yang et al. 2006).

Discussion

Until recently, one of the main criticisms of the canalisation hypothesis (auxin-flux-related localisation of PIN to cell membrane) was its inability to produce spot patterns in auxin distribution in a plant tissue without any additional assumptions on cell types (e.g. source/sink cells) (Stoma et al. 2008), despite its accurate capturing of passage patterns (Feller et al. 2015). Recent results (Cieslak et al. 2015; Hayakawa et al. 2015) have shown that it is possible to obtain both spot and passage patterns in auxin distribution considering the canalisation hypothesis, provided an extra mechanism of either auxin-mediated PIN degradation or auxin self-induced production is considered. This indicates the importance of intracellular processes for auxin transport in a plant tissue, which in the models in Cieslak et al. (2015); Hayakawa et al. (2015) were defined phenomenologically, without considering main biological mechanisms underlying those processes. Many of the mathematical models, e.g. in Hayakawa et al. (2015), and some results presented here, are restricted to periodic boundary conditions, which do not always provide best description of biological systems, however are useful for the analysis of auxin flux in homogeneous domains with no sources or sinks.

In the studies presented here we considered a detailed description of the auxin-related signalling pathway and its influence on PIN dynamics, with a key assumption of the regulation of PIN biosynthesis by ARF and PIN degradation by TIR1. This allowed us to identify that the rate of auxin-signalling-dependent PIN degradation, here represented by binding of PIN and auxin-TIR1, is key to determining the patterns of auxin distribution in plant tissues, as shown in Fig. 4.

For our model, we assumed that the mechanism of auxin-dependent degradation of PIN is similar to the degradation of Aux/IAAs via the auxin-TIR1 signalling pathway, despite limited biological evidence (Abas et al. 2006). Although we believe that including auxin-dependent PIN degradation is important for the results obtained by our model, it might be possible that this specific mechanism is not essential; for example, auxin-dependent PIN degradation has previously been modelled and shown to have a primary role in determining eventual auxin distribution pattern without the need to consider the full auxin signalling pathway (Hayakawa et al. 2015). Our model does not seem to be able to generate single-cell spots when considering flux-induced PIN localisation, and this is a key qualitative difference between our model and the model presented in Hayakawa et al. (2015), which is able to generate patterns composed of single-cell spots. Bifurcation analysis for the auxin transport and signalling pathway model presented in this paper suggests that it is possible to obtain spots and passage patterns for the same parameter values, where the likelihood of the emergence of passage patterns decreases and the likelihood of the emergence of spot patterns increases as sensitivity of PIN degradation to auxin increases. This differs from results obtained in Hayakawa et al. (2015) where a clear transition from passage-generating to spot-generating parameter regimes as sensitivity of PIN degradation to auxin increases is shown. It may be possible that adopting a similar mechanism of PIN degradation as in Hayakawa et al. (2015) in the model presented here would result in more similar bifurcation behaviour.

Our results on interactions between signalling and transport processes showed that the oscillatory dynamics in auxin concentration are obtained only when considering modified parameter values in the model equations for signalling pathway in all cells in the simulated tissue. This suggests that experimentally observed oscillations in auxin responsiveness are due to an oscillatory Aux/IAA negative feedback loop (Middleton et al. 2010) and that both the oscillatory feedback loop and PIN-mediated auxin transport through tissue are necessary for the formation of auxin distributions with local oscillations, as shown in Fig. 6. In cells other than those spots which have oscillatory dynamics, there are either damped or no oscillations, which is likely due to the very low concentrations of membrane-bound PIN in the membranes of oscillatory cells bordering non-oscillatory cells, i.e. the oscillatory dynamics of spot cells exert negligible influence on the dynamics of their non-spot neighbours. It would be of interest to investigate the dynamics of solutions of model (1)–(3) in the oscillatory parameter regime when solved on a realistic plant root geometry.

It has been observed that mechanical strain of a plasma membrane enhances PIN localisation to the corresponding membrane (Homann 1998). One model to consider such contributions was proposed in Hernández-Hernández et al. (2018), which modelled PIN localisation on the single-cell level using a discrete Boolean model, approximating continuous dynamical system, and predicted that mechanical forces could dominate molecular factors during PIN polarisation. Our numerical simulation results for the coupled auxin flux and tissue growth model (1)–(4a) indicate that mechanical forces could dominate the molecular activity since PIN is preferentially polarised, leading to the formation of auxin gradients along the axis of growth, as in Figs. 8, 9, whereas the strain-induced localisation of cytoplasmic PIN to the membrane, for membrane strain above a certain threshold, had a qualitative effect on the dynamics of auxin and PIN in a growing plant tissue. Numerical simulation results for mathematical model (1)–(4a), (5) also showed that the balanced contribution of chemical activities and mechanical forces to the PIN dynamics does not affect the type of patterns in the auxin distribution in a growing tissue.

Auxin is transported from shoot to root through the stele to the root tip where it is reorganised and then transported back up towards the shoot in the outer cell layers (Grieneisen et al. 2007). This directed auxin flux is commonly known as ‘reverse-fountain’ and has been observed to be essential for root development (Doerner 2008), for example in specifying the quiescent centre (Sabatini et al. 1999) and root responses to gravitropism (Swarup et al. 2005). Previous mathematical models described reverse flow in auxin patterns by prescribing polarisation of membrane-bound PIN (Band et al. 2014; Mironova et al. 2010; Stoma et al. 2008). Our new mathematical model for auxin transport in a plant tissue, that includes the dynamics of PIN coupled to the auxin-related signalling pathway, auxin flux, and tissue growth, is able to generate reverse flow patterns in the auxin distribution from an initial condition that does not have pre-established PIN polarity, as shown in Fig. 11. Our results suggest a plausible mechanism for the emergence of the ‘reverse-fountain’ auxin pattern observed at the root tip: the establishment of the PIN polarity that generates this characteristic auxin distribution is mechanically generated due to the dilution of PIN along growing membranes since when dilution is outweighed by strain-induced localisation the reverse-fountain patterns do not emerge. This suggests that in growing tissues strain-induced PIN localisation must be carefully balanced against other mechanisms of PIN localisation to ensure that the correct auxin distributions are established. This hypothesis opens an exciting avenue for further experimental and theoretical investigations of relations between reverse auxin flow and growth processes in plant tissues, especially in plant roots.

For a model considering auxin flux through the apoplast (1), (6)–(8), we compared different mechanisms of transmembrane auxin flux and PIN localisation to examine their influence on the formation of auxin patterns. We found that for flux-based PIN localisation both passage and spot patterns were able to be produced; however, for concentration-based PIN localisation only spots were able to emerge when combined with saturating auxin flux. When PIN was overexpressed compared to AUX1 homogeneous auxin disruptions were observed, agreeing with experimental results of reduced auxin accumulation in cells and pooling in the apoplast (Okada et al. 1991; Yang et al. 2006). Together this suggests some balance between the expressions of PIN and AUX1 is important to facilitate heterogeneous auxin distributions required for stable plant growth. These results also suggest that auxin transport through the apoplast has an effect on the dynamics and distribution of auxin and PIN in plant tissues. Further experimental and theoretical studies of relations between auxin transport through plasmodesmata and through apoplast are important for a better understanding of auxin dynamics and distribution in plant tissues.

We recognise that our model has many components and is more complicated than many other mathematical models which describe the emergence of auxin patterns, e.g. (Feller et al. 2015; Feugier et al. 2005; Hayakawa et al. 2015). We chose to include a good level of biological detail in our model so that the dynamics of all components could be predicted and compared with experimental data. However, our model can be simplified significantly by recognising that the dynamics of TIR1-containing components are faster than other reactions and so can be solved for those variables, reducing Eqs. (1), (2) towhere , , , , and . The reduced model has four fewer variables and four fewer parameters compared to Eqs. (1), (2), and combined with Eqs. (3) demonstrates similar behaviour as the full model (1)–(3) for the same parameter values, as shown in Fig. 14. Since it is highly likely that other hormone signalling networks interact with the auxin signalling pathway to maintain auxin distribution patterns (Bishopp et al. 2011), it would be feasible to combine this simplified model with mathematical models of other signalling pathways. For example, mathematical modelling of interactions between auxin and cytokinin in plant roots has been investigated in Mellor et al. (2017); Muraro et al. (2013) and it would be of interest to examine the dynamics of an extension of the models presented in this work to include cytokinin.

Fig. 14

Reduced model (9) has similar dynamics to the full model (1)–(3) for the same parameter values. New parameters in (9) are calculated from previous parameters as detailed. a Zones of pattern formation and homogeneous distribution for parameters h and resemble those in Fig. 4a. Note . b Model (9) forms passage patterns with similar characteristics as in Fig. 5a with the same parameter values. c For the same parameter values as in Fig. 11c, model (9) generates similar reverse flow patterns at the root tip, where auxin is transported down the root from the central source cells and is then redirect to the outer layers of cells where it is transported back up the root; however, the exact alignment of PIN proteins on membranes is different (Color figure online)

The influence of auxin on root growth is highly complex and still not fully understood (Sengupta and Reddy 2018). The heterogeneous distribution of auxin and its flux through plant tissues are responsible for the development of tissues including root hair (Zhang et al. 2018), vasculature (Marhava et al. 2018), and lateral roots (Li et al. 2018), shoot branching (Ongaro and Leyser 2008), flowering (Cheng et al. 2006), and of course primary root growth (Pelagio-Flores et al. 2016; Wakeel et al. 2018; Wan et al. 2018). Auxin influences these developmental processes through the Aux/IAA signalling pathway modulating transcription of relevant proteins; however, recent results are also exposing control through non-transcriptional effects downstream of the signalling pathway (Fendrych et al. 2018). Auxin transport depends on the dynamics of PIN polarity (Abas et al. 2006), whereas dynamics of PIN depend on auxin-related cellular signalling processes (Vieten et al. 2005). Considering this nonlinear coupling between signalling processes, auxin transport, and PIN dynamics, we hope that our mathematical model and analysis of nonlinear interactions between auxin flux, cellular signalling pathway, PIN dynamics, and growth, as well as hypotheses resulting from our numerical simulation results will contribute to a better understanding of the role of auxin in root development. To our knowledge, the results presented in this paper on oscillatory auxin transport, comparisons between flux-induced and strain-induced PIN localisation, and the formation of the reverse fountain without prescribed PIN polarisation patterns are novel. Our mathematical model for interactions between signalling processes and auxin transport can also be generalised to address possible direct effect of auxin-related signalling processes on the polarisation of PIN (Sauer et al. 2006), once more information about these direct interactions has been found. Further research will also include generalisation of our symplast–apoplast model to include auxin transport through plasmodesmata and to analyse the effect of auxin on transport through plasmodesmata of various signalling molecules (Han et al. 2014).