As permutations have been notably reported in Arabidopsis (Besnard et al., 2014; Guédon et al., 2013; Landrein et al., 2015; Refahi et al., 2011) and in sunflower (Couder, 1998), we sampled a variety of unrelated species in the wild and searched for permutations. We could easily find permutations in several other Brassicaceae showing spiral phyllotaxis as well as in either monocotyledonous or dicotyledonous species from more distant families such as Asparagaceae, Sapindaceae or Araliaceae (Figure 2) (Appendix section 1). As suggested by the results on Arabidopsis, these observations raise the possibility that these organ permutations result from a noise on the plastochron and that such perturbation could be a common feature of phyllotactic systems that occurs in meristems with different geometries. In addition this disorder is probably under complex genetic control as different unrelated genetic modifications (including the Arabidopsis ahp6 mutant) can modulate its intensity.

Figure 2

Download asset Open asset

To understand how the timing of organ initiation could be affected during meristem growth, we first analyzed the relative stability of inhibitory field minima in the classical deterministic model. For this, we implemented a computational version of the classical model based on (Douady and Couder, 1996b) (Appendix section 2). Primordia are created on the meristem surface at the periphery of the meristem central zone, at a fixed distance ${\displaystyle R}$ from the meristem center. Once created, primordia drift away radially from the central zone at a speed proportional to their distance from the meristem center. As soon as created, a primordium ${\displaystyle q}$ of radius ${\displaystyle r_0}$, generates an inhibitory field, ${\displaystyle E^{(q)}}$ in its neighborhood such that at a point ${\displaystyle x}$, at a distance ${\displaystyle d(q,x)>r_0}$, the inhibition due to ${\displaystyle q}$ decreases with the distance to ${\displaystyle q:E^{(q)}(x)=(\frac{r_0}{d(q,x)}{)}^s}$, where ${\displaystyle s}$ is a geometric stiffness parameter (${\displaystyle E^{(q)}}$ is often regarded as an inhibitory energy emitted by primordium q, e.g. Douady and Couder, 1996b). As a result, at any moment and for any point ${\displaystyle x}$ of the meristem surface, the existing primordia create altogether a cumulated inhibition ${\displaystyle E(x)}$ that is the sum of the individual primordium contributions: ${\displaystyle E(x)=\sum _{q\in Q}{E}^{(q)}(x)}$, where ${\displaystyle Q}$ denotes the set of all preexisting primordia. High inhibition levels on the peripheral circle prevent the initiation of new primordia at corresponding locations. However, as the preexisting primordia are moving away from the central zone during growth, the inhibitory level tends to decrease at each point of the peripheral circle. A new primordium initiates where and when the inhibitory field is under a predefined threshold on the peripheral circle.

We then performed a systematic analysis of the inhibition profiles and their dynamics along the peripheral circle and observed the following properties (Figure 3, Video 1,2):

• Property 1: let ${\displaystyle i}$ and ${\displaystyle j}$, ${\displaystyle i<j}$, denote the numbers of contact spirals that one can observe on a stem in clockwise and counter-clockwise directions (Appendix section 1), the number of local minima ${\displaystyle l}$ in the inhibition profile is bounded by ${\displaystyle i}$ and ${\displaystyle j:i\le l\le j}$: (Figure 3A–B).

• Property 2: The angular distance between local minima is a multiple of the divergence angle (Figure 3A–B).

• Property 3: The difference between the inhibition values of local minima decreases monotonically with the control parameter ${\displaystyle \mathrm{\Gamma }=r_0/R}$ (Figure 3C)

In addition, we noticed that depending on time, the difference in inhibition level between consecutive local minima may markedly vary. In some cases this difference is so small that a biological noise may lead the biological system to perceive the ordering between two or more local minima differently from the ordering of the actual inhibition levels. Such errors would lead to initiate several primordia together or to change the temporal order of their initiation, thus inducing perturbations in the sequence of divergence angle. Also, as suggested by property 3, the number of primordia initiation events affected by these errors would decrease with the ${\displaystyle \mathrm{\Gamma }}$ parameter and would thus depend on the geometry of the meristem.

Figure 3 with 1 supplement see all

Download asset Open asset

Video 1

Download asset

Video 2

Download asset

To investigate this possibility, we started to induce a perturbation in a stationary spiral pattern by forcing at a given time the system to initiate a primordium at the site of second local minimum instead of that of the global minimum (i.e. at an angle ${\displaystyle 2\varphi }$, Figure 3G, red arrow). The system was then left free to self-organize. We observed that the next primordium was always initiated at the site of the original global minimum, resulting in a divergence angle ${\displaystyle -\varphi }$ and a quasi-null plastochron (Figure 3H). The system was then able to recover from the perturbation by initiating the next primordium at the originally expected site (i.e. with a divergence angle ${\displaystyle 2\varphi }$), with a long plastochron, leading to a M-shaped pattern (Figure 3G, Video 3). The rest of the divergence angle sequence was then not affected and remained at a value close to ${\displaystyle \varphi }$ while long oscillatory perturbations were observed on plastochrons (Figure 3H). Stronger perturbations induced by forcing various other local minima to initiate instead of the global minimum (Figure 3—figure supplement 1 and Appendix section 6.6 for related control simulation experiments) similarly demonstrated that the system spontaneously makes a short distorted pattern and then returns to the normal ${\displaystyle \varphi }$ baseline in every case. We concluded that i) a noise in the perception of the local primordium order does not propagate far in the divergence angle sequence and that angle specification in the classical model patterning system is highly robust to perturbations in local minima initiation ordering and ii) however, the plastochron itself is affected during a much longer time span.

Video 3

Download asset

