Self-Organization of Plant Vascular Systems: Claims and Counter-Claims about the Flux-Based Auxin Transport Model

The plant hormone auxin plays a central role in growth and morphogenesis. In shoot apical meristems, auxin flux is polarized through its interplay with PIN proteins. Concentration-based mathematical models of the flux can explain some aspects of phyllotaxis for the L1 surface layer, where auxin accumulation points act as sinks and develop into primordia. The picture differs in the interior of the meristem, where the primordia act as auxin sources, leading to the initiation of the vascular system. Self-organization of the auxin flux involves large numbers of molecules and is difficult to treat by intuitive reasoning alone; mathematical models are therefore vital to understand these phenomena. We consider a leading computational model based on the so-called flux hypothesis. This model has been criticized and extended in various ways. One of the basic counter-arguments is that simulations yield auxin concentrations inside canals that are lower than those seen experimentally. Contrary to what is claimed in the literature, we show that the model can lead to higher concentrations within canals for significant parameter regimes. We then study the model in the usual case where the response function Φ defining the model is quadratic and unbounded, and show that the steady state vascular patterns are formed of loopless directed trees. Moreover, we show that PIN concentrations can diverge in finite time, thus explaining why previous simulation studies introduced cut-offs which force the system to have bounded PIN concentrations. Hence, contrary to previous claims, extreme PIN concentrations are not due to numerical problems but are intrinsic to the model. On the other hand, we show that PIN concentrations remain bounded for bounded Φ, and simulations show that in this case, loops can emerge at steady state.


