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(Chen et al. , 2012Robert 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 , root hair formation , 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 mecha-nism 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 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).

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 Here V i represents the volume of cell i, and S m i j 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 V i, j represents the volume of apoplast compartment bordering cell i between cells i and j, and S jk i represents the size of the border between apoplast compartments (i, j) and (i, k) 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.
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

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) 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 equations completed 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, 1 ≤ i ≤ N , and N is the total number of cells. Here, mRNAs are denoted by m i , cytosolic PIN is denoted by p i , auxin is denoted by a i , Aux/IAA is denoted by r i , TIR1 is denoted by s i , auxin-TIR1 complex is denoted by c i , PIN-auxin-TIR1 is denoted by e i , Aux/IAA-auxin-TIR1 is denoted by v i , ARF monomers are denoted by f i , ARF-Aux/IAA complexes are denoted by g i , and ARF 2 dimers are denoted by w i . 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.
Parameter α m is the rate of mRNA production, μ m is the rate of mRNA degradation, φ m is the ratio of ARF-dependent mRNA production to ARF 2 -and double ARFdependent mRNA production, and θ f , θ w , θ g , ψ g , and ψ f 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 α r , whereas β r and γ r are the binding and dissociation rates of Aux/IAA and auxin-TIR1, β g and γ g are the binding and dissociation rates of Aux/IAA and ARF, and μ r is the degradation rate of Aux/IAA from Aux/IAA-auxin-TIR1. By β a and γ a , the binding and dissociation rates of auxin and TIR1 are denoted, whereas β f and γ f are the binding and dissociation rates of two ARF proteins, β p and γ p are the binding and dissociation rates of PIN and auxin-TIR1, and μ p 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 a i , the PIN-mediated flux of auxin between neighbouring cells, and the dynamics of cellular p i and membrane-bound P i j PIN. The index i j denotes the membrane of cell i between two neighbouring cells i and j, e.g. S m i j denotes the size of the portion of the membrane of cell i between cells i and j.
Auxin a i is produced inside the cells with rate α a , degraded with rate μ a , and transported between cells by membrane-bound PIN P i j . Cellular PIN p i is translated from mRNAs with rate α p 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 P i j in that part of the cell membrane (Fig. 3). 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 ODEs where i ∼ j is short notation for j ∈ {k | cell i neighbours cell k} and V i denotes the volume of the cell i. The flux of auxin J i j a between neighbouring cells i and j and the localisation of cytosolic PIN p i from cell i to membrane portion i j facing cell j together with dissociation of membrane-bound PIN P i j back to the cell J i j p are given by Here 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 J i j a , and i∼ j H J i j a ; λ = 1. Parameter φ A denotes the rate of PIN-mediated auxin transport, λ is the maximal rate of PIN localisation to the membrane, δ p 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 1/ 1 + exp{−h J i j a /λ − θ } were chosen such that for fluxes smaller (greater) than a threshold value θ , i.e. J i j a /λ < (>) θ , 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 x i denotes the length of either the horizontal or vertical wall of cell i, χ is the maximum growth rate, and θ x is the threshold for half-maximal auxin-dependent growth rate. Here i j denotes that if x i is the horizontal (vertical) length of cell i then x j 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 y i (Y i j ) denotes the concentration of a chemical in cell i (membrane i j), and F i (F i j ) denotes the reaction terms in the corresponding equations. Since in our model the cell shapes are simplified to be rectangular, S m i j is taken to be the cell length along the appropriate axis, and V i is the product of the length and width of the cell, where dV i dt and d S m i j dt 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 membrane in 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 membranebound 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 . 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 a i , p i , and P i j , includes the dynamics of auxin A i j in apoplast, cellular AUX1 u i , and membrane-bound AUX1 U i j : where J i j a is the flux of auxin between cell i and i j-part of the apoplast, and J i j u denotes the localisation of AUX1 from cell i to membrane portion i j together with dissociation of U i j back to the cell. Parameter α u is the translation rate of AUX1 from mRNA, and μ u 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 i j-part of plasma membrane (part of membrane between cell i and i j-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 J i j A and J jk i for different parts of apoplast, where S w i j denotes the size of the interface between apoplast compartments i j and ji, S jk i denotes the size of the interface between apoplast compartments i j and ik, and V i j denotes the size of apoplast compartment i j. The passive fluxes of auxin through the apoplast and of AUX1 localisation to the membrane are given by: Here φ A is the rate of passive flux of auxin through the apoplast, ω u is the rate of AUX1membrane localisation, and δ u 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 J i j a features saturating auxin transport Jönsson et al. 2006),J i j a is an extension of the flux considered in (3) by including the presence of apoplast, J i j p is a mechanism for PIN localisation, proposed in, e.g. 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, andJ i j p 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 φ a is the rate of passive flux of auxin through the cell membrane, φ p is the rate of PIN-dependent saturating auxin flux, φ u is the rate of AUX1-dependent saturating auxin flux,φ p is the rate of PIN-dependent non-saturating auxin flux, andφ u is the rate of AUX1-dependent non-saturating auxin flux. Parameters κ e f a , κ e f p , κ e f u denote the passive, PIN-dependent, and AUX1-dependent efflux of auxin, respectively, and κ in , κ in p , κ in u denote the passive, PIN-dependent, and AUX1-dependent influx of auxin, respectively. Parameters θ p a , θ u a denote the concentration of auxin for half-maximal transport by PIN and AUX1, respectively. In localisation processes, ω p is the rate of PIN membrane localisation, and δ p is the rate of PIN membrane dissociation. Parameter κ p denotes the proportion of PIN localisation that is auxin-dependent and θ a p is the half-maximal concentration of auxin for auxin-dependent PIN localisation. In J i j p , 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 (s i ) and ARF ( f i ) 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 α a is doubled and in sink cells the rate of auxin degradation μ a 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.
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 J AB a /λ > θ. 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 bound-ary, 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 h − β p 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 α m and β p .
To investigate the likelihood of occurrence of different pattern types when pattern formation can occur, model (1)-(3) are solved for various values of β p 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 β p , 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 (β p 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 β p , 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 β p 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 3 × 3 cells with periodic boundary condition. Specifically, we found that as β p increased the ratio of occurrences of passage patterns to spot patterns reduced dramatically, from as high as ≈ 95% for β p = 1 to ≈ 25% for β p = 100, with fixed value h = 50. 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 ≈ 40% for β p = 1 and ≈ 20% for β p = 100. Since the value of h determines the range of values of β p for which heterogeneous patterns can emerge the probability of certain pattern types emerging for specific values of β p will also vary with h. Although β a and β p 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 β a from 0.5 to 50 for fixed values of β p between 1 and 250 the probability of emergence of specific pattern types was unchanged, with fixed value h = 50. Numerical simulation results for the cases of both passage and spot patterns in auxin distribution are included in Fig. 5 for λ = 0.5. We tested a small set of values of λ, and for higher values of 6 ≥ λ 1.5 the concentration of membrane-bound PIN is increased leading to much higher con-  (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).
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. λ = 0.05, 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 λ = 0.5, as shown in Fig. 6. We tested a small set of values of 0 < λ ≤ 0.5, and observed oscillations in numerical simulations for 0 < λ 0.35 but did not observe oscillations for 0.35 λ ≤ 0.5. 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 (≈ 0.25μ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  Table 2 that generates oscillatory dynamics in a single-cell model for auxin-related signalling pathway, a for λ = 0.05 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 λ = 0.5 the concentration of PIN on cell membranes rises above physically realistic ranges of ≈ 0 − 5 μM and numerical solutions display no oscillatory dynamics (Color figure online) 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, ≈ 0.3 μM for auxin and ≈ 0.1 μ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 λ + ν = 0.5. 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 (ν = 0), are presented  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 λ + ν = 0.5 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 0.4 < λ ≤ 0.5 auxin still flows from source cells to the base of the tissue as well as branching out and back up the outer cells, when 0.15 ≤ λ < 0.4 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 λ < 0.15, 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.
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 fluxinduced 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 κ p , 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 0.5 ≤ κ p ≤ 1, and homogeneous distributions were generated for 0 ≤ κ p < 0.5 (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 β p 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 0 < β p ≤ 100 and found that auxin pooling in apoplast compartments occurred for 0 < β p 50, whereas auxin concentrations were more realistic for 50 β p ≤ 100. 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, α u = 0, or PIN is overexpressed compared to AUX1, α p > α u = 1, (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-fluxrelated 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 distri-bution 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 auxinrelated 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 spotgenerating 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 'reversefountain' 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 straininduced 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) to whereμ r = μ r S tot ,μ p = μ p S tot , θ a = β a /γ a , θ r = β r /(γ r + μ r ), and θ p = β p /(γ p + μ p ). 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. 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 , 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;. 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 straininduced 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).

Table 2
Parameter values for the model of auxin-related signalling pathway (1)    Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.