Together, these results suggest that time and space in primordium initiation are largely decoupled in inhibitory field-driven self-organization as locations of primordia are strongly pre-specified and relatively stable, while plastochrons are not. This observation is in line with the observations from live-imaging of Arabidopsis shoot apical meristem that demonstrated that despite variability in the plastochron (with almost 30% of organs co-initiated) the specification of initiation sites was extremely robust (Besnard et al., 2014).

Our previous results suggest that high variability in the timing of organ initiation could result from the joint effect of noise in the perception of inhibitory fields and of the decoupling between space and time in this self-organizing system. We therefore decided to revisit the inhibitory field models to integrate locality and stochasticity as central components in the patterning system. For this, we kept from the classical deterministic model the assumptions related to the movement of primordia at the meristem through growth and to the definition of inhibitory fields. However, we completely reformulated the way primordia are initiated as local stochastic processes.

At any time ${\displaystyle t}$, the ${\displaystyle K}$ cells that make up the periphery of the central zone potentially may take the identity of a primordium, depending on the local value of the inhibitory field (signaling) in each cell. We assume that this switch in cell identity does not depend on a threshold effect as in the classical model but, rather, depends on the cellular perception of the inhibitory signal which, in essence, is stochastic (Eldar and Elowitz, 2010). We thus assumed that each cell ${\displaystyle k}$ reads out the local inhibitory field value, ${\displaystyle E_k(t)}$, and switches its state to primordium identity with a probability that (i) depends on the level of inhibition ${\displaystyle E_k(t)}$; (ii) is proportional to the amount of time δt the cell is exposed to this level of signaling (Gillespie, 1976; 1977):

where ${\displaystyle X_k(t,\delta t)}$ denotes the number of primordia initiated at cell ${\displaystyle k}$ in the time interval ${\displaystyle [t,t+\delta t]}$ (${\displaystyle X_k(t,\delta t)}$ will typically have a value 0 or 1), and ${\displaystyle \lambda }$ is a rate parameter that depends on the local inhibition value ${\displaystyle E_k(t)}$ at cell ${\displaystyle k}$ and that can be interpreted as a temporal density of initiation. To express the influence of inhibitory fields on the probability of initiation, the dependence of ${\displaystyle \lambda }$ on the local inhibition ${\displaystyle E_k(t)}$ must respect a number of general constraints: i) ${\displaystyle \lambda }$ must be a decreasing function of the inhibition as the higher the inhibition level at one site, the lower the probability to observe an initiation at this site during ${\displaystyle \delta t}$; ii) for small ${\displaystyle \delta t}$, for ${\displaystyle P(X_k(t,\delta t)=1)}$ to be a probability, we must have ${\displaystyle 0\le \lambda (E_k(t))\delta t\le 1}$; iii) the ratio of the probabilities to trigger an initiation at two different sites is a function of the difference in inhibition levels between these sites. Under all these constraints, the rate parameter takes the following form (see Appendix section 3 for a detailed derivation of the model):

where ${\displaystyle E^{\ast }}$ is a parameter controlling the sensitivity of the system to inhibition, and ${\displaystyle \beta }$ is a parameter controlling the ability of the system to discriminate between inhibition levels, i.e. to respond differently to close inhibition levels (acuity). Therefore, for each cell ${\displaystyle k}$, the probability to initiate a primordium during a small time interval ${\displaystyle \delta t}$ can be expressed as:

If we now assume that the probabilities to observe an initiation at a site ${\displaystyle k}$ in disjoint time intervals are independent, the process described by Equation 3 is known as a non-homogeneous Poisson process of intensity ${\displaystyle \lambda (E_k(t))}$ (e.g. Ross, 2014). Therefore for each cell ${\displaystyle k}$ of the periphery, our model assumes that the probability to initiate a primordium is a non-homogeneous Poisson process, whose parameter is regulated by the local level of the inhibitory field at that site.

This stochastic formulation of the model at the level of cells, called SMPmicro (Stochastic Model of Phyllotaxis at microscopic level), makes it possible to develop the calculus of different key quantities or properties of the system. For example, if we assume that recruitments of cells for organ initiation are stochastically independent from each other (the probability to draw an initiation at a site ${\displaystyle k}$ only depends on the value of the local parameters at site ${\displaystyle k}$ but not on what may be drawn at other places), then we can estimate the expected number of cells independently recruited for organ initiation during the timespan ${\displaystyle \delta t}$. Let us denote ${\displaystyle X(t,\delta t)}$ the number of cells initiated along the peripheral circle, and for a given time ${\displaystyle t}$ and a small time span ${\displaystyle \delta t}$, ${\displaystyle X(t,\delta t)=\sum _{k=1}^K{X}_k(t,\delta t)}$. Its expectation is simply the sum of the expectations of the individual independent Poisson processes:

therefore giving us an estimate of the expected number of peripheral cells initiated during time ${\displaystyle \delta t}$ and for the inhibition profile ${\displaystyle E(t)}$.

So far we have considered peripheral cells as independent sites that may independently switch to primordium identity with a probability that depends on their local level of inhibition. As a local inhibition valley may span over several cells, one might expect that the probability to trigger a primordium initiation in a valley is increased by the fact that several founder cells can potentially contribute to this initiation during time ${\displaystyle \delta t}$. To formalize this, we then upscaled our stochastic model at the level of valleys where a stochastic process is now attached to each local minimum ${\displaystyle l}$ instead of to each cell ${\displaystyle k}$. This upscaled model is called SMPmacro. The idea is that the stochastic processes of all the cells ${\displaystyle k}$ spanned by a local inhibition valley ${\displaystyle l}$ sum up and together define a stochastic process ${\displaystyle N_l}$ that a primordium is initiated in ${\displaystyle l}$ at a higher level. Being the sum of independent Poisson processes, this upscaled process is also a Poisson process with intensity ${\displaystyle {\mathrm{\Lambda }}_l=\sum _{k\in K_l}{\lambda }_k(E_k(t))}$, where ${\displaystyle k}$ varies over the set of cells spanned by the valley of the ${\displaystyle l}$ th local minimum and indexed by ${\displaystyle K_l}$. If ${\displaystyle L}$ denotes the number of local minima, then the expected number of primordia initiations during a small time laps ${\displaystyle \delta t}$ is thus:

Therefore, at microscopic scale, SMPmicro couples i) a deterministic part inherited from the classical model and related to the geometry and the dynamics of the fields and ii) a new stochastic part related to the perception of this inhibitory field, i.e. representing signal perception capacities. It relies on the assumption that the perception sites, corresponding to the cells surrounding the central zone periphery, are stochastically independent of each other. Decisions regarding primordium initiation are taken in a cell-autonomous manner, thus reflecting more realistically the outcome of the initiation signaling pathway in each cell. At a macroscopic level, in each inhibition valley several cells may trigger primordium initiation. The probability to trigger an initiation increases in SMPmacro with the size of the valley when more than one founder cell are likely to contribute to initiation. A variant of this upscaled model consists of defining the probability for a valley to initiate a primordium by the probability of the cell with lowest level of inhibition in this valley. In this variant, called SMPmacro-max, ${\displaystyle {\mathrm{\Lambda }}_l=\underset{k\in K_l}{max}{\lambda }_k(E_k(t))}$.

To study the emergent properties of this system, we implemented a computational version of SMPmacro-max and tested its sensitivity to parameter changes. In addition to the geometrical parameters of the classical model, two new parameters, ${\displaystyle \beta }$ and ${\displaystyle E^{\ast }}$ now reflect the ability of the system to perceive the inhibitory signal from the fields (or equivalently the initiation signal).

As expected from a phyllotaxis model, the stochastic model is able to produce both spiral and whorl modes, from either imposed initial distributions of organs (Figure 4A–C), or random starting points (Appendix section 6.1). However, the great majority of the sequences of divergence angle generated for different values of ${\displaystyle \mathrm{\Gamma }=r_0/R}$ displayed divergence angle perturbations (i.e. angles different from those predicted by the classical model with the same ${\displaystyle \mathrm{\Gamma }}$) of the type observed in Figure 4D–E. As suggested by their typical distributions (Figure 4F–H left), these perturbations correspond to permutations very similar to those observed on real plants in previous studies (Besnard et al., 2014; Landrein et al., 2015) with the appearance of secondary modes at multiples of 137, Figure 1F–I. The amplitude of these secondary modes is correlated with the amount of perturbations in the sequences.

Figure 4

Download asset Open asset

To complete this analysis, we looked at plastochron distributions in the simulated sequences (Figure 4F,G,H right). In all cases, distributions were displaying a single mode largely spread along the ${\displaystyle x}$-axis. Interestingly, the more sequences were perturbed, the more negatively skewed the distributions, showing thus a higher occurrence of short plastochrons (Figure 4H). This is reminiscent of the observation of co-initiations in growing meristems associated with perturbed phyllotaxis (Besnard et al., 2014; Landrein et al., 2015). The stochastic model is thus able to produce perturbed sequences with realistic series of divergence angles and corresponding realistic distributions of plastochrons.

We then aimed to quantitatively assess the complexity level of permutations as a proxy for plastochron noise. For this we focused on spiral phyllotaxis modes and used two measures: the density of permuted organs respectively involved in 2 and 3-permutations ${\displaystyle {\pi }_2=2.{\sigma }_2/\mathrm{\Sigma }}$ and ${\displaystyle {\pi }_3=3.{\sigma }_3/\mathrm{\Sigma }}$ where ${\displaystyle {\sigma }_2}$ and ${\displaystyle {\sigma }_3}$ are respectively the number of 2- and 3-permutations in the sequence and ${\displaystyle \mathrm{\Sigma }}$ the total number of organs in the sequence. The quantities ${\displaystyle {\sigma }_2}$ and ${\displaystyle {\sigma }_3}$ were estimated a posteriori from simulated sequences using the algorithm described in (Refahi et al., 2011). The total density of permuted organs involved either in 2- or 3 permutations is denoted by ${\displaystyle \pi ={\pi }_2+{\pi }_3}$.

We first explored the intensity of permutations of different natures (2- and 3-permutations) in simulated phyllotaxis sequences. For this, we carried out simulations for a range of values for each parameter ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$ (Figure 5—source data 1). These values can be regarded as predictions of the model for each triplet ${\displaystyle (\mathrm{\Gamma },\beta ,E^{\ast })}$. Remarkably, the simulations show that the values of ${\displaystyle \pi ,{\pi }_2}$ and ${\displaystyle {\pi }_3}$ are not linearly related to each other. To make this relationship explicit, we plotted the proportion of organs involved in 2-permutations ${\displaystyle {\pi }_2}$ as a function of the total proportion of perturbed organs ${\displaystyle \pi }$ (Figure 5). Surprisingly, the points, when put together on a graph, were organized in a narrow crescent showing a convex curve-like relationship between ${\displaystyle \pi }$ and ${\displaystyle {\pi }_2}$, revealing a remarkable property of the stochastic model: the more perturbations there are in the simulated sequences, the higher the proportion of 3-permutations, and this independently of the model parameters. The fact that this non-linear relationship emerges from a fairly large sampling of the parameter space suggests that it can be considered as a key observable property characterizing the model’s underlying structure.

We tested this by gathering all measured values of ${\displaystyle {\pi }_2}$ and ${\displaystyle {\pi }_3}$ published in the literature for Arabidopsis (Besnard et al., 2014; Landrein et al., 2015) and plotted the corresponding points on the original graph showing the model’s simulations (Figure 5, red crosses). The measured points fall within the range of predicted values and show that the measured values follow the same non-linear variation as the one predicted by the model: the larger the total percentage of permutations, the larger the proportion of 3-permutations in the sequences. This confirms a first prediction from the stochastic model and indicates that disorder complexity increases non-linearly with the frequency of disorders.