Introduction
Among its many functions, the plant hormone auxin is known to induce the formation of primordia in the shoot apical meristem, at locations where it accumulates at sufficient levels [1].
As these primordia become plant organs, the distribution of auxin at the microscopic scale of meristems is a fundamental determinant of the final shape of the plant at the macroscopic scale. The molecular mechanisms underlying auxin patterning are not fully known, but it is clearly established that auxin is actively transported between cells by transporter proteins, including the family of PIN proteins. Although it is known that the distribution of these transporters on cell membranes is asymmetric and depends on the distribution of auxin in neighbouring cells, the molecular details of this dependency remains a crucial unknown in the biology of auxin.
PIN proteins have been identified as efflux carriers (i.e., they transport auxin out of the cell to the intercellular space), and are either oriented towards or against the auxin gradient, depending on the developmental stage of the primordium. More precisely, PIN proteins initially orient themselves towards a precise region of the meristem, away from older primordia. Auxin accumulates at the PIN convergence point, and induces the formation of a new primordium; during this phase, PIN are oriented against the auxin gradient. Incipient primordia then act as auxin sinks and induce auxin depletion in surrounding cells. Soon after initiation, PIN orientation switches in adaxial cells (i.e. cells placed on the side towards the axis of a primordium, usually the upper side) from "towards the center" to "towards the outside" (see Fig. 1), whereas PIN in the primordium orient themselves basally towards the inner tissue of the meristem (see Fig. 1 a). These processes lead to canalization, that is, to the development of vascular strands which transport auxin downward from the primordium to inner tissues. As reviewed by [2], PIN polarization reversal is enhanced by the increase of auxin synthesis in incipient primordia soon after their initiation [3,4]. To date, it is not established how PIN is oriented in response to auxin gradients, see e.g. [5] for a recent review.
The gaps in knowledge and the complexity induced by the large numbers of molecules at both intra-and inter-cellular scales make it difficult to treat the process by intuitive reasoning alone; mathematical models are therefore required to understand these phenomena and propose biologically plausible scenarios [6]. Based on experimental results, many authors including [7,8,9,10,11,12] have built computational models relying on the following hypotheses: auxin is synthesized in all cells of the shoot apical meristem at rates that may vary, and is transported from cell to cell passively by diffusion and actively through efflux and influx carriers. The different models can be classified in two main classes, depending on the key hypothesis underlying PIN orientation (see Fig. 2): • The flux-based hypothesis is based on the work of [13], and was first formalized by [14,15] and later by [7,10,12,16]. The assumption is that PIN concentration increases in membranes according to the strength of the auxin flux. This induces a positive feedback loop on the auxin flux, which is then amplified. It was first introduced as an intuitive explanation for the venation process. Models based on this hypothesis have been used to reproduce correctly the orientation of PIN proteins towards and against the auxin gradient [10], leading to the canalization/venation process. However, flux-based models have been criticized because of the lack of evidence for the existence of auxin flux sensors [17].
• The concentration-based hypothesis is based on the experimental result that PIN proteins orient towards incipient primordia, and thus presumably toward auxin maxima. In these models it is assumed that the rate at which PIN proteins accumulate on a membrane depends on the auxin concentration in the neighbouring cell [8,9,11]. How does a cell detect the auxin concentration in neighbouring cells? [9] cites an experimental result of [18] showing the negative feedback of auxin concentration on PIN endocytosis (a process which removes PIN from the membrane to the cytosol). Moreover, [19] proposed the following explanation: auxin accumulation in cells induces growth and expansion, which can be detected by the neighbouring cells due to mechanical forces applied on the cell membrane.
Most published studies have investigated models relying on one or both of the above hypotheses by means of computer simulations, aiming to reproduce realistic auxin distribution patterns in 'virtual meristems'. However, any mathematical implementation of these hypotheses involves a large number of parameters, few of which are experimentally known to date. The exact patterning abilities of each model, sometimes the source of scientific controversy, thus cannot be fully characterized by these simulations: the fact that a particular type of pattern has not been observed in simulations might simply be explained by unexplored parameter regimes.
Until new experimental data more precisely explain the underlying molecular mechanisms and parameter values, the only possible definitive assessment of a model is through rigorous mathematical deduction and the derivation of formal proofs. As all existing models are highdimensional and non-linear, only a few such results have been obtained to date.
In [11], concentration-based models have been shown to produce auxin patterns in a way that is similar to the well-known reaction-diffusion models, through transport-induced bifurcations from a steady state corresponding to a homogeneous distribution of auxin. In particular, these models can produce a variety of spotted and striped patterns, reminiscent of primordia and veins, respectively. The same class of models were investigated more recently in [20], where it was proved that typical patterns took the form of isolated blocks within a background of auxin-depleted cells. Using rigorous numerical bifurcation techniques, [21] were able to show that a typical bifurcation from a homogeneous auxin steady state could lead to either regularly-spaced spikes of auxin, or to stable oscillations. For flux-based models, [16] have shown that different distributions of PIN can lead to a homogeneous auxin distribution, and that these PIN patterns can be characterized in graph-theoretical terms. These authors also give conditions under which homogeneous auxin patterns are stable, and show that stable oscillations of auxin can occur in 1-D rows of cells through flux-based transport.
In this paper, we continue the study of flux-based models. We focus on the equations proposed by [10] as they have been able to reproduce PIN polarization toward incipient primordia (and thus against the auxin gradient), the reversed polarization (along the auxin gradient), phyllotactic patterns and venation. [10] proposed that the two different PIN polarization behaviours are consequences of different cellular responses to auxin. In the cells of the epidermal layer, the response function linking auxin flux to PIN insertion on the membrane is assumed to be linear, whereas the inner tissue response function is assumed quadratic.
However, several simulation studies raised problems for the model [22,23], the canals were predicted to contain reduced auxin concentrations, while experimental studies show that auxin concentrations within canals can be higher than the background concentration, see, e.g., [24] or [25]. The model has thus been extended in various ways. [26] modeled auxin influx carriers. [23] combined both (the concentration-and flux-based hypotheses) in a hybrid model, each cell having both (flux-and concentration-based PIN cycling) depending on the auxin concentration. [12] modified the model of [10] by also considering PIN cycling and the concentration within the cytosol of each cell. They argued that it is then able to predict canals with higher auxin concentration than in the surrounding tissues.
Here we will show that the model of [10] is itself able to generate auxin canals and PIN polarization with and against auxin gradients, and prove that it is in fact able to create canals having higher auxin concentrations than those in the surrounding cells. Mathematical proofs are given that lead to precise parameter regimes that yield such behaviours. Hence, contrary to what is claimed in the literature, the model reproduces patterns of auxin concentrations like those seen in nature. Previous studies dealing with an unbounded response function introduced a cut-off [12] or halted simulations [7] when auxin and PIN concentrations became too high; the authors argued that this behaviour was due to numerical problems. We prove that an intrinsic property of the model is divergence in finite time, and thus that high values of PIN concentration are not due to numerical problems but are a consequence of the computational model itself. We finally study the steady state patterns. These are, at least for quadratic response functions, composed of directed forests of directed trees rooted at sinks, in the absence of cycles. They thus ressemble the plant vascular system. Special emphasis is given to source-and sink-driven systems, for which exact solutions are provided.

