Spatial regularity control of phyllotaxis pattern generated by the mutual interaction between auxin and PIN1

Phyllotaxis, the arrangement of leaves on a plant stem, is well known because of its beautiful geometric configuration, which is derived from the constant spacing between leaf primordia. This phyllotaxis is established by mutual interaction between a diffusible plant hormone auxin and its efflux carrier PIN1, which cooperatively generate a regular pattern of auxin maxima, small regions with high auxin concentrations, leading to leaf primordia. However, the molecular mechanism of the regular pattern of auxin maxima is still largely unknown. To better understand how the phyllotaxis pattern is controlled, we investigated mathematical models based on the auxin–PIN1 interaction through linear stability analysis and numerical simulations, focusing on the spatial regularity control of auxin maxima. As in previous reports, we first confirmed that this spatial regularity can be reproduced by a highly simplified and abstract model. However, this model lacks the extracellular region and is not appropriate for considering the molecular mechanism. Thus, we investigated how auxin maxima patterns are affected under more realistic conditions. We found that the spatial regularity is eliminated by introducing the extracellular region, even in the presence of direct diffusion between cells or between extracellular spaces, and this strongly suggests the existence of an unknown molecular mechanism. To unravel this mechanism, we assumed a diffusible molecule to verify various feedback interactions with auxin–PIN1 dynamics. We revealed that regular patterns can be restored by a diffusible molecule that mediates the signaling from auxin to PIN1 polarization. Furthermore, as in the one-dimensional case, similar results are observed in the two-dimensional space. These results provide a great insight into the theoretical and molecular basis for understanding the phyllotaxis pattern. Our theoretical analysis strongly predicts a diffusible molecule that is pivotal for the phyllotaxis pattern but is yet to be determined experimentally.