Figure 5

Download asset Open asset

Figure 5—source data 1

Source files for simulated permutation intensities.

This file contains a table showing the variation of permutation intensities with model parameters. We ran simulations using the SMPmacro-max model for different parameter values where the local minima of inhibition profile indicate the potential initiation sites (for this, the inhibitory field values were estimated at 360 sampling points regularly distributed around the periphery of the central zone). We then used the combinatorial model (Refahi et al., 2011) to detect permutations in the simulated sequences. The model has mainly three parameters, ${\displaystyle \beta }$, ${\displaystyle E^{\ast }}$ and ${\displaystyle \mathrm{\Gamma }}$. For each parameter value of ${\displaystyle \beta }$, ${\displaystyle E^{\ast }}$ and ${\displaystyle \mathrm{\Gamma }}$, we run sixty simulations. Each simulation generated a sequence of 25 divergence angles. We then analyzed the sequences using the combinatorial model to detect permutation patterns.

https://doi.org/10.7554/eLife.14093.012

Download elife-14093-fig5-data1-v1.pdf

Based on our simulations, we then investigated the variations of these perturbations as a function of the model parameters ${\displaystyle \mathrm{\Gamma },\beta .E^{\ast }}$ and s. As a general trend, for a value of s fixed to 3 as in (Douady and Couder, 1996b), we noticed that for a given value of ${\displaystyle \mathrm{\Gamma }}$, an increase in ${\displaystyle \beta }$ was roughly counteracted in terms of permutation intensity by a decrease in ${\displaystyle E^{\ast }}$. Likewise, an increase of ${\displaystyle \mathrm{\Gamma }}$ was canceled by a decrease in ${\displaystyle \beta }$. We gathered all these observations in one graph and plotted the global amount of perturbation ${\displaystyle \pi }$ in a sequence as a function of a combination of the three original model parameters and reflecting the observed trend: ${\displaystyle {\mathrm{\Gamma }}_P=\mathrm{\Gamma }\beta E^{\ast }}$. In the resulting graph (Figure 6A), each point corresponds to a particular instance of the three model parameters. The points form a narrow and decreasing band associating a small set of possible perturbation intensities π with each value of ${\displaystyle {\mathrm{\Gamma }}_P}$. For combinations of the three model parameters leading to a small ${\displaystyle {\mathrm{\Gamma }}_P}$, the perturbations can affect up to 50% of the organs whereas for high values of ${\displaystyle {\mathrm{\Gamma }}_P}$, there may be no perturbation at all. Interestingly, ${\displaystyle {\mathrm{\Gamma }}_P}$ being a combination of the three elementary model parameters, plants having identical geometrical parameter ${\displaystyle \mathrm{\Gamma }}$ may show substantially different intensities of perturbations if their perception for different values of parameters ${\displaystyle \beta ,E^{\ast }}$. Consequently, since the intensity of perturbation in the system is more simply reflected by ${\displaystyle {\mathrm{\Gamma }}_P}$ than by values of ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$ taken independently, ${\displaystyle {\mathrm{\Gamma }}_P}$ can be considered to be a control parameter for perturbations.

Figure 6 with 1 supplement see all

Download asset Open asset

To further investigate the structure of the stochastic model, we then studied how the usual observable quantities of a phyllotaxis system, i.e. divergence angles and plastochrons, depend on the model parameters ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$ and s.

In the classical deterministic model, divergence angles are a function of a unique control parameter ${\displaystyle \mathrm{\Gamma }=r_0/R}$ (Douady and Couder, 1996b), meaning that the same divergence angle can be obtained in the model for different couples of ${\displaystyle r_0}$ and ${\displaystyle R}$ provided that their ratio is unchanged. We thus checked whether ${\displaystyle \mathrm{\Gamma }}$ could also serve as the control parameter for the divergence angles and plastochrons in the stochastic model. For this, we simulated various spiral phyllotaxis sequences (from the Fibonacci branch where divergence angles are close to 137°) by varying ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$ and estimated the corresponding divergence angle ${\displaystyle \alpha }$ and plastochron T, Figure 6B,D. We observed that, although the point clouds evoke the corresponding curves in the standard deterministic model (Appendix 1—figure 2), significantly different values of the divergence angles (or plastochrons) can be observed for many particular values of ${\displaystyle \mathrm{\Gamma }}$. This phenomenon suggests that ${\displaystyle \mathrm{\Gamma }}$ is not a satisfactory control parameter in the stochastic model: for a given value of ${\displaystyle \mathrm{\Gamma }}$, varying ${\displaystyle \beta }$ or ${\displaystyle E^{\ast }}$ can significantly modify the divergence angle or the plastochron. We therefore tested various combinations of these parameters (Figure 6—figure supplement 1), and finally found that for ${\displaystyle {\mathrm{\Gamma }}_D={\mathrm{\Gamma }}^{s/3}\frac{1}{{\beta }^{1/6}{E}^{\ast 1/2}}}$ both cloud of points collapse into a single curves (Figure 6B–C, D–E and Appendix section 6.3). For this definition of ${\displaystyle {\mathrm{\Gamma }}_D}$, each value of the control parameter can be associated with quasi-unique divergence angle and plastochron, i.e. defines a precise observable state of the system.

We then checked whether parastichies were also controlled by ${\displaystyle {\mathrm{\Gamma }}_D}$. For this, we computed the parastichies ${\displaystyle (i,j)}$ corresponding to each simulated sequence and plotted the sum, ${\displaystyle i+j}$, as a function of ${\displaystyle {\mathrm{\Gamma }}_D}$ (Figure 6F). The resulting points were arranged on a stepwise curve, where each step corresponds to a Fibonacci mode: (1,2), (2,3), (3,5), (5,8). Each value of ${\displaystyle {\mathrm{\Gamma }}_D}$ thus defines a precise value of the mode.