The flux model
We suppose that the set of cells i 2 V = {1, Á Á Á, M} is arranged as a graph G ¼ ðV; EÞ. The edge set is composed of pairs of cells (i, j) which are nearest neighbours, denoted as i * j in what follows. The flux model describes the time evolution of the concentrations of auxin molecules a i (t) and of PIN proteins p i (t) within each cell i. We thus need the following variables, for each cell i and for each interface (i ! j), where some of the PIN proteins which are contained in cell i are positioned, facing cell j, see Fig. 3: S ij ¼ surface area of the membrane between cell i and j in m 2 N i :¼ set of neighbours of cell i: All the parameters and variables used in this work are summarized in Table 1.  Table 1. Notations.

G ¼ ðV; EÞ
Graph of node set V and of edge set E.
Oriented sub-graph of node set V* and of oriented edge set E*. Auxin flux through the membrane between cell i and j due the diffusion.
Auxin flux through the membrane between cell i and j due the PIN active transport.
γ A Auxin active transport rate.

D
Normalized auxin diffusion rate.
T Normalized active transport rate.  Following [10], we will consider the flux associated with diffusion and active transport.
¼ diffusive auxin flux through the membrane between cell i and j in mol m À2 s À1 ; ¼ auxin flux associated with active transport; is a permeability constant reflecting the capability of auxin to cross the membrane. The flux due to active transport across a membrane separating cells i and j is modelled as characterizes the transport efficiency of the PIN pumps.
[10] proposed the following set of ordinary differential equations where the coefficients a i a (mol m −3 s −1 ) and b i a (s −1 ) are auxin production and degradation factors of cell i. ρ 0 (mol m −2 s −1 ) is the basal PIN insertion rate into the membrane, μ (s −1 ) is the basal removal rate from the membrane. The first equation (1a) gives the rate of variation of auxin concentration in cell i as a function of protein production and degradation, and of the flux J i!j through the membrane ij. The second equation reflects [13] original concept that canalization is induced by a positive feedback between flux and transport. This canalization hypothesis was then formalized in [14,15]. According to the second equation (1b), the variation of the concentration p ij of PIN proteins transporting auxin to cell j is a consequence of 1) insertion in the membrane induced by the flux and 2) basal insertion and removal of PIN protein from the membrane. The intensity of PIN insertion into the membrane i is modelled using a non-linear response function h. The positive feedback mechanism is modelled by assuming that h is increasing and non-negative. It is furthermore assumed that no PIN proteins are removed from the membrane when the number of incoming auxin molecules is larger than the number of outgoing auxin molecules; in mathematical terms, this means that Typical examples of such functions are where κ is a positive parameter in mol m −2 s −1 and J ref is an arbitrary reference flux in mol m −2 s −1 . [7] also consider functions of the form It is important to note that the properties of the steady states associated with the o.d.e. (1) strongly depend in general on the nature of h [7,27]. We will provide more precise results on this topic in what follows. In [7,12,16] a modified model is also considered where PIN transport (between the membrane and the cytosol) is distinguished from PIN degradation/creation in each cell i. We follow the version proposed in [16] which includes the time evolution of the PIN concentrations P i (t) in the cytosol of each cell i, and takes the form: where μ is the removal rate of PIN from membrane. For this new model, the global PIN concentration inside each cell is preserved when α p = β p = 0, that is, Let a = (a i ) i 2 V , P = (P i ) i 2 V and p = (p ij ) (i, j) 2 E be the vector containing auxin and PIN concentrations in the cells and on the membranes.
[10] studied the system of differential equation (1) numerically, letting vary parameter values including the auxin synthesis and degradation rates in each cell i. They were able to reproduce PIN polarization with and against auxin gradients. The system with h(x) = x 2 can also reproduce the canalization phenomenon where canals formed by PIN polarization are surrounded by cells without transporter. With a linear feedback function h(x) = x, [10] reproduced phyllotactic patterns and an inhibiting field. [12] showed using simulations that [10]'s model is able to generate phyllotactic patterns and canalization without changing the feedback function h, simply by altering the auxin synthesis/ degradation rates. Indeed, if incipient primordia acts as sinks by degrading auxin at higher rates than other cells, they attract and accumulate auxin molecules. The primordia then become sources as they grow, inducing the initiation of the inner plant vascular system.
From the fact that h(x) = 0 for x 0, the equilibrium (a Ã , P Ã , p Ã ) requires that for every pair of nearest neighbour cells i * j, either p Ã ij ¼ 0 or p Ã ji ¼ 0. This property allows an orientation of the graph G to be associated with each steady state, where an edge (i, j) is oriented as i ! j if and only if p Ã ij > 0. Note that with this procedure p ij = p ji = 0 implies that the corresponding edge is actually removed in the oriented graph; the procedure thus leads to an oriented sub-graph denoted by G in what follows. This oriented graph was defined in [16], where it was used to characterize further properties of steady states with homogeneous auxin, i.e. a i = a Ã 8i by the condition that G must be well-balanced.

