Emergence of complex patterns in a higher-dimensional phyllotactic system

A hypothesis commonly known as Hofmeister’s rule states that primordia appearing at the apical ring of a plant shoot in periodic time steps are formed in the position where the most space is available with respect to the space occupation of already-formed primordia. A corresponding two-dimensional dynamical model has been extensively studied by Douady and Couder, and shown to generate a variety of observable phyllotactic patterns indeed. In this study, motivated by mathematical interest in a theoretical phyllotaxis-inspired system rather than by a concrete biological problem, we generalize this model to three dimensions and present the dynamics observed in simulations, thereby illustrating the range of complex structures that phyllotactic mechanisms can give rise to. The patterns feature unexpected additional properties compared to the two-dimensional case, such as periodicity and chaotic behavior of the divergence angle.


Introduction
In the history of phyllotaxis and its mathematical analysis, a hypothesis commonly attributed to botanist Wilhelm Hofmeister [1] on the dynamical formation of primordia at the shoot's apical ring has received particular attention.In the course of the shoot's growth, existing primordia radially move away from the tip of the axial-symmetric apex.The hypothesis considers the scenario of a regular time interval between the formation of two successive primordia, the plastochrone, and states that the formation of the newest primordium on the apical ring occurs in the largest available space left by the existing primordia.Kirchoff [2] pointed out that this common presentation of Hofmeister's rule in the contemporary literature, albeit biologically reasonable, does not exactly correspond to Hofmeister's original work.Here, we use the term "modified Hofmeister rule" in order to connect to related works while accommodating historical fact.
Following the early geometric work by van Iterson [3], the notion of space availability is often understood as the space between hard discs representing already-formed primordia.Atela et al. [4] presented a rigorous mathematical analysis of a corresponding dynamical model based on the modified Hofmeister rule, in which the apex is represented as a cylinder.A related approach is based on the assumption that each already-formed primordium generates a repulsive energy that decreases with increasing distance from it.The location offering the largest available space is then defined as the position in which the cumulative energy potential generated by all existing primordia is minimal.In particular, this accounts for the influence exerted by primordia other than the immediate adjacents on the position of formation of the newest.Several authors [5][6][7][8][9] (and references therein) exemplified this idea.Notably, Douady and Couder [10] applied it to the modified Hofmeister rule and presented a seminal study on quantitative and qualitative dynamical properties of their model.
Along extensive research focusing on the analysis of static phyllotactic patterns, hypotheses on the dynamics of their formation alternative to the modified Hofmeister rule have been put forward, notably by Schwendener [11] and Snow and Snow [12], and subsequently incorporated into mathematical models (e.g., [10,13]).We refer to Adler et al. [14] for a comprehensive review.
The following model can be seen as a minimalist version of the one by Douady and Couder [10].Let the two-dimensional apex surface be represented by ℝ 2 , whose origin (0,0) marks the apex center.The unit circle S 1 represents the apical ring on which new primordia are formed before they continuously and radially move away from the origin with velocity equal to 1. Existing primordia are characterized in terms of the time t and the position x ∈ S 1 of their appearance.Each of them accounts for an energy potential of E(d) in a distance of d length units away from it.As in [10], we use E(d) = d −p for p > 0. For n existing primordia, characterized by (t i ,x i ) i=1,…,n , and a given plastochrone λ > 0, the n+1st primordium is thus formed at time and position thus, among all points on S 1 , x n+1 is chosen as the one where the combined energy potential from all previously formed primordia is minimal.This model differs from the one by Douady and Couder [10] essentially only in the geometry of the apex, which is conic in their model, and the velocity of the primordia.On a qualitative level, the dynamics are nevertheless identical.
The smaller the value p, the more influential older primordia are for the position of the newest.Conversely, as p → ∞, the dynamics converge to those in the case of hard discs [15].
We recapitulate some of the experimental results obtained by Douady and Couder [10] that hold for the model (1) as well: for any plastochrone λ and well behaved initial conditions, the divergence angle δ n = ∠x n ,x n−1 , confined by two successive primordia, converges to a stable equilibrium δ = δ(λ) as n becomes large.In the course of many plants' ontogeny, the plastochrone decreases continuously [16][17][18], motivating to initiate the model iteration at large λ and then to quasi-statically decrease it during the iteration process.The latter notion, introduced by Douady and Couder [10], refers to the iteration being run until an equilibrium state is attained, before decreasing the parameter to a minimally lower value and repeating the procedure.One observes a continuous dependency of the stable equilibrium divergence angle on λ, forming the Fibonacci branch (dotted black line in Fig. 2a).Moreover, as λ approaches 0, the appropriate divergence angle tends oscillatory to the golden angle φ = (1 + √5)/2 ≈ 137.508°,where φ = 2π/(1 + φ) denotes the golden ratio -a quantity that, among others, has been analytically related to the optimization of packing density [19] and light interception by leaves [13].
When, for any stable configuration (i.e., δ n = const.),each primordium is connected to the nearest two among those formed after it, two sets of spirals in ℝ 2 , the parastichies, emerge, whose members turn clockwise and anticlockwise, respectively.The numbers of spirals of the two sets are always two successive elements of the Fibonacci sequence, growing beyond all bounds as λ approaches 0.
In this paper, we generalize the model (1) to three spatial dimensions and study the emerging patterns by means of numerical simulations (given that the system is too complex for a rigorous mathematical analysis), in particular, with regard and in comparison to the above-mentioned planar features.Motivated by mathematical interest rather than a concrete biological problem, this aims to explore the range of structures generated by a theoretical system that is inspired by phyllotactic mechanisms.

