Toward a 3D model of phyllotaxis based on a biochemically plausible auxin-transport mechanism

Polar auxin transport lies at the core of many self-organizing phenomena sustaining continuous plant organogenesis. In angiosperms, the shoot apical meristem is a potentially unique system in which the two main modes of auxin-driven patterning—convergence and canalization—co-occur in a coordinated manner and in a fully three-dimensional geometry. In the epidermal layer, convergence points form, from which auxin is canalized towards inner tissue. Each of these two patterning processes has been extensively investigated separately, but the integration of both in the shoot apical meristem remains poorly understood. We present here a first attempt of a three-dimensional model of auxin-driven patterning during phyllotaxis. We base our simulations on a biochemically plausible mechanism of auxin transport proposed by Cieslak et al. (2015) which generates both convergence and canalization patterns. We are able to reproduce most of the dynamics of PIN1 polarization in the meristem, and we explore how the epidermal and inner cell layers act in concert during phyllotaxis. In addition, we discuss the mechanism by which initiating veins connect to the already existing vascular system.


Author summary
The regularity of leaf arrangement around stems has long puzzled scientists. The key role played by the plant hormone auxin is now well established. On the surface of the tissue responsible for leaf formation, auxin accumulates at several points, from which new leaves eventually emerge. Auxin also guides the progression of new veins from the nascent leaves to the vascular system of the plant. Models of auxin transport have been developed to explain either auxin accumulation or auxin-driven venation. We propose the first threedimensional model embracing both phenomena using a unifying mechanism of auxin transport. This integrative approach allows an assessment of our present knowledge on how auxin contributes to the early development of leaves. Our model reproduces many observations of auxin dynamics. It highlights how the inner and epidermal tissues act together to position new leaves. We also show that an additional, yet unknown,