Model simplification
Consider a graph G ¼ ðV; EÞ with regular cells, i.e., the surface areas between membranes S ij and the cell volumes V i are all equal. Let W S ij V i ; 8(i, j) 2 E. Set: We will focus on a simplified model which is obtained by assuming that for a small parameter % 0. This regime assumes that (passive) diffusion is negligible compared to active transport. Without loss of generality, we assume that T = 1. In this regime, where J i ! j = (a i p ij − a j p ji ). We show moreover in the Supporting Text that if a solution xðtÞ ¼ ðaðtÞ; pðtÞÞ of (St0) converges towards an equilibrium x Ã = (a Ã , p Ã ), then, for small enough , the solution (a (t), p (t), P (t)) of (3) having the same initial condition remains in a small neighbourhood of xðtÞ for all future time.

Divergence in finite time
Most of the papers dealing with (1) or (St0) with unbounded F introduce a cut-off or stopped simulations when the values of PIN concentrations p ij became too large (see, e.g., [7,12]). These papers argue that this divergence is due to numerical problems. We show that these extreme values are not due to numerical problems, but that the p ij (t) can diverge toward 1 in finite time; this is hence an intrinsic property of the computational model. We prove in S1 Text that the slow equation (St0) has an unique solution, which is defined for all t ! 0 when F is a bounded function if a i a ; b i a ; m; l; a p =b p > 0, 8i 2 V. When F is not bounded, either the solution is defined for all t ! 0, or only for t in a bounded interval [0, δ); in the latter case, the solution diverges toward 1 as t ! δ. An example where PIN concentrations on membrane diverge in a finite time is provided in the supporting information.

Steady state topologies
As discussed above, every steady state leads to an oriented sub-graph G ¼ ðV; E Ã Þ of the basic graph G ¼ ðV; EÞ. For every pair of nearest neighbours i * j of E, E Ã either does not contain an oriented edge linking these two cells, or contains at most one of the two possible directed edges (i ! j) and (j ! i). Hence, Let I Ã be the set of all isolated cells of G (i.e.those cells i 0 for which ]{k i 0 } = 0 and ]{k ! i 0 } = 0). For clarity, denote by G Ã the sub-graph of G where the isolated cells have been removed, i.e., G Ã is the directed graph of edge set E Ã , and of node set Let G Ã be an oriented sub-graph of G. The following notions will be useful in what follows, see Fig. 4

a):
• The out-degree (resp. in-degree) of a cell i in G Ã is the number of cells j such that i ! j (resp. j ! i). The degree is the sum of the out-degree and the in-degree.
• A sink of G Ã is a cell with out-degree 0 and in-degree > 0.
• A source of G Ã is a cell with in-degree 0 and out-degree > 0.