Material and methods
The model (1) considered hitherto extends straightforwardly from ℝ 2 to higher dimensional spaces ℝ m : instead of the apical ring represented by S 1 , primordia are then formed on the sphere S m−1 = {x ∈ ℝ m : ||x|| = 1}.The active sphere S m−1 thus takes on the role of the active ring S 1 sensu Smith et al. [9] in two dimensions.Analogous to the way they move away from (0,0) in the previous two-dimensional case, they subsequently move away from the m-dimensional origin (0,…,0) on a straight line.All existing primordia account for an energy potential on S m−1 (cf.Fig. 3), which, at the time set by the phyllochrone, determines the position where the newest primordium is formed in the same manner as in the two dimensional model.In this scenario, the model equations ( 1) are readily modified by replacing S 1 by S m−1 , and considering x, x 1 , …, x n+1 ∈ S m−1 .
We aim at characterizing the structures arising from this model.A first step is to generalize the divergence angle.Refer, for example, to Strang [21] for details on the linear algebra used here.In m dimensions, we denote where, in general, denotes the vector subspace spanned by vectors v 1 , …, v n ∈ ℝ m , and where ∠U 1 ,U 2 denotes the angle between two vector subspaces U 1 and U 2 in the subsequent sense.The angle δ n (j) for j ∈ {2, …, m−1} is well-defined if x n+j and x n are each linearly independent of the j-1-dimensional vector subspace U j = span(x n+j−1 , …, x n+1 ).In this case, δ n (j) is given by the angle drawn by two vectors that are each orthogonal to U j , and whose linear span with U j equals span(x n+j , …, x n+1 ) and span(x n+j−1 , …, x n ), respectively.In the opposite case, it is reasonable to define it as π.For a given plastochrone λ, the angles δ n (1) , …, δ n (m−1) fully describe an arrangement modulo rotations.Fig. 1 illustrates δ n (1) and δ n (2) in three dimensions.
Already in three dimensions the system proves to be considerably less robust than in 2D, making a very fine grid of S 2 necessary.In most simulations we worked with ca.2.6 × 10 7 quasi-regularly distributed grid points (in 2D a grid size of order 10 3 is largely sufficient).At this resolution, the particular geometry of the grid did not create Fig. 1 The two generalized divergence angles δn (1) and δn (2) in three dimensions.Like in the two-dimensional case, δn (1) is the angle between the locations of the two most recent primordia, with respect to (0,0,0).δn (2) is given by the angle between the two planes that are spanned by xn, xn−1 and (0,0,0), and xn−1, xn−2 and (0,0,0), respectively.As such it is equal to the spherical angle between the arcs (on S unwanted artifacts as demonstrated by simulation results proving robust against small random perturbations of the grid.Using an even finer grid in simulations for selected, small values plastochrones yielded identical results to those obtained with the standard grid.The step size of the quasi-static decrease of λ towards 0 was chosen sufficiently small (down to 0.001) in order for the resulting patterns to always respond smoothly to any change, especially near the transitions of qualitatively different patterns.In this article, we present results for m = 3 with perspectives on higher dimensional spaces in the "Discussion".We will first focus on the case of p = 2 for the energy function, and discuss other values later.The particular choice p = 2 allows computationally considerably faster simulations.This is because computing ||x − (1 + t n + λ − t i )x i || −p in (1) and its higher-dimensional analogue for any p, involves, by definition of the norm, computing ||x − (1 + t n + λ − t i )x i || 2 first.The potential additional exponentiation by −p/2 is computationally very intensive.We thus ran baseline simulations using p = 2, and then modified specific aspects from there, as detailed below.

Results
In this section we describe the model dynamics of the three-dimensional analogue of (1) in the course of decreasing the plastochrone λ quasi-statically from an initially high value.

Linear and planar arrangements
Just as in the two-dimensional case, the energetic influence on S 2 of all but the youngest primordium is negligible for all large enough λ at the moment of the formation of the newest primordium.Hence, it appears opposite to the last one, so that δ n (1) = δ (1) = π in the equilibrium (Interval I in Fig. 2a, Fig. 2c).In this arrangement, δ n (2) is not well-defined in (2), and thus set to δ (2) = π.With decreasing λ, the influence of the second-to-last primordium increases, whereas the influence of even older primordia is at first still negligible.It is easy to show that in this scenario, starting from λ ≈ 5, the newest primordium will form on the circle given by the intersection of S 2 and the plane through (0,0,0) that is spanned by the last two primordia (Fig. 2d).Thus, we observe the quantitatively identical planar spiral structures as in the two-dimensional case, and, without loss of generality, δ n (1) = δ (1) < π in the equilibrium, while unchangedly δ (2) = π by definition (Interval II in Fig. 2a).
The spatial arrangement of primordia in ℝ 3 for given λ in Interval II is identical to the one in two dimensions for the same parameter.In particular, the arrangement at some time t n coincides with the arrangement at t n+1 modulo a rotation by δ (1) about the normal of the plane in which primordia lie, and appropriate (planar) parastichy spirals can be drawn (cf.[10]).  (1(green) and δ' (2) (blue; see text) against those λ for which they converge (a).The dotted black line correspond to δ in 2D, which converges for all λ > 0 and coincides which δ (1) in Intervals I and II.b infn δn (1) and supn δn (1) (green) and infn δ'n (2) and supn δ'n (2) (blue), as well as, dotted, the respective means in the appropriate λ range (see text).Bottom: the set {xn}n in the modified Eq. ( 1), i.e., the set of all points on S 2 where primordia are formed in the equilibrium arrangement for λ in Interval I (c), II (d), III (e), IV (f).g Top view of f.
Geometrically, this implies that primordia appear alternately on two parallel and equal of size small circles of S 2 (Fig. 2e), before moving away from the origin on a twodimensional circular conical surface.
In two dimensions, a so-called parastichy change is observed around λ ≈ 1.8, where the newest primordium appears no longer closest to the last and second-to-last primordium, but to the second-to-last and third-to-last since, as λ decreases, the energetic influence of the third-to-last primordium eventually predominates that of the last primordium [10].This is different in three dimensions: by appearing alternately on the two small circles, before subsequently moving along the conical surface, new primordia maintain their distance to all three previously formed primordia.Fig. 3a illustrates the corresponding energy potential on S 2 .Fig. 3 Illustration of the cumulative energy potential on S 2 for the stable arrangement for λ = 1.7 (a) and λ = 0.198 (b) (both viewed from two different sides).Yellow corresponds to the highest energy potential on S 2 , dark blue to the lowest (for clarity, the color scale corresponds to the logarithm of the actual potential).Small black circles illustrate previously formed primordia (connected to (0,0,0) to help clarify their position in space) that jointly generate the potential.The resulting position of the formation of the newest primordium is colored in red.
The conical arrangement is maintained as the energetic potential generated by the fourth-to-last primordium increases, whereas the fifth-to-last primordium remains negligibly far away.Fig. 4 illustrates the neighborhood relations for λ in the Intervals I to II.
Regarding the spatial arrangement of primordia in ℝ 2 for λ in Interval III, we observe spiral patterns similar to the parastichies in 2D when connecting primordia formed on either one of the small circles in a suitable way (Fig. 5a).This results from the fact that points corresponding to either even or odd iteration steps (points on red and blue lines in Fig. 5, respectively) feature a phyllotactic pattern of their own on the appropriate single cones.
Consecutive arrangements in equilibrium for λ in Interval III are equal up to a rotary reflection, i.e., a reflection through the plane separating the two cones on which primordia lie and an appropriate rotation about the axis perpendicular to this plane.Analogously, in the two-dimensional model (and for all λ > 0), an equilibrium arrangement at some iteration step n is equal to the arrangement at step n−1 after a rotation about the divergence angle.

Oscillating divergence angles
The dissolvement of the conical arrangement is triggered by the approach of the fifthto-last primordium starting at λ ≈ 0.2 down to λ ≈ 0.197 (Interval IV in Fig. 2b).The corresponding energy potential on S 2 is shown in Fig. 3b.The two sequences δ n (1) and δ' n (2) , hitherto convergent, begin to oscillate (Fig. 6a) in the equilibrium states, with amplitudes increasing with decreasing λ as shown in Fig. 2b.
In Fig. 7, successive elements are not connected by lines as in Fig. 6a (the general oscillatory pattern is the same though), which highlights interesting patterns, the perhaps most remarkable of which is observable between λ ≈ 0.2 and λ ≈ 0.198 and readily visible in Fig. 7b,c: the elements of the sequence δ n (1) appear to be uniformly arranged along a set of periodic smooth functions differing by translation.Fig. 6b shows the appropriate section of Fig. 7b    for the appropriate parameter range, where a = sup n δ n (1) − inf n δ n (1) (increasing with decreasing λ), c = 0.5(sup n δ n (1) + inf n δ n (1) ), and b are functions of λ. a and c can be inferred from Fig. 2b, while the scaling function b is in fact nearly constant in this range of λ (not shown).Note that the number seven of distinct subsequences remains the same throughout this parameter range.We have not yet been able to explain the occurrence of this particular number.In Fig. 7, the sequence δ' n (2) (λ) features a distinct pattern as well, resembling an overlap of several sine-related functions.We omit its categorization, which is more extensive and bears little additional value, at this point.
The number 7 appears also in the spatial distribution of primordium origins on S 2 (Fig. 2f): In comparison to the small circles observed for λ in Interval III, waved lines emerge, each with seven distinct crests and troughs.Indeed, we also observe a 7-periodicity in the neighborhood relations (i.e., the analogue of Fig. 4, not shown).As for the arrangement of primordia in ℝ 3 for λ in Interval IV, spiral structures can still be identified, however the lines are "wavy" (Fig. 5b).

Fluctuations and chaos
As can already be seen in Fig. 7d, the regularity of the divergence angle pattern begins to dissolve around λ ≈ 0.197.A certain periodic structure is still visible, but less distinct than for larger λ.This trend continues with decreasing λ, until eventually no structure whatsoever is apparent and the sequence appears to fluctuate chaotically (not shown).
During this process, the subset of S 2 where new points are formed (Fig. 2f,g) loses its regularity as well by continuously "spreading" from the one-dimensional sine structure in Fig. 2f,g towards the equator and the poles of the sphere, as λ decreases (Fig. 8a-c), eventually covering all of S 2 (Fig. 8d).
We propose an informal argument as to why there does not seem to be a spiral-like pattern in the spatial distribution of primordia for all sufficiently small λ.In the 2D version of the model, (1), for an arbitrary but fixed Fibonacci number F n , F n smooth, congruent and disjoint spirals originating from (0,0), intersecting S 1 at the vertices Fig. 7 δn (1) (blue) and δ'n (2) (green) for λ = 0.2 (a), λ = 0.199 (b), λ = 0.198 (c), λ = 0.197 (d).

Fig. 6 a
The oscillating pattern of δn (1) in the stable state for λ = 0.198, where lines connect successive elements.b The same sequence with seven appropriate subsequences connected and colored, highlighting an underlying sinelike pattern.
of an S 1 -inscribed regular F n -sided polygon, and passing through all existing primordia can be drawn for all sufficiently small λ.The set of unit tangent vectors of these spirals at the points of their intersection with S 1 can be seen as a finite subset of an appropriate vector field on S 1 (Fig. 9).For any fixed F n , these vectors tend to become tangent to S 1 as λ → 0. Choosing increasingly large F n in this process, we approach indeed a tangent vector field on S 1 .Now, if, in the three-dimensional case, primordia arrangements for any arbitrarily small λ could be captured by a set of spirals with the same properties as listed above, this would be in contradiction to the hairy ball theorem which rules out the existence of non-vanishing tangent vector fields on S 2 .This brings us to reason that an analogue of the particular 2D structure for sufficiently small λ cannot exist in three dimensions.This argument does not exclude the possibility of a different underlying regularity, which we have not succeed to identify.
Other primordia speeds, energy potentials, and initial conditions Primordia speed.We reran simulations under the premise that primordia do not move away from the origin with constant speed, but with their velocities increasing linearly as argued to be biologically reasonable by Douady and Couder [10].Thus, a primordium is located exp(t) length units away from the origin after existing for t time steps.Fig. 10a summarizes the analogue of Fig. 2a,b.The results are qualitatively identical to the ones for primordia moving at constant speed, including the 7-periodicity of the divergence angle series (not shown).This may be argued to be expectable given that for sufficiently small plastochrones λ, the existing primordia that determine the position of the newest one, i.e., those relatively close to S 2 , will be at distances from the origin that are similar to those in the case of constant speed since in that case exp(t) ≈ 1 + t.

Energy potential.
In two dimensions, model dynamics are qualitatively identical for different exponents p of the energy potential [10], which we would also expect for the three-dimensional case.Indeed, the stages passed in the process of decreasing an initially high plastochrone λ towards 0 as described hitherto, namely the distichous, spiral planar, and conic arrangements as well the transition to chaos, were equally observed for a sample of values p ≠ 2.More quantitatively, we do not observe the increasing corrugation of the small circle structure into sine-like curves shown in Fig. 2d; the forced expansion of the domain where primordia are formed occurs in different, somewhat less striking ways (Fig. 11b,d).By contrast, the sine-related pattern of the divergence angles was observed for larger p (Fig. 11c), as opposed to irregular behavior for low values (Fig. 11a).In all cases, we observe chaotic fluctuations for all sufficiently small values of λ.We also reran our baseline simulations with the additional assumption of a time-dependent decrease of the energy potential induced by primordia (cf.[9]).Specifically, the energy potential of a primordium t time units after it was formed is multiplied by exp(−t).Fig. 10b illustrates the results, which are again, qualitatively identical to the previous simulations.

Initial conditions.
In their twodimensional simulations, Douady and Couder [10] observed additional branches, i.e., continuously varying pairs (λ,δ(λ)) of plastochrones and corresponding stable divergence angles, beyond the main Fibbonacci branch that is obtained from a quasi-static decrease of λ starting from a sufficiently high value (dotted black line in Fig. 2a).These secondary curves, such as the so-called Lucas branch, only exist for suitably small λ, and are obtained using suitable initial primordia arrangements.We ran our simulations (with constant primordium velocity and potential E(d) = d −2 ) for a range of suitably small plastochrones using a series of random initial conditions, but did not find equivalents of the secondary branches observed in two dimensions.In the majority of cases, the pattern converged to the appropriate arrangement in Fig. 2a,b.In the remaining cases, we observed periodic repetitions in the divergence angle series that we non-smooth and without particular structure (not shown).Douady and Couder [10] (their Fig. 5b) mention the same feature in two dimensions.

Discussion
A series of unexpected complex patterns emerged in the simulations of our generalization of the phyllotactic model (1) to three dimensions.Whereas in the 2D case the divergence angle converges for all plastochrones λ > 0, here, we observed a series of bifurcations as λ decreases, yielding convergence, periodicity and chaos of  the generalized divergence angles.This illustrates some of the remarkable dynamical properties that the modified Hofmeister rule can generate from a mathematical viewpoint.
While computational intensity currently prevents us from investigating higher dimensional cases, the previous observations already allow certain preliminary insights.In general ℝ m , all sufficiently large λ imply the alternate formation at opposite points of S m−1 , i.e., the intersection of S m−1 with a one-dimensional subspace.For smaller λ, the increasing influence of the second-to-last primordium brings about the spiral arrangement in a two-dimensional subspace.This process of decreasing λ leading to the gradual expansion or the spatial arrangements into higher dimensions continues until the dimension of the space considered is reached.As of this point, locations of primordia formations begin to gradually cover the whole of S m−1 (not just lowerdimensional subsets).Of particular interest are model dynamics for sufficiently small λ.Our conjecture for the occurrence of chaos in the three-dimensional case based on the hairy ball theorem would not apply to higher even dimensions, since the theorem does not hold there, which leads us to expect a regular behavior as λ approaches 0 for even dimensions m ∈ 2ℕ.The emergence of the golden angle φ, which, in our simulations we mentioned only briefly for λ at the lower boundary of Interval III in Fig. 2, will be particularly interesting to investigate in these scenarios.
By construction, models based on the modified Hofmeister rule cannot generate decussate, whorled, bijugate and multijugate patterns, which are characterized by several primordia forming at the same time.The alternative hypothesis by Snow and Snow [22] that a new primordium does not necessarily form after a particular time interval, but instead at the time and position on the apical ring where sufficient space is available proved to allow for these patterns (cf.[23] for a mathematical analysis).Douady and Couder [15,24] presented a dynamical model based on Snow and Snow's hypothesis, in which a new primordium forms at the time and position when and where the cumulative energy potential emitted from previously formed primordia falls below a given energy threshold, which replaces the phyllochrone as the key system parameter.Since the correct time and position of primordium formation can only be determined by means of iterative evaluations of the total energy potential along the apical ring, simulations of this approach are computationally considerably more intensive.This has prevented us from simulation the otherwise straightforward analogue of the model by Douady and Couder [15,24] in higher dimensions.These authors argued that for every n ∈ ℕ there are is an interval of energy potential parameters such that in this range the whorled mode in which n primordia form at the same time, i.e., in the shape of a regular n-gon, is a stable arrangement.The fact that there are only finitely many regular polyhedrons in higher dimensions, in particular the five platonic solids in three dimensions, shows that a direct analogue of this observation cannot exist in higher dimensions.Nevertheless, studying the interaction of the existing regular polyhedrons in this case would be very interesting.
Fig.1The two generalized divergence angles δn(1) and δn(2) in three dimensions.Like in the two-dimensional case, δn(1) is the angle between the locations of the two most recent primordia, with respect to (0,0,0).δn(2) is given by the angle between the two planes that are spanned by xn, xn−1 and (0,0,0), and xn−1, xn−2 and (0,0,0), respectively.As such it is equal to the spherical angle between the arcs (on S 2 ) from xn to xn−1, and from xn−1 to xn−2, as shown in the figure.

Fig. 4
Fig. 4 Distances of existing primordia (at positions xn, (1 + λ)xn−1, (1 + 2λ)xn−2, … ∈ ℝ 3 ) to the newest one (at xn) at the moment of its formation in equilibrium for λ ∈ [0.2, 5](resized by λ −1 for the purpose of visualization).Below λ ≈ 0.2 these stable relations begin to oscilate.Whereas the increasing proximity of the fourth-to-last primordium does not alter the conic arrangement, the influence of the sixth-to-last at λ ≈ 0.2 induces oscillations in the divergence angles and the spatial arrangement.

Fig. 9
Fig. 9 2D primordia arrangement for λ = 0.01 along with eight spirals and the corresponding vector field.

Fig. 10
Fig. 10 Analogue of Fig. 2a,b in the case of (a) linearly increasing primordium velocity and (b) time-dependent decay of the energy potential generated by primordia.The parallel branches for small λ show, as in Fig.2b, the supremum and infimum of the divergence angles.For all sufficiently small values, the pattern is, again, chaotic.

Fig. 11
Fig. 11 Locations of primordia formations for simulations using p = 1 and p = 4 for the exponent of the energy potential, for parameters in the analogues of Interval IV, specifically λ = 0.255 and λ = 0.185, respectively.