Introduction
In plants, most developmental processes are driven by the spatiotemporal distribution of the growth regulator auxin. The versatility of the morphogenetic role played by auxin relies on its self-regulated polar transport, in which auxin transport feedbacks on auxin efflux carriers, the PIN proteins [1]. This process can lead to various distribution patterns, depending on the specific geometry of the organ considered. Although auxin transport and patterning occur fundamentally in three-dimensional tissues, they are usually explored in one-or two-dimensional models. These restrictions are commonly justified by geometrical considerations. For instance, models of phyllotactic patterning in the shoot apical meristem (SAM) of Arabidopsis assume that the formation of convergence points of PIN1 polarization at primordia takes place in the single epidermal cell layer L1 [2][3][4][5], based on observations by Reinhardt et al. [6]. However, it is known that inner tissues are essential for positioning primordia [7][8][9]. Bayer et al. [10] proposed an auxin transport model integrating both the formation of convergence points in the L1 and the patterning of vascular strands in the subepidermal layers. Their model is implemented on a 2D cellular template representing a longitudinal section through the meristem. Thus, primordia positioning and midvein development are simulated in, respectively, one and two dimensions. These restrictions in dimensionality limit the range of potential behaviors displayed by the model, and therefore its capacity to integrate and assess current knowledge on phyllotaxis. Although fully three-dimensional models would be a significant step forward, they still present technical challenges [11]. In a notable effort in this direction, Gruel et al. [12] modeled gene expression dynamics and cell-to-cell diffusion of signals on a 3D cellular template to explore positioning and maintenance of the stem cell niche in the meristem. Modeling auxin polar transport is arguably more challenging, since it requires in addition a description of PIN exo-and endocytosis between cytoplasm and membranes, and, ideally, auxin diffusion in the apoplast. Several technical difficulties pertaining to the modeling of 3D plant tissues could be alleviated by using the topological notion of cell complexes. It provides a framework for locally describing how the components of a plant tissue (cell interiors, membrane elements, cell walls) are connected to each other. Based on such a description, efficient algorithms can be developed for simulating changes in the structure of the tissue or flows within the tissue [13]. An implementation of this paradigm has been utilized by Yoshida et al. [14] to investigate the control of division orientation in early Arabidopsis embryogenesis. Cell complexes are likely to turn out even more useful when fluxes between tissue components are involved.
Another defying peculiarity of the SAM is the co-occurrence of convergence and canalization at midvein initiation. In the L1, PIN1 polarize towards convergence points, where auxin accumulates. From these points, strands of cells with high PIN1 expression extend into the subepidermal layers, with PINs displaying a canalization pattern. This co-occurrence suggests that these two modes of auxin transport regulation do not rely on completely different mechanisms, but instead relate somehow to each other.
Convergence and canalization are commonly conceptualized using two distinct polarization models, respectively referred as "up-the-gradient" and "with-the-flux". There has been interest to go beyond this dichotomy. Some models attempted to explain phyllotaxis with either purely up-the-gradient [15] or with-the-flux [5] polarization, but they lack biological plausibility or contradict experimental data [16]. Following another approach, Bayer et al. [10] put forward the concept of "dual polarization", whereby both polarization mechanisms operate concurrently, with a continuous transition from up-the-gradient to with-the-flux polarization depending on local auxin concentration. Although the dual-polarization model displays good agreement with observations, it remains difficult to explain how two qualitatively different mechanisms of PIN1 localization can coexist within the same cells. Shifting attention from Arabidopsis to Brachypodium, O'Connor et al. [17] took advantage of the existence of the PIN1 duplicate sister-of-PIN1 (SoPIN1) proteins in grasses by attributing different roles to the two proteins. However, no such duplication exists in Arabidopsis, where both convergence and canalization are apparently performed by PIN1 protein alone.
Moreover, there are fundamental problems with most polarization models [16]. In withthe-flux polarization, it seems unrealistic that cells directly sense net fluxes of auxin through their membranes. Regarding up-the-gradient polarization, it is not clear how a cell could react to the auxin concentration in their neighbors. Several alternative mechanisms have been proposed. Wabnik et al. [18] assumed that auxin gradients in the apoplast are informative enough to drive polarization. Although their model is capable of transitioning between up-the-gradient and with-the-flux regimes, and might thus also account for phyllotaxis, it is yet to be determined whether significant auxin gradients can form in the very narrow spaces between meristematic cells. Abley et al. [19] proposed a model for tissue polarity based on cell-cell coupling through a diffusive mediating molecule. They used it to reproduce several polarization behaviors [20], but phyllotactic patterns and midveins initiation seem out of reach.
The central role of auxin transport as a polarity patterning factor has been questioned in the light of several findings (see a review in [21]). For instance, it has been shown that PIN1 polarities can be oriented by the activity of a transcription factor [9], or that PIN1 polarities are not disrupted when polar auxin transport is blocked [22,23]. As an alternative mechanism, it has been proposed that PIN polarities are influenced by mechanical perturbations [24,22,25]. A model based on this approach can generate a whorled pattern of auxin maxima [22], but it is not yet known whether it could reproduce the range of polarizations observed in phyllotaxis. In the present work, we adopt the point of view that the auxin/PIN1 system is the central regulator, with feedbacks from downstream effectors (transcription factors, cytokinins, mechanics).
In all up-the-gradient models, PIN polarization is assumed to be fast compared to the production and turnover of PINs, as well as to changes in cellular auxin concentration, so that PIN concentrations at membranes are set to their steady-state values [11]. Since the cycling rates of PINs are completely unknown [16], this assumption may not be valid. If there is indeed significant latency between changes in auxin concentration and PIN polarization, upthe-gradient models could fail to capture some aspects of convergence point formation.
In an attempt to explain how cells can measure the direction and magnitude of auxin fluxes, Coen et al. [26] hypothesized the existence of "tally molecules" produced or consumed at the membrane when auxin enters/exits the cell. The concentration of these molecules at a membrane would act as a proxy for the magnitude of auxin influx and efflux through the membrane. Cieslak et al. [27] presented several biochemically plausible implementations of this concept, assuming that tally molecules modulate PIN allocation to membrane, thus giving rise to a feedback of auxin fluxes on PIN localization. When a local increase in auxin influx decreases the abundance of PIN proteins in the corresponding part of the membrane, canalization patterns emerge. Conversely, when auxin influx locally increases PIN allocation, convergence points form. The transition between these two regimes is controlled by variations in intracellular auxin concentration: An excess of auxin, relative to a threshold concentration, causes a sharp switch of the rate of a reaction occurring inside the membrane. This, in turn, provokes the transition from up-the-gradient to with-the-flux polarization. In addition, thanks to its detailed description of the underlying biochemical reactions, the model does not set PIN concentrations at membranes to their steady-state values, and can thus capture potential transient states during fast polarization events. However, it has been tested only on two-dimensional square grids.
In order to explore three-dimensional aspects of auxin patterning during phyllotaxis, we adapt Cieslak's model to arbitrary irregular 3D tissue geometries, and implement it using the cell complex paradigm. We begin by running the model on a single layer of cells to reproduce convergence point formation. We then investigate how this process is affected by inner cell layers. Finally, we determine under which conditions midveins initiated at convergent points in the epidermal layer can progress in the inner tissue and eventually connect to the already existing vascular system. We show that Cieslak's mechanism produces dual polarization in three dimensions, and discuss the conditions under which it can happen. Furthermore, in line with previous studies, we highlight the necessity of an additional mechanism for guiding developing veins towards preexisting vasculature.