Stable steady state topologies
This section focuses on the stability of the system (St0) with a quadratic F(x) = x 2 . Let E 0 ¼ ða 1 a =b 1 a ; Á Á Á ; a M a =b M a ; 0; Á Á Á 0Þ be the configuration with every p ij equal to zero, i.e., the associated directed graph is composed of isolated nodes only. For our choice of the function F, this configuration is always asymptotically stable, see [16]. From the latter, we also know that steady states with homogeneous auxin and non-zero PIN exist for well-balanced orientations G Ã , but their stability is unknown in general. We now prove that for a quadratic F, the stability of arbitrary steady states (i.e., where auxin is not necessarily homogeneous) can in fact be characterized using some simple necessary conditions on the sub-graph G Ã .

Instability criterion
Theorem 1 Consider the system (St0) with F(x) = x 2 . Assume that a i a > 0, b i a ! 0, 8i 2 V, and that μ, c > 0. Let (a Ã , p Ã ) 6 ¼ E 0 be a bounded equilibrium of the system (St0) of associated oriented sub-graph G Ã . Then, • If G Ã contains no sink cells, then (a Ã , p Ã ) is unstable.
• If (a Ã , p Ã ) is stable, then G Ã is an oriented sub-graph of G composed of trees directed from leaves to roots such that every cell i 2 V has at out-degree 1.
We hence observe that every stable equilibrium is obtained by considering an oriented graph G Ã composed of directed trees pointing to sinks, see Fig. 4 b). The stable patterns do not contain loops, as it has been observed through simulations in most of the papers dealing with the flux model. Simulations show that loops can obtained by using a bounded function F.
To obtain further confirmation of Theorem 1, we performed numerical simulations on a regular grid of cells. In most simulations with homogeneous auxin production rates a i a ¼ a a 8i, the simulations ended on the steady-state E 0 , where all cells are isolated. We thus considered a grid of 8 × 11 cells within which two groups of three cells have a higher production rate, see Fig. 5 a) and b). We observe that simulation seems to converge to a state composed of rooted trees as predicted by Theorem 1. The S1 Text provides also precise formulas for computing locally asymptotically stable configurations for quadratic response functions.
Note that the result above applies at a limit where the diffusion constant vanishes D = 0, and for regular tissues where all cells have the same volume and surface area. Neither of these two conditions is likely to be strictly satisfied in real systems, but they greatly simplify mathematical analysis. Also, we expect Theorem 1 to remain true for systems where the two conditions are nearly satisfied, i.e. where D is non-zero but small, and where cells do not vary too much in size. This intuition was confirmed by numerical simulations, as seen in Fig. 5 c) where the previously observed steady-state is qualitatively unchanged by slightly relaxing the two conditions. We next focus on source and sink driven systems, which play a fundamental role in plant growth, see, e.g., [12].

Sink-driven system
As explained in the Introduction, experiments show that the primordia of the L1 layer act first as auxin sinks, inducing auxin depletion in surrounding cells. We study the stable patterns associated with this phase within the flux-based modelling framework, in the extreme case where there is a single primordium, which is located at some site i 0 such that b i a % 0; 8i 6 ¼ i 0 and b i 0 a > 0: i 0 is evacuating auxin at positive rate while the other cells (of the L1 layer) have low auxin degradation rates. Each cell i produces auxin at rate a i a ! 0. We will in this way obtain exact solutions when b i a ¼ 0, i 6 ¼ i 0 , that permit to determine if really, within the modelling framework given by (St0), auxin is depleted in surrounding cells. When i 6 ¼ i 0 is such that b i a ¼ 0 and a Ã i < þ1, then necessarily i = 2 I Ã . We will show in the S1 Text that i 0 must be the unique sink cell. Hence, the oriented forest G Ã reduces to a directed spanning tree rooted at i 0 . Let G Ã i be the directed rooted sub-tree of G Ã that points to i, of node set V Ã i (see Fig. 6), and let a a ðiÞ ¼ X be the global auxin production rate associated with the sub-tree G Ã i . We show in the S1 Text that the steady state auxin concentrations are given by the exact formulas where we c = λα p /(μβ p ). Starting from a source node of the rooted tree G Ã , the sums α a (i) increase along the unique path to i 0 : one deduces then that the auxin concentration decreases along the paths as long as i 6 ¼ i 0 , so that the paths are directed against the auxin gradients (see Figs. 6 and 7). This also implies the existence of auxin depleted zones in the neighbourhood of the primordium i 0 .