We concluded that ${\displaystyle {\mathrm{\Gamma }}_D}$ can be considered as a second control parameter of the stochastic model, relating the system’s state to the observable variables α and T through a unique combination of the system’s parameters.

According to the stochastic model, a particular phyllotaxis system is characterized by a particular set of values of the parameters ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$. Upon growth, a specific spatio-temporal dynamics emerges that is characterized by observable variables: divergence angle, plastochron and frequency of permutations. In the current state of our knowledge and measuring means, the parameters ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$ are not directly observable. Therefore, we investigated what can be learned about them from the observable variables.

For this, we can use the relationships established above between the control parameters ${\displaystyle {\mathrm{\Gamma }}_P}$ and ${\displaystyle {\mathrm{\Gamma }}_D}$ and the observable variables ${\displaystyle \alpha ,T,{\pi }_2,{\pi }_3,}$ .... For a genotype ${\displaystyle G}$, we have:

Using the characteristic curves of Figure 6C,E, both measurements of ${\displaystyle \alpha }$ and ${\displaystyle T}$ give possible estimates for ${\displaystyle {\mathrm{\Gamma }}_D}$. In a consistent model, these estimates should be compatible. Similarly, an estimate of ${\displaystyle {\mathrm{\Gamma }}_P}$ can be derived from the observation of ${\displaystyle \pi }$. In the case of the WT for example, ${\displaystyle {\mathrm{\Gamma }}_D}$ is estimated to be in [0.425, 0.5] and ${\displaystyle {\mathrm{\Gamma }}_P}$ in [8.30, 11.49] according to the observed value of the plastochron, divergence angle and permutation intensity. This value is itself derived from the measurement of the plastochron ratio ${\displaystyle \rho }$ (Appendix section 4.1.3).