Model of biochemical reactions and transport
We adapted to 3D tissues a network of reactions and transport initially formulated by Cieslak and colleagues for two-dimensional grids of square cells (see figure 9 in [27]). The resulting network is represented as a Petri net in Fig 1A (the system of differential equations derived from this net is detailed in S1 Appendix). It is based on the following hypotheses: 1. Auxin efflux is a two-step process. In the first step, cellular auxin (A c ) forms a complex with membrane-bound PIN proteins (PIN m ): where T out1 is the rate of APIN formation. The APIN complex is bound to the membrane. In the second step, APIN spontaneous dissociation releases auxin in the apoplast (A a ): where T out2 is the rate of APIN spontaneous dissociation. As a result of this two-step process, the local concentration of APIN is proportional to the local magnitude of auxin efflux. In this sense, the APIN complex plays the role of a tally molecule measuring auxin efflux. The proportionality coefficient (or conversion factor) between APIN concentration and auxin efflux is 1/T out2 .
2. Auxin efflux locally increases exocytosis of cellular PIN (PIN c ) through the tally molecule APIN: where σ apin is the rate of PIN exocytosis by APIN. Two APIN molecules are needed to transport one PIN from the cytoplasm to the membrane. This ensures quadratic feedback of auxin flux on PIN exocytosis [28][29][30].
3. Auxin influx also occurs in two steps. First, apoplastic auxin forms a complex with AUX/ LAX efflux carriers: where T in1 is the rate of AAUX formation. The AAUX complex is bound to the membrane. In the second step, AAUX spontaneously dissociates: where T in2 is the rate of AAUX spontaneous dissociation. After dissociation, auxin is released in the cytoplasm. Paralleling APIN, the local concentration of AAUX is proportional to the local magnitude of auxin influx, and thus plays the role of a tally molecule measuring auxin influx. The proportionality coefficient (or conversion factor) between AAUX concentration and auxin influx is 1/T in2 . The symbol � represents reactants or products outside the modeled system (auxin precursors, for instance). Reactions are pictured as rectangles labeled with their rate constants (see main text for an explanation of each term). Arrow orientation indicates whether a species is a reactant or a product in a given reaction. A double-headed arrow indicates a catalytic activity. The two double-headed red arrows represent the impact of the tally molecules APIN and AAUX on PIN exocytosis. They are involved in exocytosis reactions as catalysts with stoichiometric coefficients equal to 2. Other stoichiometric coefficients are equal to 1. An exception is auxin-dependent PIN biosynthesis (rectangle with an S-shaped curve), which follows a more complex kinetics with two parameters (see main text). (B) Graphical convention for simulation results (simplified representation in 2D). Cell interiors are colored in green, with a level of opacity proportional to auxin concentration; for instance the bottom left cell has a much higher concentration than the top cell. Red rectangles on membrane elements picture PIN proteins bound to these elements, with a width proportional to PIN concentration. Arrows make local polarization more apparent; they are deduced from PIN localization. (C) Lateral diffusion of PIN proteins on a cell membrane. Membrane-bound PINs are not fixed at some position on the membrane but freely diffuse on it. This idea is pictured on a 2D section of a cell, with membrane-bound PINs represented as red ellipses. On average, PINs move from regions with higher concentrations to regions with lower concentrations. Since we discretized each cell membrane into a set of polygonal faces (referred to as membrane elements), PIN diffusion was modeled as fluxes between adjacent faces of a given cell membrane, using a discrete Fick's law. https://doi.org/10.1371/journal.pcbi.1006896.g001

Auxin influx locally increases PIN exocytosis through the tally molecule AAUX:
where σ aaux is the rate of PIN exocytosis by APIN. Again, the feedback of auxin flux on PIN exocytosis is quadratic.