Source-driven system
As stated in the Introduction, experiments show that the primordia of the L1 layer later act as auxin sources from where the internal vascular system initiates. We study here the vascular patterns predicted by the flux-based model when there is an isolated primordium. This model has been criticized because of a perceived inability to produce high auxin concentrations inside initiating veins, see, e.g, [23] and the references therein [9,28]. Contrary to this claim, we will prove that the model can correctly predict auxin concentrations. In the following, we assume that there is a primordium located at i 0 , with auxin production rate higher than that of the other cells. Mathematically, we assume here that and that As G Ã is composed of oriented trees pointing to roots, it contains at least one source, which must be i 0 . Furthermore, we prove in the S1 Text that the fact that G Ã is a forest composed of directed trees, with a distinguished cell i 0 , imply that the vein G Ã must be a linear chain i 0 À!i 1 À!i 2 À! Á Á Á À!i n ; that is, must be a directed line with no vascular strands ending in a sink i n .  We can distinguish two main parameter regimes, leading to different kinds of flux. In the first case, the auxin concentrations satisfy ðAÞ : a b a i 0 a i 1 Á Á Á a i n ; (see Figs. 6, 8 and 9), so that the auxin flux inside the linear vein goes along the auxin gradient.
In the second case it flows against the gradient (see Figs. 8 and 9). We provide precise mathematical formulas defining the parameter regimes leading to (A) and (B) in the S1 Text.
Simulations of these two kinds of source-driven system are provided in Fig. 10. We consider a line of L = 20 cells (without periodic boundary conditions) with i 0 = 1 and take To recover the two regimes (see Fig. 8), we take a) a 1 a ¼ 0:9 leading to PIN polarization with the auxin gradient and a higher auxin level in the canal than in the surroundings Thus, contrary to what has been claimed in previous studies, see, e.g., [23] and the references therein, the flux-based model does not necessarily lead to auxin concentrations within canals that are lower than the background concentration α/β associated with isolated cells. Instead, the resulting configuration depends on the source-strength of i 0 (a i 0 a ) relatively to the other cells (i.e α). One also observes the following counter-intuitive and surprising result: the flux increases along the vein when the production rate α 0 of the source is low, whereas the flux decreases for highly productive sources with large α 0 .

Discussion
We have studied various aspects of the flux-based auxin transport computational model. The model has been criticized because numerical simulations indicate that auxin concentrations within initiating veins are lower than those seen experimentally. Various extensions have been proposed to overcome this problem. We have proved that contrary to what is claimed in the literature, the flux-based model is able to correctly reproduce auxin concentrations within veins. A mathematical analysis permits us to give in a precise way the related parameter regimes. Our analyses show that for quadratic and unbounded response functions F, the steady state vascular patterns are formed of directed trees and thus do not contain loops. Moreover, we have proven that PIN concentrations can diverge in finite time, explaining in this way why previous simulation studies were forced to introduce cut-offs to ensure that PIN concentrations remained bounded. On the other hand, we showed that PIN concentrations remain bounded for bounded F; simulations suggest that in this case, loops can emerge at steady state. Thus, according to this model, the self-organization of the auxin flux leads to different patterns as a function of the boundedness of F. Plants having a tree-like vascular system should therefore be modelled using unbounded F, at the risk of recovering unrealistically high values of PIN  concentration. Finally, our analysis of source-and sink-driven systems reveals new and surprising results: for sink-driven systems, the auxin flux is directed against auxin gradients, leading to an auxin depleted zone in the neighbourhood of the primordium, while for source-driven systems, a linear vein emerges which can have high auxin concentration and be directed with the auxin gradient. Surprisingly, the latter situation occurs when the auxin production rate of the source is relatively low. This model is thus able to reproduce phylotactic pattern and canalization by varying auxin synthesis and degradation rate in each cell. The remaining question is why a cell becomes a sink or source.