Based on estimated values of ${\displaystyle {\mathrm{\Gamma }}_D}$ and ${\displaystyle {\mathrm{\Gamma }}_P}$, equations 6 define a system of 2 equations and 3 unknowns (s was fixed to 3 in the simulations). This system is underdetermined and does not allow us to identify exactly the values of ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$. However, being underdetermined in one dimension only, this system plays the role of a generalized control parameter: a given value of the generalized control parameter (${\displaystyle {\mathrm{\Gamma }}_D}$, ${\displaystyle {\mathrm{\Gamma }}_P}$) determines the parameters ${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$ up to one degree of freedom. If one of the parameters is given, then the others are automatically determined according to equation 6.

Several recent works demonstrated that mutations or changes in growth conditions could alter phyllotaxis and disorder patterns. Our model predicts correlations between main observable phyllotaxis variables. According to the model, plastochron positively correlates with ${\displaystyle {\mathrm{\Gamma }}_D}$ (Figure 6E). If signaling is not altered by the experimental setup (${\displaystyle \beta }$ and ${\displaystyle E^{\ast }}$ are unchanged) then ${\displaystyle {\mathrm{\Gamma }}_D}$ positively correlates with ${\displaystyle \mathrm{\Gamma }=r_0/R}$ (Equation 6). Therefore, the model predicts that, like in the classical model, plastochron positively correlates with r₀ (size of the primordia inhibitory fields) and negatively correlates with R (size of the central zone). However, the model makes also in this case the new prediction that plastochron negatively correlates with the frequency of observed permutations (Figure 6A).

A series of recent observations support this prediction. By changing growth conditions (plants first grown in different day-length conditions and then in identical conditions) or by using different accessions or mutants with markedly different meristem sizes from that of the wild type (Landrein et al., 2015), changes in the size of the meristem could be induced. The authors hypothesized that this change affected the size of the central zone only and not the size of primordia inhibitory fields. Corresponding changes in Γ were observed to positively correlate with the frequency of organ permutations and negatively correlate with plastochron, as predicted by our model. In addition, the stochastic model makes it possible to quantitatively estimate the changes in central zone sizes from the measured phyllotaxis disorders with an error less than 5% (Appendix section 4.3).

In a previous study on ahp6 mutants, Besnard et al. (2014) showed that the frequency of disorders could markedly augment while the size of meristems did not significantly change like in (Landrein et al., 2015). As discussed above, this change in disorder intensity could theoretically be due to an alteration of initiation perception in the mutant. However, the stochastic model suggests that it is not the case. Indeed, we re-analyzed plastochrons of the mutants and could observe that, although the change is limited, mutant plastochrons are significantly smaller than those of the wild type (Appendix section 4.2). According to Figure 6E, this means that ${\displaystyle {\mathrm{\Gamma }}_D}$ is reduced in the mutant. If this reduction was due to an increase in either ${\displaystyle \beta }$ or ${\displaystyle E^{\ast }}$, then according to the model, one should expect a corresponding increase of ${\displaystyle {\mathrm{\Gamma }}_P}$ (Equation 6) and, thus, a decrease of disorders (Figure 6A). On the contrary a significant increase of disorders was actually observed, suggesting that perception is not altered and that, rather, a decrease in ${\displaystyle \mathrm{\Gamma }=r_0/R}$ could be the source of variation. Since the size R of the meristems did not change, the model suggests that ahp6 is altered in the size ${\displaystyle r_0}$ of the primordia inhibitory fields.

In both the previous analysis and in related works (Besnard et al., 2014; Guédon et al., 2013; Landrein et al., 2015; Refahi et al., 2011), the permutation detection was restricted to 2- and 3-permutations. However, the stochastic model potentially predicts the existence of higher order permutations, i.e. 4- and 5-permutations, in Arabidopsis thaliana especially for small values of the control parameter ${\displaystyle {\mathrm{\Gamma }}_P}$. Following this prediction, we revisited the measured divergence angles in (Landrein et al., 2015) on the WS4 mutant grown in standard conditions (${\displaystyle {\pi }_2}$= 29.45%, ${\displaystyle {\pi }_3}$= 14.91%) and for which 7.7% of angles were left unexplained when seeking for 2- and 3-permutations (Figure 7A). When higher order permutation are allowed in the detection algorithm (Appendix sections 4.1 and 5), most of the unexplained angles for WS4 can be interpreted as being part of 4- and even 5-permutations (Figure 7B).

Figure 7

Download asset Open asset

Previous observations (Landrein et al., 2015) point to the existence of positive correlations between meristem size and intensities of perturbations. As discussed above, the stochastic model explains this: if a mutation, or a change in growth conditions only affects the geometry of the system ${\displaystyle \mathrm{\Gamma }=r_0/R}$, then ${\displaystyle {\mathrm{\Gamma }}_D}$ and ${\displaystyle {\mathrm{\Gamma }}_P}$ are both affected in the same sense (Equation 6). However, it also allows predicting new observable facts.

Assume that a mutation or a change in growth conditions affects the ability of the plant to perceive initiation signals without modifying the geometry of the system, i.e. ${\displaystyle \beta }$ and/or ${\displaystyle E^{\ast }}$ are modified while ${\displaystyle \mathrm{\Gamma }}$ is left unchanged. Then, according to equation 6, this induces opposite variations in ${\displaystyle {\mathrm{\Gamma }}_D}$ and ${\displaystyle {\mathrm{\Gamma }}_P}$ that can be detected with the observed variables. For example, assuming that a mutation decreases the sensitivity of the system to the initiation signal (decrease of ${\displaystyle E^{\ast }}$), ${\displaystyle {\mathrm{\Gamma }}_P}$ is decreased while ${\displaystyle {\mathrm{\Gamma }}_D}$ is increased. The decrease of ${\displaystyle {\mathrm{\Gamma }}_P}$ induces an increase of the perturbation intensity ${\displaystyle \pi }$ while the increase of ${\displaystyle {\mathrm{\Gamma }}_D}$ induces an increase of the plastochron. The model thus predicts that, in such a case, it is possible to expect an augmentation in the disorder correlated with a decrease in organ initiation frequency. To date, such a fact has not yet been observed and constitutes a testable prediction of our model.

While the stochastic model is able to reproduce the major spiral and whorl phyllotaxis patterns, stochasticity induces alterations in the patterning process mainly affecting the plastochron i.e. the timing of organ initiation. These alterations take the form of permutations of the order of organ initiation in the meristem. If the plastochron is small as frequently observed in real sequences of permuted organs, permutations in the model can be considered equivalent to co-initiations that have been identified in the Arabidopsis meristem as the main source of permutations observed on the inflorescence stem. Less frequently, simulated permutations can have longer plastochrons. In this case, they can be interpreted as true permutations of the order of organ initiation in the meristem, consistently with the low frequency of such meristematic permutations observed also in Arabidopsis (Besnard et al., 2014). These results are in line with a previous attempt at introducing stochasticity in the classical deterministic model that could also induce permutations (Mirabet et al., 2012). However in this latter work only a limited frequency of defects could be induced (even when the noise was fixed at high levels), while requiring time discretization and post-meristematic randomization of organ order when more than one organ initiation were detected in the same simulation time step. By contrast, the capacity of the stochastic model to reproduce faithfully perturbed sequences as observed in Arabidopsis indicates that the model captures accurately the dynamics of the phyllotactic system. Taken with the fact that permutations are observed in a variety of species and genotypes from a given species, our theoretical results identify noise on the plastochron as a common characteristic of phyllotaxis systems that may generate disorders in these developmental systems.

It is important to note that our work points at stochasticity in signaling mechanisms allowing perception of inhibitory fields as the most likely origin of this developmental noise but does not entirely rule out the idea that other phenomenon might contribute (see Appendix section 6). A major contribution of stochasticity in signaling is supported by the robustness of the model to changes in different assumptions and parameters. However we have observed that spatial discretization for example can modify the frequency of permutations in the model, although this effect is limited (see Appendix section 6.2). A possible interpretation is that changes in the size of cells could also contribute to a certain point to noise on the plastochron, an idea that could be further explored.

Importantly, phyllotaxis dynamics in the stochastic model relies not only on the geometry of the inhibitory fields captured by the ${\displaystyle \mathrm{\Gamma }}$ parameter as in previous deterministic models (Douady and Couder, 1996b; Richards, 1951) (Figure 8B), but also on two new parameters ${\displaystyle E^{\ast }}$ and ${\displaystyle \beta }$ (Figure 8C). ${\displaystyle E^{\ast }}$ and ${\displaystyle \beta }$ describe respectively the sensitivity of cells to the inhibitory signal and the acuity of their perception, i.e. their capacity to differentiate close signal values. Our work thus suggests that a robust self-organization of a 3D developmental system driven by lateral inhibition depends both on the geometr of inhibitory fields in a tissue but also on the signaling capacities of cells in tissues. These theoretical observations are consistent with the key role of signaling in phyllotaxis (Vernoux et al., 2011), and with the setting of patterning dynamics in animal systems, downstream of morphogenetic signals (Kutejova et al., 2009). This pinpoints the interplay between global information provided by signal distribution and local interpretation of the information as a general principle for patterning emergence. In addition, we predict that due to pre-specification of initiation sites, noise in phyllotaxis is expected mainly on the timing of patterning. This might explain the selection through evolution of genetic mechanisms, such as the one recently described implicating the AHP6 protein (Besnard et al., 2014), to diminish noise on plastochron and disorders in phyllotaxis.

In biological systems, disorders are frequently viewed as a result of biological or environmental noise that mainly alters systems function or development. It is in this sense for instance that noise on phyllotaxis patterns had previously been analyzed (Itoh et al., 2000; Jeune and Barabé, 2006; Peaucelle et al., 2008). Here we show that biological noise at microscopic scale may be revealed at macroscopic scale in the form of organ disorders, the permutations. Our stochastic model of phyllotaxis suggests that these disorders bear information on the more profound, hidden variables that control the phyllotaxis patterning. Much like digital watermarks that represent a copyright or any information to be hidden in images or audio signals, the actual variables (${\displaystyle \mathrm{\Gamma },\beta ,E^{\ast }}$) representing the state of the phyllotaxis system are not directly apparent in the plant phenotype (i.e. in the sequence of lateral organ angles and its dynamics). However, by scrutinizing carefully the image or, here, the phyllotaxis pattern and their perturbations with adequate decoding algorithms, it is possible to reveal the hidden information that was “watermarked” in the original signal. In this way, permutations together with divergence angle α and plastochron ${\displaystyle T}$ reveal key information about the state of the system that has produced them, as their knowledge drastically reduces the set of possible ${\displaystyle \mathrm{\Gamma }}$, ${\displaystyle \beta }$ and ${\displaystyle E^{\ast }}$ values. Reciprocally, any experimental alteration of these values modifies the biological watermark. Our model suggests that this change is reflected in macroscopic alterations of the phyllotaxis patterns that convey information about their possible molecular origin.

To illustrate this, we used our stochastic model to confirm that changes in permutation frequencies due to changes in growth conditions are most likely explained by a specific modulation of ${\displaystyle \mathrm{\Gamma }}$, as previously proposed (Landrein et al., 2015) (Appendix section 4.1.5). Such a biological watermarking also suggests a different interpretation of the function of AHP6 in phyllotaxis (Besnard et al., 2014). Movement of AHP6 from organs has been proposed to generate secondary inhibitory fields that filter co-initiation at the meristem, decreasing the frequency of permutations. Based on permutation modifications, our theoretical framework suggests that AHP6 effect on the phyllotactic system does not need to be viewed as an additional specific mechanism acting on plastochron robustness but could be simply interpreted as a mechanism increasing ${\displaystyle \mathrm{\Gamma }}$ slightly. As no differences in auxin-based inhibitory fields could be detected between ahp6 mutant and wild-type plants (Besnard et al., 2014) (Appendix section 4.2), our model thus leads to a vision with composite inhibitory fields resulting at least from the combined effect of auxin-based and AHP6-based subfields. Conversely, this predicts that inhibitory fields in general cannot be explained only by auxin-based mechanisms as previously proposed.

Combined with data on divergence angles and on the plastochron, we further predict that the phyllotactic disorders could be used to identify mutants affected in biological mechanisms that controls β and/or E*. The model indicates that such mutations would have an opposite effect on frequency of permutations and on plastochron. Mutants behaving accordingly would allow not only to test this prediction of the model but also to dissect the molecular mechanisms at work. Precise and automated quantifications of the permutation, divergence angles and plastochron would allow for screening for such mutants and should become feasible with the fast development of phenotyping tools (Dhondt et al., 2013; Granier and Vile, 2014).

Our model only takes into account stochasticity in the perception of inhibitory fields by cells and is based on two biologically plausible assumptions: that this perception is mostly cell autonomous and that it only depends on the local level of the inhibitory signal. This provides a reasonable abstraction of local stochastic fluctuations in i) hormonal concentrations related to inhibition produced by each primordium ii) in the activity of the signal transduction pathway leading to initium creation. The detailed molecular mechanisms controlling organ initiation are for the moment only partially known. However the capacity of the stochastic model to capture accurately phyllotaxis suggests that it also captures plausible emergent properties of the underlying molecular mechanisms. This model thus not only provides a framework to understand the dynamics of patterning in the meristem but also the properties of the signaling mechanisms that process the different signals involved. Note also that the predictive capacities of our model suggest that noise on perception could be the most influential source of noise in the system. However demonstrating this would require further exploration of other potential sources of stochasticity acting at different scales, such as growth variations, spatial discretization of the peripheral zone (to account for the real size of plant cells), in order to assess their relative contribution to disorders. Moreover, similarly to the classical deterministic model of phyllotaxis, our stochastic model does not explicitly account for the cascade of molecular processes that participate to the establishment of new inhibitory fields at the location of incipient primordia. This might limit the ability of these models to fully capture the dynamics of the self-organization of the system. To do so, more mechanistic versions of this stochastic model could be developed in the future, combining more detailed cellular models of hormone-based fields, e.g. (Jönsson et al., 2006; Smith et al., 2006a; Stoma et al., 2008), and stochastic perception of these hormonal signals in 2D or 3D models with cell resolution.