AAUX catalyzes dissociation of APIN:
After dissociation of APIN, auxin is released back into the cytosol. ν apin is the rate of APIN dissociation by AAUX and its value depends on auxin concentration in the cell. For low auxin concentration, ν apin is close to zero, so that auxin influx acts only positively on auxin efflux. This leads to convergence point formation. Conversely, for high auxin concentration, ν apin has a high value, so that AAUX rapidly breaks up APIN. Therefore, auxin influx hinders auxin efflux, which leads to canalization. The transition between the two extremal values of ν apin follows a sigmoid function (see S1 Appendix).
6. Besides exocytosis mediated by tally molecules, there are constitutive exo-and endocytosis of PIN proteins, with respective rates σ p and μ p .
7. In Cieslak's original scheme, a fixed pool of PIN proteins is assumed in every cell. This assumption does not held for the full SAM, in which most inner cells contain much less PIN proteins than epidermal cells. Furthermore, there is no evidence for a fixed cell pool of PINs, while auxin-dependent PIN biosynthesis is strongly supported by experiments [31]. Therefore, we describe changes in the total PIN concentration of a cell by the equation where [A c ] is the concentration of auxin in the cell, [PIN c ] is the concentration of cytoplasmic (non-exocytosed) PIN in the cell, ρ p captures the up-regulation of PIN biosynthesis by auxin, κ p controls saturation of PIN biosynthesis, and μ p � is the rate of PIN decay. This equation is inspired by the model of phyllotaxis of Smith et al. [4]. We want to stress, however, that the hypothesis of auxin-dependent PIN biosynthesis is not indispensable for canalization and convergence to occur. The same behaviors can be displayed with fixed pools of PINs.
8. In each cell, auxin is synthesized and degraded at fixed rates (respectively σ a and μ a ). In the apoplast, auxin diffuses freely.
9. To account for the fact that tissues surrounding the SAM act more as auxin sinks than as sources, cells at the border of the template do not produce auxin. The five most central cells are also non auxin-producers, to avoid convergence point formation in this region. This is in agreement with modeling and experimental data [2,32].

Parameter determination and correlations
Since the reaction network used in the model is putative, no experimental values are available for most of the parameters. PIN exo-and endocytosis are established processes, but the associated rates are not known. For parameters already present in Cieslak's model, we took the same values, except for auxin production (σ a ), auxin turnover (μ a ), and the threshold auxin concentration (a th ) for the shift in the value of ν apin . In our 3D model, we had to decrease auxin production, increase auxin turnover, and raise the threshold auxin concentration. Because simulations on 3D templates take very long time, it was not possible to run thousands of simulations to evaluate parameter uncertainties and make a sensitivity analysis. However, we did such an analysis on a 2D square grid to know whether there were correlations between parameters before moving to 3D templates. It turned out that strong correlations were numerous. For instance, the constitutive rate of PIN endocytosis (μ p ) anticorrelates with the rate of APIN spontaneous dissociation (T out2 ), while the latter correlates with the rate of PIN exocytosis by APIN (σ aaux ). Or the auxin threshold (a th ) for the shift in the value of ν apin anticorrelates with auxin decay (see S1 Appendix for a full description of correlations). However, as explained above, a few parameters had different values in the 3D simulations. Thus, 2D analysis does not fully replace 3D analysis. To complement the 2D analysis, we assessed the sensitivity of the 3D model to two critical parameters: a th and the transition steepness (see S1 Appendix).

3D model of tissue
We constructed a 3D model of meristematic tissue. Each tissue template was built in two steps. First, we tessellated a given volume with truncated octohedra. The truncated octohedron has the advantage of filling space, while having a more complex shape than other convex space-filling polyhedra such as cube and prisms, which makes it more realistic as a meristematic cell. In a second step, the vertexes of our mesh were moved by random amounts to introduce irregularities in the template. Each face of the polyhedral cells represents a discrete element of cell membrane. When two cells share a face, this face defines two discrete elements of membrane (one for each cell) and a discrete element of apoplast (shared by both cells).

Implementation and visualization of simulations
The model was implemented in C++ using the VVe modeling environment, an extension of the VV system [33]. The tissue is represented as a 3D cell complex. The representation and associated topological operations have been adapted from the work of Brisson [34,14]. The simulations are visualized using the following graphical convention (see Fig 1B): • The opacity of a cell is set proportional to its cellular auxin concentration. A cell devoid of auxin is transparent.
• Cells with a high auxin biosynthesis rate (most of L1 cells) are colored in green.
• Cells with a low auxin biosynthesis rate (inner cells plus a few L1 cells) are colored in white.
• Cells with a high auxin turnover rate (i.e. sink cells) are colored in purple and are fully opaque.
• Membrane-bound PINs on a membrane element are represented as a red plate whose thickness is proportional to PIN concentration.