Introduction
Living organisms often form periodic patterns with spatial regularity in a self-organizing manner [1,2]. One such well-known example is phyllotaxis, the arrangement of leaves on a plant stem. The phyllotaxis exhibits various types of patterns depending on the plant species and this has attracted many scientists because of its beautiful geometric configuration [3]. Phyllotaxis is originated at the shoot meristem, in which leaf primordia are periodically formed by maintaining a constant distance from each other [4][5][6][7]. This spatial regularity is established by the mutual interaction between a mobile plant hormone auxin and its efflux carrier membrane protein PIN1, which cooperatively generate small regions with high auxin concentrations called auxin maxima that are involved in leaf primordia. In the process of the auxin maxima formation, auxin accumulates at the presumptive position of a future primordium while PIN1 is polarized toward the center position [8][9][10][11]. According to this experimental observation, the auxin maxima pattern is often explained by the concept of "up-the-gradient" in which auxin is transported by PIN1 against its own gradient while PIN1 is polarized toward higher auxin [11][12][13].
Based on this concept, a class of mathematical models (corresponding to Model O in this paper) has been proposed, in which PIN1 is localized to a cell membrane depending on the auxin concentration of neighboring cells. Because these models can successfully reproduce the spatial regularity of auxin maxima and various phyllotaxis types, they are excellent models for understanding the nature of phyllotaxis pattern formation [11,12,[14][15][16]. On the other hand, they are highly simplified and abstract models, and could have various problems when considering the molecular mechanism. With respect to spatial structure, these can be distilled into two major points. First, extracellular space (i.e., apoplast space) is absent in these models, in which auxin moves directly between cells by PIN1 and diffusion. However, in plant tissues, auxin is transported between cytoplasm and apoplast because cells do not contact each other but are separated by apoplast space. Second, it remains unclear how cells sense auxin concentrations of neighboring but separated cells for PIN1 polarization.
Sahlin et al. [16] showed that, despite the extracellular space, self-organized patterns can be generated by the "up-the-gradient" concept that PIN1 polarization depends on the auxin concentration of neighboring cells. However, it has still not been clarified how the information on auxin concentration is transmitted between neighboring cells. Conversely, although the apoplast space is also considered by Webnik et al. [17] and Cieslak et al. [18], these reports describe the canalization pattern generated by the "with-the-flux" concept, which is different from auxin maxima by the "up-the-gradient" concept.
In this paper, therefore, we investigated how the spatial regularity of phyllotaxis pattern is controlled under realistic conditions, that is, in the presence of extracellular space. We first confirmed that the introduction of extracellular space has a disruptive effect on the spatial regularity in the conventional model (Model O), even in the presence of direct auxin diffusion between cytoplasm or apoplast spaces. This result strongly suggests that an unknown molecular mechanism is required for phyllotaxis pattern formation. We also found that the spatial regularity can be restored by assuming a diffusible molecule that mediates the feedback signaling from auxin to PIN1 polarization. This theoretical analysis strongly predicts a diffusible molecule that is critical for phyllotaxis pattern but remains to be found.

Models
Cells are tightly arranged in a one-or two-dimensional space. Auxin (Models O, A, and B) and a hypothesized molecule X (Model B) are uniformly distributed in a cell, and their concentrations in cell i are denoted by a i and x i , respectively (Fig 1). Auxin efflux carrier PIN1 is unevenly distributed to the cell membrane, and its density in the membrane of cell i toward neighboring cell j is denoted by p i,j . Models A and B consider apoplast (i, j), the extracellular In Model A, auxin is transported between cytoplasm and apoplast by PIN1 and influx carrier, while PIN1 is polarized depending on auxin concentration of neighboring apoplast spaces. (C) In the framework of Model B, an assumed molecule X is incorporated into Model A. Molecule X is expressed in response to cytosolic auxin and diffuses freely between cytoplasm and apoplast. (D) Model B considers various feedback regulations from molecule X to auxin-PIN1 dynamics. (E) In Model B6, PIN1 is polarized depending on X concentration of neighboring apoplast spaces, instead of auxin. (F) In addition to simple diffusion between cytoplasm and apoplast, Models A and B consider direct diffusions between neighboring cells (symplast diffusion) and between neighboring apoplast spaces (apoplast diffusion). (G) Schematic representation of L 1 and L 2 : indices for spatial scale of auxin maxima pattern in numerical simulations. a i (or x i ) and a 0 i;j (or x 0 i;j ) are auxin (or X) concentrations in cell i and apoplast (i, j), respectively. p i,j is PIN1 density of the membrane toward cell j in cell i. space between neighboring cells i and j, and concentrations of auxin and molecule X in apoplast (i, j) are denoted by a 0 i;j ð¼ a 0 j;i Þ and x 0 i;j ð¼ x 0 j;i Þ, respectively.

Model O
Change of auxin concentration of cell i (a i ) is described by where cell j is a neighbor of cell i, A is related to the synthesis rate, G a is the degradation rate, D a is the diffusion coefficient, E p is the efficiency of PIN1 efflux carrier, and f i,j (= −f j,i ) is the net flow of auxin by PIN1 from cell i to cell j, consisting of auxin efflux and influx. Auxin is constantly synthesized and degraded at a constant rate (the first term of the right-hand side of Eq 1), is transported by PIN1 (the second term), and diffuses between neighboring cells (the third term). On the other hand, change of PIN1 density (p i,j ) is described by where G p is the degradation rate, K is the number of neighboring cells, p is a constant related to PIN1 density, and φ 0 (a j ) is the regulatory function for PIN1 polarization (Fig 1A). PIN1 is localized to cell membrane depending on the auxin concentration of neighboring cells and is degraded at a constant rate. The total PIN1 amount of cell i, P i ∑ j p i,j , satisfies dP i /dt = G p (Kp − P i ), indicating that Kp is the stable equilibrium of P i . Thus, equilibria of a i and p i,j are given respectively by a eq ¼ A and p eq ¼ p ð4Þ When G p is sufficiently large, p i,j quickly approaches equilibrium: Therefore, Eq 5 can be used instead of Eq 3 in a simplified version of Model O. Model O in this study is equivalent to models reported previously in references [11] and [12]. Model equations of the simplest form in [11] can be described by where T, D, and P are constants. Eqs 6 and 7 are identical to the simplified version of Model O (Eqs 1, 2 and 5) with G a = 0, D a = D, E p = T, p = P/K, and φ 0 (a j ) = a j . Conversely, equations used in [12] are somewhat complicated compared to those in [11]. However, phyllotactic patterns can be generated under the condition of fixed total PIN1 concentration (i.e., [PIN] i is constant in Eq 2 of [12]), no saturation of auxin synthesis (i.e., κ IAA = 0 in Eq 5), and the linear dependence of the flux on auxin concentration (i.e., replacement of ½PIN 2 i and ½PIN 2 j by [PIN] i and [PIN] j , respectively, in Eq 3). In addition to these conditions, by considering a regular cell lattice (i.e., cell side length l i!j is constant in Eq 2), model equations can be simplified by where ρ IAA , μ IAA , D, T, κ T , P, and b are constants. Eqs 8-10 are identical to the simplified Model O (Eqs 1, 2 and 5) with G a = μ IAA , A = ρ IAA /μ IAA , D a = D, E p = T, p = P/K, and φ 0 ða j Þ ¼ b a j , except for the saturation effect of auxin from neighboring cells on the flux in Eq 9. This effect would negatively affect the pattern formation in a manner that the saturation effect becomes strong and accordingly the pattern tends to disappear as κ T increases. Therefore, this effect is not essential for generating a phyllotactic pattern, indicating that Eqs 8-10 are equivalent to Model O.

Model A
Incorporation of extracellular region. Model A was constructed by incorporating apoplast (i, j), the extracellular space between neighboring cells i and j, into Model O (Fig 1B). Auxin concentration of apoplast (i, j) is denoted by a 0 i;j ð¼ a 0 j;i Þ. Changes of cytosolic auxin (a i ) and apoplast auxin (a 0 i;j ) are described by where G a , A, and E p have the same notations as in Eqs 1 and 2, D a is the diffusion coefficient between cytoplasm and apoplast, V is the volume ratio of apoplast to cytoplasm, E q is the efficiency of influx carrier function, q is influx carrier density of a cell side, and f i,j is the auxin flow from cell i to apoplast (i, j), consisting of efflux by PIN1 and influx by auxin influx carrier. Auxin is synthesized constantly in cytoplasm and degraded at a constant rate (the first terms of the right-hand sides of Eqs 11 and 12), is transported by carriers (the second terms), and diffuses between cytoplasm and apoplast (the third terms). On the other hand, the change of PIN1 density (p i,j ) is described by where G p , K, and p have the same notations as in Eq 3, and φ a ða 0 i;j Þ is the regulatory function for PIN1 polarization. PIN1 is localized to cell membrane depending on the auxin concentration of neighboring apoplast spaces and is degraded at a constant rate ( Fig 1B). As in Model O, Kp is the stable equilibrium of the total PIN1 amount of a cell. Equilibria of a i , a 0 i;j , and p i,j are given respectively by where a 0 A/(KV(E p p + D a ) + 2(E q q + D a ) + VG a ). Effect of symplast or apoplast diffusion of auxin. In addition to simple diffusion between cytoplasm and apoplast as described above, signal molecules in plants have two major diffusion types ( Fig 1F). One is symplast diffusion, which is the direct diffusion between cells via narrow tube-like structures called plasmodesmata, through which small molecules including small RNAs and transcription factors can migrate between neighboring cells [19][20][21]. Auxin is a small signal molecule and is reported to pass through plasmodesmata [22][23][24]. The other type is apoplast diffusion by which signal molecules such as secreted peptides can freely move among extracellular spaces because they are connected to each other in plant tissues [25,26]. Thus, we investigated the effect of the symplast or apoplast diffusion of auxin on pattern formation.
In Model A, the symplast diffusion of auxin is incorporated by replacing Eq 11 with where D a1 is the diffusion coefficient between neighboring cells. On the other hand, the apoplast diffusion of auxin is incorporated by replacing Eq 12 with where apoplast (k, l) is a neighbor of apoplast (i, j), and D a2 is the diffusion coefficient between neighboring apoplast spaces.

Model B
Incorporation of diffusible molecule. Model B was constructed by incorporating an assumed diffusible molecule X into Model A ( Fig 1C). As with auxin, concentrations of molecule X in cell i and apoplast (i, j) are denoted by x i and x 0 i;j ð¼ x 0 j;i Þ, respectively. Changes of x i and x 0 i;j are described respectively by where G x is the degradation rate, D x is the diffusion coefficient between cytoplasm, V is the volume ratio of apoplast to cytoplasm, and θ(a i ) is the regulatory function of auxin on X synthesis. Molecule X is synthesized depending on cytosolic auxin and degraded at a constant rate (the first terms of the right-hand sides of Eqs 18 and 19) and diffuses between cytoplasm and apoplast (the second terms) (Fig 1C). We used Eqs 11-14, 18 and 19 as the framework of Model B. Equilibria of x i and x 0 i;j are given respectively by where x 0 θ(a eq )/((KV + 2)D x + VG x ).

Feedback regulations by diffusible molecule X.
We incorporated various feedback regulations from molecule X to auxin-PIN1 dynamics into the Model B framework (Fig 1D; Eqs 11-14, 18 and 19) by replacements as follows: (Model B1) Effect of cytosolic X (x i ) on influx carrier amount (q) is incorporated by the replacement of q ! ψ 1 (x i )q in Eq 13: (Model B6) Effect of apoplast X (x 0 i;j ) on PIN1 localization to cell membrane is incorporated by the replacement of φ a ða 0 i;j Þ ! φ a ða 0 i;j Þφ x ðx 0 i;j Þ in Eq 14: and φ x ðx 0 i;j Þ are regulatory functions that depend on molecule X. The equations and regulatory functions used in numerical simulations are summarized in S1 Table. Effect of symplast diffusion of molecule X. In Model B, we examined the symplast diffusion of molecule X (i.e., direct diffusion between neighboring cells; Fig 1F), instead of the simple diffusion between cytoplasm and apoplast, by replacing Eqs 18 and 19 with where G x and θ(a i ) have the same notations as in Eqs 18 and 19, and D x1 is the diffusion coefficient between cells. Molecule X is synthesized depending on cytosolic auxin and degraded at a constant rate (the first terms of the right-hand sides of Eqs 27 and 28) and diffuses between cells (the second term of Eq 27). Equilibria of x i and x 0 i;j are given respectively by x eq ¼ yða eq Þ and x 0 eq ¼ 0 ð29Þ Effect of apoplast diffusion of molecule. We also examined the apoplast diffusion of molecule X (i.e., direct diffusion between neighboring apoplast spaces; Fig 1F), instead of the simple diffusion, by replacing Eqs 18 and 19 with where G x , V, and θ(a i ) have the same notations as in Eqs 18 and 19, K is the number of neighboring cells, S is the secretion coefficient, D x2 is the diffusion coefficient between neighboring apoplast spaces, and apoplast (k, l) is a neighbor of apoplast (i, j). Molecule X is synthesized depending on cytosolic auxin and degraded at a constant rate (the first terms of the right-hand sides of Eqs 30 and 31), is secreted to apoplast spaces (the second terms), and diffuses between apoplast spaces (the third term of Eq 31). Equilibria of x i and x 0 i;j are given respectively by

Numerical simulations
We used one-dimensional arrays of N = 200, 50, or 40 cells and two-dimensional sheets of 20 × 20 or 14 × 14 hexagonal cells in the numerical simulations. Initial values of auxin, PIN1, and molecule X are given by their equilibrium with 1.0% fluctuation. The numerical simulations were performed using the Euler method with time step Δt = 0.001 under the periodic boundary condition. Equations and regulatory functions used are summarized in S1 Table. Parameter values used are described in figure legends.

Index of spatial scale
To evaluate the spatial scale of auxin patterns generated by the numerical simulations, we used wavelength of auxin maxima (L 1 ) and average size of auxin maximum (L 2 ) as indices of the spatial scale (Fig 1G). In a one-dimensional array of N cells in the periodic boundary condition, auxin concentration of cell n or apoplast (n, n+1) is denoted here by c n , where n = 0,1,Á Á Á, N and c 0 c N . Wavelength of auxin maxima pattern (L 1 ). After discrete Fourier transform of auxin concentration c n : F k ð Þ ¼ 1 showing the largest spectral intensity of |F(k)| 2 and then used the corresponding wavelength defined by as an index of the spatial scale. Average size of auxin maximum (L 2 ). A cell or an apoplast space with c n > " c is called an "auxin spot" and a successive string of auxin spots is denoted by an "auxin cluster" where " c is the average of c n . The average size of auxin maximum (L 2 ) is defined as the total number of auxin spots divided by that of auxin clusters.

Model O
The phyllotaxis pattern of auxin maxima has been explained by a class of simplified mathematical models based on the feedback dynamics between auxin and PIN1. Because these models do not consider extracellular region (i.e., apoplast), auxin moves directly between cells by PIN1-dependent directional transport and passive diffusion to change its distribution ( Fig  1A). On the other hand, PIN1 is asymmetrically localized to the cell membrane, preferentially toward neighboring cells with high auxin concentrations. We used Model O (Eqs 1-3), which is one of the simplest representations of such dynamics, to examine the spatial regularity control of auxin maxima pattern. As in previous reports [11,12], we confirmed that Model O can form auxin maxima patterns with spatial regularity, an essential characteristic of phyllotaxis, focusing on its spatial scale.
Spatial regularity control in Model O. We considered a one-dimensional array of N cells under the periodic boundary condition in Model O (Eqs 1-3 with K = 2), and performed a linear stability analysis of the equilibrium (see S1 Text (ii) for detail). This theoretical analysis shows that eigenvalue λ k , which is associated with the growth rate of the pattern with wavenumber k(= 0,1,Á Á Á,N − 1), is given by where ν cos(2πk/N) 2 [−1,1], c 0 −(G a + 2c 1 + 2c 2 ), c 1 E p p + D a , and c 2 À E p pa eq φ 0 0 ða eq Þ= 2φ 0 ða eq Þ. The condition for non-uniform patterns is described by where νÃ = −c 1 /4c 2 . When Eq 35 is satisfied (i.e., spatial homogeneity is broken), the spatial scale (i.e., wavelength LÃ) of the pattern with the highest growth rate depends on parameter values and is given by where LÃ increases as νÃ becomes large.
In the case of the regulatory function for PIN1 polarization φ 0 (a j ) = a j n , Eq 35 and νÃ become D a < ð2n À 1ÞE p p À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2nG a E p p q and ð37Þ respectively, where R a D a /E p p corresponds to the strength ratio of auxin diffusion to transport by PIN1. Eigenvalue λ k (ν) (Eq 34) and the condition for pattern formation (Eq 37) are shown in Fig 2A and 2B, respectively. This result predicts that the spatial scale of formed patterns LÃ becomes large as the diffusion coefficient (D a ) increases or auxin transport by PIN1 (E p p) and the regulatory strength of PIN1 polarization (n) decrease. This prediction is supported by numerical simulations, in which non-uniform auxin distribution is observed in the parameter area corresponding to Eq 37 (Fig 2B), and the wavelength of auxin pattern (L 1 ) and average size of auxin maximum (L 2 ) increase by decreasing p or increasing D a (Fig 2C-2H). Besides, Model O can also generate regular patterns of auxin maxima in the two-dimensional space (Fig 3A). These results are consistent with previous reports [11,12]. Therefore, Model O can reproduce the regular distance between auxin maxima, which is an essential characteristic of phyllotaxis pattern, suggesting that "up-the-gradient" is a central concept of this pattern formation.

Model A
Incorporation of extracellular space. While Model O is suitable for understanding the "up-the-gradient" concept as described in the previous section, it is not appropriate when considering cellular spatial structures because plant cells are separated from each other by the extracellular region (i.e., apoplast). Thus, we investigated how the introduction of extracellular region affects the spatial regularity of auxin maxima. In the revised model called Model A (Eqs [11][12][13][14], auxin moves between cytoplasm and apoplast; outwardly by PIN1, inwardly by influx carriers, and non-directionally by passive diffusion, instead of the direct migration between cells in Model O (Fig 1B). On the other hand, PIN1 is asymmetrically localized to the cell membrane depending on the auxin concentration of neighboring apoplast spaces, instead of that of neighboring cells in Model O. Spatial regularity control in Model A. We consider a one-dimensional array of N cells that are separated from each other by apoplast space under the periodic boundary condition in Model A (Eqs 11-14 with K = 2), and we performed a linear stability analysis of the equilibrium as in Model O (see S1 Text (iii) for detail). This theoretical analysis shows that eigenvalue λ k , which is associated with the growth rate of the pattern with wavenumber k(= 0,1,Á Á Á,N −1), is given by where ν cos(2πk/N) 2 [−1,1], c 0 α + β − (2(E q q + D a )/V + G a ), c 1 (α − β)/2, α 2(E p p + D a )(E q q + D a )/V(2E p p + 2D a + G a ), and b E p pð2ðE q q þ D a Þa 0 eq þ G a AÞφ 0 a ða 0 eq Þ= Vð2E p p þ 2D a þ G a Þφ a ða 0 eq Þ. It is also shown that the condition for generating spatial heterogeneity is described by When Eq 40 is satisfied (i.e., spatial homogeneity is broken), the spatial scale (i.e., wavelength LÃ) of the pattern with the highest growth rate (i.e., the largest λ k ) is independent of parameter values and always given by This theoretical result indicates that the pattern of apoplast spaces alternating between high and low auxin concentrations always grows fastest compared to that with longer spatial scales and is consistent with the result reported by Sahlin et al. [16]. This also indicates that the extracellular space has a destructive effect on the spatial regularity of auxin maxima.
In the case of regulatory function for PIN1 polarization φ a ða 0 i;j Þ ¼ ða 0 i;j Þ n , Eq 40 becomes Eigenvalue λ k (ν) (Eq 39) and the condition for pattern formation (Eq 42) are shown in Fig 4A  and 4B, respectively. These theoretical results are consistent with numerical simulations, in which a non-uniform distribution of auxin occurs in the corresponding parameter region of Eq 42 (Fig 4B-4D). On the other hand, opposite to Model O, the spatial scale of formed patterns becomes extremely small in most parameter conditions examined. In addition, similar Spatial regularity control of phyllotaxis pattern results are observed in the case of two-dimensional space, in which spatially regular patterns of auxin maxima seen in Model O are completely eliminated in Model A (Fig 3B). Effect of symplast diffusion of auxin. In Model A described, we consider auxin diffusion between cytoplasm and apoplast (Fig 1F, simple diffusion). However, because auxin can directly move between cells through plasmodesmata [22][23][24] (Fig 1F, symplast diffusion), we examined the effect of the symplast diffusion on pattern formation by replacing Eq 11 with Eq 16. Numerical simulations show that patterns with larger spatial scales cannot be recovered under this diffusion condition (Fig 5A and 5B).
Effect of apoplast diffusion of auxin. In plant tissues, the apoplast spaces are connected to each other and auxin can move freely among them (Fig 1F, apoplast diffusion). Thus, we also examined the effect of the apoplast diffusion by replacing Eq 12 with Eq 17. Linear stability analysis predicts that this change causes no essential effects on the spatial regularity of the auxin maxima pattern (S1 Text (iv)). In fact, this theoretical prediction is supported by numerical simulations, in which the spatial scale of the generated patterns remains extremely small as in Fig 4 (Fig 5C and 5D).

Model B
Incorporation of diffusible molecule. The previous section showed that extracellular space has a disruptive effect on the spatial regularity of auxin maxima. This result strongly suggests that phyllotaxis pattern cannot be explained by the "up-the-gradient" framework alone but requires an unknown molecular mechanism. On the other hand, the spatial regularity can be generated even in the presence of extracellular space if PIN1 polarization depends on the auxin concentration of neighboring cells [16], suggesting a mechanism that transmits the auxin concentration between neighboring cells. This could be fulfilled by considering a factor  Table. https://doi.org/10.1371/journal.pcbi.1006065.g005 that is induced by auxin and diffuses freely. To verify this possibility, we assumed a molecule X that is expressed depending on auxin concentration and diffuses between cytoplasm and apoplast (Fig 1C). The framework of the revised model called Model B is obtained by incorporating molecule X into Model A, and then we considered various feedback effects of X on auxin-PIN1 dynamics as follows (Fig 1D):  Table. To examine whether each feedback regulation can restore the spatial regularity of the auxin pattern as in Model O, we carried out numerical simulations with systematically changing regulatory strengths of apoplast auxin on PIN1 polarization (n) and of molecule X (m) (Fig 6). In the one-dimensional space, we found three conditions for regular pattern formation: negative values of m in Models B2 (Fig 6C, 6D and 6N) and B5 (Fig 6I, 6J, and 6N) and positive values of m in Model B6 (Fig 6K, 6L, and 6O). However, the former two do not appear to be plausible because no regular patterns are formed under the corresponding conditions in the two-dimensional space (S1 and S2 Figs). By contrast, in both one-and two-dimensional spaces, only Model B6 generates auxin maxima with relatively large wavelengths, which depend on parameter values (Fig 5K and 5L). This result suggests that the feedback via PIN1 polarization is essential for the auxin maxima pattern.

Effect of diffusion variation in Model B.
In addition to the diffusion between cytoplasm and apoplast, signal molecules in plants have two major diffusion types (Fig 1F): (i) direct diffusion between cells (symplast diffusion) [19][20][21] and (ii) diffusion among extracellular spaces (apoplast diffusion) [25,26] as described in the above. We examined how these diffusion types affect the spatial regularity. Symplast or apoplast diffusion can be incorporated into Model B by replacing Eqs 18 and 19 with Eqs 27 and 28 or with Eqs 30 and 31, respectively (S1 Table). These diffusion types provide numerical simulations similar to those in the simple diffusion (Fig 6), in which normal auxin maxima patterns can be restored only in the feedback regulation of PIN1 polarization (Model B6) but not in the other conditions (Models B1-B5) (S3 Fig). This result reconfirms that the feedback from auxin to PIN1 polarization is crucial for phyllotaxis pattern whereas the diffusion type of the diffusible molecule is not. Next, we investigated Model B6 in detail.
Spatial regularity control in Model B6. In Model B6, regular patterns of auxin maxima can be generated even in the absence of the feedback from apoplast auxin to PIN1 polarization (Figs 5K and 5L and S3K and S3L; n = 0). Accordingly, this feedback regulation is not essential and is not considered in the following analysis (Fig 1E; Eqs 11-13, 18, 19 and 26 with K = 2 and φ a ða 0 i;j Þ ¼ 1). As with Model A, we consider a one-dimensional array of N cells and performed a linear stability analysis of the equilibrium (see S1 Text (v) for detail). By a similar argument to that of Model O, we can obtain approximate equations corresponding to Eqs 34-36: jn Ã j < 1 and l k ðn Ã Þ > 0 ð44Þ where ν cos(2πk/N) 2 [−1,1], νÃ % −c 1 /4c 2 , c 0 2c 1 − 2c 2 − (2E p p + 2D a + G a ), c 1 γ(E p p + D a ), c 2 À gkE p pa eq y 0 ða eq Þφ 0 x ðx 0 eq Þ=2φ x ðx 0 eq Þ; γ (E q q + D a )/(2E q q + 2D a + VG a ), and κ D x G x /(2D x + G x )(2D x + VG x ) Eqs 43-45 are associated with the pattern growth rate, condition for generating spatial heterogeneity, and spatial scale of generated patterns, respectively. This theoretical result indicates that spatial scale LÃ changes depending on parameter values as in Model O, and thus the spatial regularity of auxin pattern can be restored in Model B6.
where R a D a /E p p and R x D x /G x . Eigenvalue λ k (ν) (Eq 43) and the condition for pattern formation (Eq 44) are shown in Fig 7A and 7B-7D, respectively. Eqs 45 and 46 suggests that the spatial scale (LÃ) increases as diffusion coefficients of auxin (D a ) and molecule X (D x ) become large (Fig 7B), and this is consistent with numerical simulations (Fig 8A, 8B and 8G-8J). However, X diffusion is essential for pattern formation whereas auxin diffusion is not, because auxin maxima can be formed in the absence of auxin diffusion (D a = 0) but cannot without X diffusion (D x = 0) (S4 Fig). In plant tissues, the extracellular region is usually a very small space compared with cytoplasm, indicating that V, the volume ratio of apoplast to cytoplasm, is sufficiently small (i.e., V ( 1). Eq 46 also suggests that νÃ, accordingly LÃ, increases as V becomes small and has the limit: lim V!0 νÃ % 2(1 + R a )(1 + 2R x )/rm (Fig 7C). This is also consistent with numerical simulations (Fig 8C and 8D).
On the other hand, r and m are related to the regulatory strengths of auxin on X synthesis and of X on PIN1 polarization, respectively. Stable patterns are theoretically predicted and numerically generated in two separate parameter areas of r,m > 0 and r,m < 0 (Figs 7D, 8E and 8F), in which the two regulations are both stimulatory and inhibitory effects, respectively, leading to in-phase and anti-phase, respectively, patterns of auxin and molecule X (Fig 8K and  8L). This suggests that regular patterns are induced by stimulating PIN1 polarization by auxin while the way in which individual reactions are controlled is not so important.

Discussion
In phyllotaxis pattern formation, auxin maxima involved in leaf primordia are formed in a self-organizing manner while maintaining a constant distance from each other. However, the molecular mechanism generating the spatial regularity remains unclear. This spatial regularity has been explained by simple mathematical models (Model O; Figs 1A, 2 and 3A), in which PIN1 is polarized preferentially toward neighboring cells with higher auxin concentrations [11][12][13]. But these models have two major problems concerning spatial structure: one is the absence of extracellular space and the other is how cells perceive auxin concentrations of neighbors for PIN1 polarization.
In this study, therefore, we intensively investigated how the spatial regularity of the phyllotaxis pattern is controlled under appropriate conditions for plant cells. We showed theoretically and numerically that auxin maxima patterns with large spatial scale are completely eliminated by introducing extracellular space (Model A; Figs 1B, 3B, 4 and 5). This strongly suggests that phyllotaxis pattern requires an unknown molecular mechanism as well as auxin-PIN1 mutual interaction. Furthermore, we found that regular patterns can be restored by the simple and plausible assumption that a diffusible molecule mediates the feedback from auxin to PIN1 polarization (Model B6; Figs 3C and 8). Although we mostly investigated in onedimensional space, the same can be applied to the case of two-dimensional space. Model O can generate regular patterns of auxin maxima (Fig 3A). This spatial regularity is completely disrupted by considering extracellular space in Model A (Fig 3B), but is restored by introducing a diffusible molecule in Model B6 (Fig 3C). This diffusible molecule plays the role of transmitting auxin concentration to neighboring cells.
Auxin reportedly enhances the PIN1 localization at the cell membrane [27][28][29]. AUXIN-BINDING PROTEIN 1 (ABP1) might act as an apoplastic auxin receptor in the signaling pathway of PIN1 polarization although the function of ABP1 has recently contended [28,[30][31][32]. However, ABP1 probably makes no contribution to the phyllotaxis pattern formation regardless of whether it is an actual auxin sensor or not, because our study strongly suggests that the auxin maxima pattern cannot be established by the direct regulation of auxin on PIN1 polarization (Figs 4-6 and S1-S3).
Although our study showed that a diffusible factor can restore regular patterns that are disrupted by the presence of extracellular space, this finding does not rule out other possibilities for the spatial communication between cells. One such possible mechanism is mechanical force, including stress and strain, which affects the morphogenesis of plants and animals [33][34][35][36]. Mechanical force could stabilize the outgrowth of leaf primordia by feedback mechanism in which mechanical stress induces alignment of microtubules, enhancing cell elongation and primordial outgrowth, which reinforces the stress field [37]. In contrast, experimental evidence that mechanical force is involved during auxin maxima formation has not yet been obtained [4,7]. However, auxin could alter the mechanical properties of the extracellular matrix by inducing cell-wall loosening [5,6,38,39], suggesting that mechanical force may contribute to the pattern formation.
Our model predicts that the spatial scale of generated patterns (LÃ) becomes large by increasing νÃ, which follows νÃ / (1 + R a )/a eq |θ 0 (a eq )|, where R a D a /E p p and θ(a i ) is the regulatory function of auxin on X synthesis (Model B6; S1 Text, Eq 44). Therefore, if molecule X predicted in this paper exists, LÃ is affected by the regulatory activity of auxin on the expression of X as well as by the amount of PIN1. That is, under the conditions of a constant amount of PIN1 and constant activity of PIN1 recycling between cytosol and cell membrane, it is expected that, as the gene expression control by auxin becomes weak, the spacing between auxin maxima gradually increases and patterns suddenly disappear under a threshold of the control strength. This prediction could be used to experimentally validate whether or not our model is correct. Auxin affects the expression of many genes by cooperating with the TRANSPORT INHIBITOR RESISTANT 1/AUXIN SIGNALING F-BOX (TIR1/AFB) F-box proteins, the AUXIN/INDOLE-3-ACETIC ACID (Aux/IAA) transcriptional coregulators, and sequence-specific binding proteins called AUXIN RESPONSE FACTORs (ARFs) [40][41][42]. Because these factors are possible candidates that control the expression of molecule X, our model could be verified using plants showing various expression activities by genetically manipulating these factors.
Our study could predict a diffusible factor that is essential for phyllotaxis pattern but remains to be found. This factor(s) X must satisfy the following requirements: (i) X expression is regulated by auxin concentration.
(iii) X affects PIN1 polarization in the cell membrane.
(iv) The complete defect of X results in that of auxin maxima and accordingly that of leaf primordia.
Although such factors like X are not yet known, several diffusible molecules affecting PIN1 polarization have been reported. Strigolactone is a mobile plant hormone and cooperates with auxin to control shoot branching of plants. Auxin positively regulates the transcript of strigolactone biosynthesis genes and, in turn, strigolactone signaling triggers PIN1 depletion from the plasma membrane [43][44][45]. On the other hand, GOLVEN (GLV) genes encode small secretory peptides that are involved in root gravitropic responses and meristem organization in Arabidopsis. Transcription of GLV genes is rapidly induced by auxin, and the GLV peptide treatment stimulates the localization of auxin efflux carrier PIN2 at the cell membrane [46,47]. Another mobile plant hormone cytokinin plays important roles in various developmental events through crosstalk with other plant hormones including auxin [48][49][50]. For example, in vascular differentiation, cytokine affects the orientation of PIN proteins in cell membrane while auxin regulates cytokinin signaling [51,52]. Besides, also during lateral root organogenesis, cytokinin enhances the PIN1 depletion from cell membrane to affect PIN1 polarization [53][54][55]. Furthermore, it is reported that the localization of PIN proteins is affected by diffusible molecules such as jasmonate and narciclasine [56,57]. It is not yet clear whether these molecules are involved in the phyllotaxis pattern or not. However, in near future, we hope that a molecule predicted theoretically in this study will be revealed experimentally.  Fig 5 (Fig 1F). Numerical simulations were carried out in a similar manner as shown in Fig 6. Equations and regulatory functions are used as in S1 Table with