Heterogeneity of biological systems at all scale has attracted an ever-growing attention in the recent years (Oates, 2011). Deterministic models do not account for the high variability that can be observed in systems behaviors, indicating that they fail to capture some key characteristics of biological systems (Wilkinson, 2009). While more demanding computationally, stochastic models are required in such cases, e.g. (Greese et al., 2014; Uyttewaal et al., 2012; Wennekamp et al., 2013), and our work illustrates how dynamic stochastic modeling can help understanding quantitatively self-organization and more broadly patterning in higher eukaryotes.

Based on the classical model of phyllotaxis (Appendix section 2), a complete and formal presentation of the stochastic model is described in the Appendix section 3. In particular, it is shown how the exponential form of the intensity law can be derived from basic model assumptions and how different observable quantities can be expressed using the model parameters.

A computational version of the stochastic model SMPmacro-max was implemented in Python programming language using Numpy and SciPy . Similarly to (Douady and Couder, 1996b), unless otherwise stated, the stiffness parameter was fixed in all simulations to s = 3. The non-homogeneous Poisson process was simulated using the algorithm described in (Ross, 2012). A pseudo-code version of the stochastic model algorithm is given in the Appendix section 3.2.

To estimate the value of the different variables characterizing phyllotaxis ${\displaystyle \alpha }$, ${\displaystyle T}$ and ${\displaystyle \pi }$, ${\displaystyle {\pi }_2}$, ${\displaystyle {\pi }_3}$, etc. in either simulated or observed sequences, we used a method based on the algorithm developed in (Refahi et al., 2011) and described in the Appendix section 4. As this algorithm is central to the identification of permutations, we additionally tested its ability to detect correctly permutations on synthetic data in which known permutation patterns were introduced (Appendix section 5). Results show that the algorithm is able to detect permutations with a success rate of 98% on average.