Cieslak's mechanism can produce auxin convergence points in a single layer of irregular 3D cells
We first run our model on a tissue template consisting of 261 cells arranged in a single planar layer with a roughly circular shape. Initial auxin concentration was set to zero in all cells.
There was initially no PIN in the cells and very little PIN on membranes. Note that the resulting patterns were independent from these initial conditions, as long as auxin initial concentration was below the auxin threshold. In the first part of the simulation, auxin concentration increased steadily until instabilities occurred and PINs started to polarize. Then, convergence points formed rapidly (Fig 2A and S1 Movie). A visible unrealistic feature in the PIN convergence pattern obtained from this simulation was that some membrane elements with a small area were favored over elements with a larger area in terms of PIN allocation. PINs tended to accumulate on them. This especially happened when two cells with very different auxin concentrations shared only a very small element of cell wall. Then, an auxin flux tended to establish, which promoted PIN allocation to the corresponding membrane element(s). Since the membrane elements of the two adjacent cells were very small, the flux kept a high magnitude, so equilibrium was reached very slowly. In the meantime, membrane-bound PIN concentration reached unrealistically high values. When the membrane elements involved were extremely small, this could even lead to numerical instabilities.
That aberration has various causes. First, discretizing cell membranes and walls, as it is commonly done in tissue modeling, introduces some distortions in the way fluxes are modeled. This does not represent such a problem in regular (square or hexagonal) cell grids, in which all discrete elements have the same size. Three-dimensionality makes the problem more acute because it spontaneously generates a broader range of areas. Second, it is unrealistic not to set an upper limit to the concentration of membrane-bound PINs since PIN proteins occupy some space on a membrane and thus can not reach too high densities due to steric limitations. A solution could have been to assume that PIN allocation follows some limiting kinetics, for instance Michaelis-Menten. But as this sounded a bit ad hoc, we chose instead to assume that membrane-bound PINs diffuse between adjacent membrane elements of the cell (Fig 1C; see mathematical details in S1 Appendix). Indeed, although some membrane proteins are tethered to the cell wall, it has been shown that PINs are actually diffusing in the membrane [35]. Through diffusion, PINs would move away from higher-density membrane elements to neighboring lower-density elements, and thus could not over-accumulate on a single element.
We ran another simulation on the same tissue template, but this time assuming significant lateral diffusion of membrane-bound PIN proteins. We obtained a similar pattern of convergence points, except for the tendency of PINs to accumulate on small membrane elements, which was no longer observed, as expected ( Fig 2B and S2 Movie). The PIN lateral diffusion hypothesis did not only eliminate over-accumulation of PINs on small parts of a membrane (since, as expected, PINs diffused from regions with higher concentrations to adjacent regions with lower concentrations), but it also mitigated the distortions caused by artificial membrane discretization, since the separation between membrane elements was made less strict. Furthermore, it made the model of PIN polarization more dynamic: membrane-bound PINs could then also be reallocated through lateral diffusion, not only through endo-and exocytosis.

Inner cell layers actively counters leak of auxin from the L1
We considered a tissue template with four layers of cells to address the question of how auxin is kept confined in the L1 until venation begins. We ran several simulations on this template, setting a high auxin biosynthesis rate in L1 cells (except from cells at the border and center, as previously) and a low auxin biosynthesis rate in inner cells (L2, L3, and L4).
We first assumed that only epidermal cells can synthesize PINs, while inner cells are devoid of PINs. Our simulations could not reproduce convergence points in the epidermal layer due to significant leaks of auxin into inner layers. This prevented auxin concentration in the L1 to reach the critical value at which convergence points form. To investigate whether PIN proteins in inner layers could influence auxin patterning in the epidermal layer, we amended our assumption by assigning auxin-dependent PIN biosynthesis to all cells. We could then recover the epidermal patterning (Fig 3 and S3 Movie). This was only possible if the rate of auxindependent PIN biosynthesis was high enough, so that PINs were quickly produced in the L2 as soon as auxin was entering L2 cells from the L1. The newly produced PINs polarized towards the L1, obeying up-the-gradient polarization, and successfully countered inward auxin flow until convergence points formed. These results point to an active contribution of L2 PINs to auxin patterning in the epidermal layer. This is in line with observations and simulations by Bayer et al, [10], who also reported upward PIN polarization in the L2. It should be noted, though, that the convergence pattern obtained with four cell layers is a bit blurred compared with the pattern obtained with a single layer.

Convergence points initiate downward veins
Using the same template and parameter values as in the previous simulation, we investigated vein initiation. We found that veins developed from convergence points (Fig 4 and S4 and S5 Movies). Very early in the formation of a convergence point, most PIN proteins were allocated to the bottom membrane of the cells in which auxin converged. Auxin flow was thus locally directed towards the L2 and a midvein was initiated.
It turned out that the initiation of veins was possible only if auxin importer concentration in inner cell layers was at least as high as in the L1. This is in contrast to experimental results by Kierzkowski et al. [8] and Bainbridge et al. [36], who did not observe auxin influx carrier AUX1 nor LAX1 outside the L1. This suggests the existence of other auxin importers than AUX1 and LAX1 in the inner tissue.

Developing veins are only weakly attracted by auxin sinks
In order to investigate how veins are progressing in the inner tissue and potentially connect to the pre-existing vasculature, we modified the previous template by assigning three evenlyspaced cells in the bottom layer (L4) as sinks. Sink cells had a high auxin turnover rate, which was intended to mimic the effect of a functional vein draining auxin. Thus, each sink cell induced a local gradient of auxin concentration centered around itself.
In our simulations, when a cell was exporting auxin towards a neighboring cell in the inner tissue, the auxin concentration in the latter cell was increasing until it started to polarize in turn towards a third cell, and so on. Veins developed this way, from convergence points towards regions with lower auxin concentration. Therefore, they progressed globally downward, away from L1 auxin-producing cells (Fig 5 and S6, S7 and S8 Movies). Often, a vein did not develop as a single sequence of cells, but instead as a bundle of such parallel sequences. Here, lateral diffusion of PIN proteins makes a clear difference. Without lateral diffusion of PINs, each vein formed as a single cell sequence, and this sequence followed preferentially paths connecting cells through very small apoplast element. This was the same effect we had seen in convergence point formation in the epidermal layer, where small membranes were  favored over larger ones. Again, this did not fit observations, so lateral diffusion of PINs appears as a necessary feature for modeling venation in a 3D irregular tissue.
Since developing veins globally progressed towards lower auxin concentrations, and since sink cells laid at the center of low auxin concentration regions, it was expected that all veins were going to eventually reach a sink. However, local auxin gradients surrounding sinks were very shallow and only had a short-range attraction power on developing veins. Veins which did not pass close to a sink missed their target and got lost in the tissue (Fig 5C and S8 Movie). Therefore, although Cieslak's mechanism does reproduce several features of vein formation in a 3D tissue, it does not provide an efficient way to find sinks.

Assuming a vein attraction factor (VAF) explains connection of developing veins to already existing vasculature
The poor capacity of developing veins to find sinks in our simulation is suggestive of an additional mechanism for guiding veins towards the already existing vasculature. Bayer et al. [10] reached a similar conclusion. This mechanism could take the form of facilitated diffusion of auxin through plasmodesmata, which has been proven by Smith and Bayer [37] to reliably and robustly connect sources to sinks. Both facilitated diffusion and polar transport could operate in parallel, with facilitated diffusion first defining a preferred route to the closest sink, and polar auxin transport then establishing polarization along this route. But there is little experimental support for facilitated diffusion of auxin in plant tissue. Therefore, following Bayer et al. [10], we resorted to an alternative model which posits a vein attraction factor (VAF) emitted by existing vasculature. In the original model by Bayer et al. [10], the vein attraction factor diffuses from cell to cell. Every cell can measure VAF concentrations in its neighbors and biases PIN allocation towards the membrane elements facing the cells with highest VAF concentrations. In Cieslak's model, however, there is no such remote sensing between cells. Its biochemical plausibility rests on the fact that all reactions occur locally, between a cytoplasm and a membrane, or between a membrane and neighboring apoplast. We propose an alternative VAF behavior, which is compatible with the locality of Cieslak's model. In our scheme, sink cells release VAF in their neighboring apoplast at a fixed rate. The VAF diffuses within the apoplast continuum, never reentering cells. However, it can bind to membranes and unbind from them, at some fixed rates. The concentration of bound VAF on a membrane element favors PIN exocytosis toward this membrane element.
We implemented this mechanism (see details in S1 Appendix) and found in our simulations that VAF gradients could establish and set a local polarity around every sink cell. The spatial range of VAF-induced polarity depends on VAF production rate and diffusion coefficient. If these values were adequately set, every initiating vein eventually headed and connected to the nearest sink cell (Fig 6 and S9, S10 and S11 Movies). The VAF hypothesis introduced an additional feature to vein development. Although the cells surrounding a sink had, at the beginning, low auxin and PIN contents, they rapidly polarized towards their neighbor sink. As time went by and VAF diffused farther away, second-order neighbor cells started to polarize towards the sink. When an initiating vein reached a cell already polarized towards a sink, this polarization got reinforced thanks to the high influx of auxin from the vein and the PIN production induced. PIN polarization was then firmly established between the source and the sink. To sum up, the complete vein was the result of the encounter of two opposed movements: the progression of an initiating vein towards a sink, and the expansion of a polarized region centered around this sink. In the initiating vein, cells had high auxin concentrations since each cell had to reach a threshold concentration before it could polarize towards another cell and made the vein progress. But, when the initiating vein reached the polarized region near the sink, auxin could flow rapidly to the sink, without accumulating. Thus, at the moment the polarized path was established, cells in the upper part of the vein had high auxin concentrations while cells closer to the sinks had low concentrations. The   former cells kept their high concentration whereas the latter cells only very slowly increased their concentration.