The parameter space of the stochastic model was explored by varying values of ${\displaystyle \mathrm{\Gamma }}$, ${\displaystyle \beta }$ and ${\displaystyle E^{\ast }}$. 60 stochastic runs have been made for each 3-tuple of the parameter values. Each simulation run generated a sequence of 25 divergence angles and corresponding plastochrons. The different observable variables have been extracted from these simulations. Results are reported in Tables 1 and 2 of the Figure 5—source data 1.

The models describing the different non-linear relationships between the observable variables and the control parameters ${\displaystyle {\mathrm{\Gamma }}_D}$ and ${\displaystyle {\mathrm{\Gamma }}_P}$ were fitted with Gauss-Newton non-linear least-squares method (Bates and Watts, 2007). Approximate 95%-prediction bands of the response variables were computed by assuming random errors of the models independent and identically normally distributed.

Phyllotactic patterns emerge from the arrangement of organs on concentric circles, each circle corresponding to a ”wave” of organ (or primordia) production, which is driven centrifugally by growth. Each circle contains the same number, called jugacy, of regularly spaced primordia, and is rotated with respect to the previous circle by a constant angle, called the divergence angle. Denoting the divergence angle by $\varphi$ and the jugacy by $j$, phyllotactic patterns are traditionally classified as follows:

Note that to completely characterize a phyllotactic pattern, one also needs to specify the radii of the concentric ”co-initiation” circles mentioned above. A first parameter needed for this is the radius of the central zone, a circular region at the centre of the apical meristem where no organ is ever produced. This radius is thus a lower bound for the radii of the co-initiation circles. Some regularity has been observed (or assumed) in these radii, and they often are specified thanks to a single additional parameter, which can be either

• the plastochron, defined as the time between two successive waves of organ formations, or

• the plastochron ratio, defined as the ratio between the radii of two successive co-initiation circles.

The use of these quantities as characteristic parameters implies that they are constant, which is only partially corroborated by the observation. In particular, permutations of organs correspond to irregularities of the plastochron. They allow to complete the description of ”ideal” phyllotactic patterns such as those depicted in Appendix 1—figure 1.

Appendix 1—figure 1

Download asset Open asset

In the spiral case (which is the most widespread in nature), as the name indicates successive organs are arranged along a generative spiral, see Appendix 1—figure 1a. In multijugate cases, $j$ copies of this spiral are superimposed, each rotated by an angle multiple of $\pi /j$. When a sufficient number of organs are present, other spirals naturally occur to the eye. These spirals are formed by series of nearest neighbours, or second nearest neighbours. Such spirals are called parastichies. Because they are more apparent than the generative spiral, these parastichies are often used as characteristics of a phyllotactic pattern. More precisely, the parastichies formed by sequences of nearest and second nearest neighbours appear as spiralling outwards in opposite directions. One can thus count the number of spirals oriented in each direction, which determines a pair of integers $(p,q)$ called the mode of the phyllotactic pattern. A well-known observation is that the divergence angle is often close to the golden mean angle (${\displaystyle \approx {137.51}^{\circ }}$), which leads to phyllotactic modes $(p,q)$ where $p$ and $q$ are successive numbers from the Fibonacci sequence ${(F_n)}_{n\ge 1}=(1,1,2,3,5,8,13,21,\mathrm{\dots })$ is defined by $F_1=F_2=1$ and $F_{n+2}=F_n+F_{n+1}$. Other irrational numbers are also observed, giving rise to sequences closely related to the Fibonacci sequence. The most frequent pattern after Fibonacci is probably the so-called Lucas phyllotaxis, where ${\displaystyle \varphi \approx {99.5}^{\circ }}$ and modes are successive terms in the Lucas sequence ${(L_n)}_{n\ge 1}=(2,1,3,4,7,11,18,29,\mathrm{\dots })$ is defined by $L_1=2,L_2=1$, and $L_{n+2}=L_n+L_{n+1}$.

To generate phyllotactic patterns two main mechanisms have been proposed. In 1868, Hofmeister proposed that a new primordium in generated periodically (the period being the plastochron defined above) in the largest available space on the central zone periphery. Later in 1962, Snow & Snow suggested that a new primordium is generated on the central zone periphery where and when there is available space, given the previous history of the system. In this view, phyllotaxis (and in particular the plastochron) emerges as a consequence of a local production rule.

Before describing models in more details, let us simply incorporate the notion of permutations in the description of phyllotaxis above. The permutations observed in plants and discussed in the text concern the order of organs along the generative spiral (hence they are defined for spiral phyllotactic patterns). We call $n$-permutation a series of $n$ organs whose order along the generative spiral is permuted. The order on the spiral corresponds to the order of appearance of primordia, and is materialized geometrically by the radial position which reflects the passage of time (see next section for more details). Hence, an $n$-permutation corresponds to a situation where $n$ successive organs as shown on Appendix 1—figure 1a have non-increasing radial positions: for instance, there is a $3$ permutation whenever the radial coordinates of organs $3$, $4$ and $5$ verify ${\displaystyle r_3>r_4}$ and ${\displaystyle r_3>r_5}$ (with the obvious notation). A more formal definition of this notion of $n$-permutations was published in (Refahi et al., 2011).