The model accounts for pin1 mutant phenotype and the effect of primordium ablation
To gain further support for the model, we tried to reproduce defective or artificially disturbed phenotypes. As a first example, the pin1 phenotype is characterized by the absence of primordia [38], In the model, we mimicked this mutation by a two-fold reduction in PIN production (both constitutive and auxin-dependent). Although this alteration seemed moderate, it turned out to be sufficient to inhibit primordium formation (S12 Movie). Most residual PINs in the L1 were oriented toward the L2, but without formation of distinct veins.
As a second example, we tried to recreate the effect of the laser ablation of an incipient primordium. In such experiments, a new primordium forms in the vicinity of the ablated site [39]. In the model, we removed L1 cells of a newly formed primordium and we could observe the rapid emergence of a new primordium next to the ablated one (Fig 7 and S13 Movie).
However, unlike what has been reported in ablation experiments [22,24], we did not observe a strong PIN polarization away from the wound. Such a wounding-induced PIN repolarization pattern is believed to be due to altered stress distribution and mediated through microtubules [22,24].

Discussion
The shoot apical meristem is a challenging system in our understanding of auxin patterning due to the overlap in time and space of different auxin transport regimes in an irreducible three-dimensional geometry. It poses both a conceptual and a technical problem. From a conceptual point of view, the co-localization of up-the-gradient and with-the-flux polarizations in the same cells, with the same efflux carriers, strongly suggests the existence of a common mechanism, while most models so far have addressed each polarization regime separately with incompatible hypotheses. From a technical point of view, building a 3D computational model of plant tissue including both symplast and apoplast, plus polar auxin transport, is a challenge which had never been taken up. As a significant first step towards that objective, we built a model of simplified meristematic tissue on which a biochemically plausible mechanism of auxin transport was implemented. According to this mechanism, the switch between the two polarization regimes depends on the auxin concentration. This approach reproduced key features of 3D phyllotactic patterning and offered new insights in the dynamics of PIN1 polarization in the shoot apical meristems.
We emphasize the importance of membrane-bound PIN diffusion in auxin patterning, for both convergence point and canalization. This has been overlooked in 2D models, in which the sizes of the discrete membrane elements are approximately of the same order of magnitude, and usually quite large. But in 3D geometry, the higher complexity of the topology results in a very wide range of areas for discrete membrane elements, with occasionally extremely small values. This compartmentalization becomes then overly artificial and ignores the intrinsic continuity of the lipid cell membrane in which bound PIN proteins can move laterally. It leads to aberrant accumulation of PINs on a few very tiny membrane elements. Assuming lateral diffusion of PIN proteins partially overcomes this issue by smoothing the barriers between neighbor discretized elements of the same cell membrane. It should be noted, however, that membrane-bound PINs display a more complex behavior, due to mechanisms limiting their lateral movements [40]. They tend to form clusters in which lateral mobility is strongly reduced. Only a minor fraction of membrane-bound PINs are unclustered and laterally mobile [41]. Future models could investigate in more details the significance of PIN clustering for cell polarization at tissue scale.
We demonstrated that inner cell layers, especially the L2, significantly contribute to auxin patterning in the L1. The higher auxin biosynthesis rate in the epidermal layer is not sufficient by itself to reach the critical concentration required for the formation of convergence points. In the early stages of phyllotactic patterning, auxin has to be kept confined in the L1. We showed that this confinement can be performed by up-the-gradient polarization in the L2, where PIN proteins polarize upward, towards auxin-rich L1 cells. A necessary condition is that auxin-dependent PIN production is fast enough to react quickly to the first auxin leaks into the L2. Bayer et al. [10] assumed reduced symplasmic communication between L1 and inner tissues, based on experimental evidence [42]. We did not make that assumption and found that PIN polarization alone can efficiently confine auxin in the L1 until midvein initiation. We also found that mutual influences between cell layers have a slight but non-negligible impact on the final pattern. This suggests that precise and realistic descriptions of surface patterning should take into account effects from the underlying tissue.
It turned out that midveins cannot develop if the concentration of auxin influx carriers is lower in the inner tissue than in the L1. This seemingly contradicts experimental results [8,36] showing that auxin influx carriers AUX1 and LAX1 are almost exclusively present in the epidermal layer. But at the same time, Bainbridge et al. [36] reported the expression of another auxin importer, LAX2, in the initiating vein. LAX2 has also been shown to regulate vascular patterning in cotyledons [43]. Our model does not exclude the possibility that auxin influx carriers present in the inner tissue could be restricted to the cells forming the initiating veins. However, in this case, they have to be expressed from the very beginning of vein initiation to make it possible. Since the mechanism inducing LAX2 expression is not known, amending the model in this direction seems both speculative and premature. Our investigation also pointed out that auxin gradients induced by the existing vascular system are too shallow to efficiently attract initiating veins. This fact had been already noted in 2D simulations performed by Bayer et al. [10], and is linked to the more general problem of establishing polarity within a tissue. Various models of tissue polarity have been developed (see a comparative study in [20]). In the shoot apical meristem of Arabidopsis, the region in which a polarity has to be established to guide an initiating vein is relatively small. Therefore, a simple mechanism based on a putative diffusive molecule is sufficient, as already outlined by Bayer et al. [10]. We proposed a detailed process for the action of such a Vein Attraction Factor. We assumed that the VAF is transported apoplastically, which seemed to us more realistic than the symplastic transport first proposed by Bayer et al. [10].
This vein attraction mechanism raises new questions about source-sink connections. Classically, an initiating vein is viewed as progressing in a tissue from a source until it reaches a sink. In our simulations, each sink builds up an expanding polarity field centered around itself. When the progressing vein meets the expanding polarity field, the connection is almost fully established. Then, the vein progresses very quickly and directly to the sink, locally increasing polarity in its wake. The relative contribution of each process-free progression of the initiating vein and expansion of the polarity field-depends on the biosynthesis rate and diffusion coefficient of the VAF, but also on the precise timing of vein initiation and VAF release. Further experimental studies are needed to improve our understanding of how primordia connect to the vascular system. Modeling efforts are hindered by the unknown nature of the putative VAF.
Another obvious limitation to our model is that it does not reproduce phyllotactic patterns. This is essentially because it does not include tissue growth, which plays a crucial role in the dynamics of phyllotactic patterning. However, simulating meristem growth with mechanical processes would add another layer of technical complexity and requires much more development and computational resources. Making the size and shape of the template closer to an actual meristem is also a necessary improvement for future models to gain more realism.
It can also be argued that the PIN polarization mechanism used in our model is questionable since there is no evidence for auxin complexes. Yet, the core of the mechanism is more generic than it looks at first sight. Cieslak et al. [27] proposed various other implementations of it, reflecting different biochemical assumptions. Since the precise reactions underlying the feedback of auxin transport on PIN localization are unknown, it is difficult to discriminate between these implementations. Other ones can be designed based on new experimental findings. For instance, the fact that convergence point formation and vein initiation seem to be associated with two distinct groups of auxin influx carriers could be exploited to amend the current reaction network. Other ingredients could be introduced, such as the inhibiting effects of auxin on PIN endocytosis [44] and PIN vacuolar degradation [45,46]. Alternatives will probably emerge in the near future from experimental studies, and 3D models like ours will be valuable tools to thoroughly assess their explanatory powers in terms of auxin pattern formation.