Universality in Phyllotaxis: a Mechanical Theory
aa r X i v : . [ n li n . PS ] J a n Universality in Phyllotaxis: a Mechanical Theory
Hyun-Woo Lee and Leonid S. Levitov
Physics Department, Massachusetts Institute of Technology,77 Massachusetts Avenue, Cambridge, MA02139 (Dated: December 2020)
Abstract:
One of humanity’s earliest mathematicalinquiries might have involved the geometric patterns inplants. The arrangement of leaves on a branch, seeds ina sunflower, and spines on a cactus exhibit repeated spi-rals, which appear with an intriguing regularity providinga simple demonstration of mathematically complex pat-terns. Surprisingly, the numbers of these spirals are pairsof Fibonacci numbers consecutive in the series 1, 2, 3, 5,8, 13, 21, 34, 55... obeying a simple rule 1+2=3, 2+3=5,5+8=13 and so on. This article describes how physicshelps to clarify the origin of this fascinating behavior bylinking it to the properties of deformable lattices growingand undergoing structural rearrangements under stress.
Contents
I. Introduction
II. Mechanical theory
III. The geometric model
IV. The energy model
V. Universality of Fibonacci numbers
VI. Robustness of the problem versus anisotropic growth 19B. The disk and cone geometries 20
VII. Summary References published in: Symmetry in Plants, Series in MathematicalBiology and Medicine, eds. R. V. Jean, D. Barab´e, (WorldScientific Pub Co Inc, 1998) I. INTRODUCTIONA. The patterns of repetitions in plants
The visual beauty of plants comes, in part, due to thehighly regular and well organized patterns of leaves, flo-rets, seeds, scales and other structural units. Widely ad-mired, the regular patterns in plants are not merely aes-thetically pleasing, they often have unexpected relationwith mathematics. Here we will be concerned with oneparticular kind of such patterns — the spirally, or helical,arrangements. Spirally patterns occurring in plants havelong been known to have a surprising connection withnumber theory. The area of botany that studies suchpatterns is called phyllotaxis (the word can be translatedas “leaf arrangement”). It has been recognized long agothat the so-called Fibonacci sequence: 1, 2, 3, 5, 8, 13,21, 34, 55, 89 ..., where every number in the sequenceappears as a sum of two preceding numbers, is of greatimportance in phyllotaxis. The prominence of Fibonaccinumbers in phyllotaxis is well accounted for in both spe-cialized and popular literature, which includes the gemssuch as “On Growth and Form” by D’Arcy Thompson [1]and “Symmetry” by H. Weyl [2].The connection between Fibonacci numbers and heli-cal packings of units in plants, which is at the heart ofthe subject of phyllotaxis, appear to be quite general,although the details depend somewhat on the plant ge-ometry. In cylinder-shaped objects, such as fir-tree conesor pineapples, the scales have a regular arrangement inwhich two families of helices can be identified, having theright and left helicity. These are known in the literatureas parastichy helices [3] (see Fig. 1a). The numbers ofhelices in each family are invariably found to be the Fi-bonacci numbers. Moreover, the numbers of the rightand left helices are non-equal and are given by the pairsof Fibonacci numbers consecutive in the sequence, suchas (2,3), (3,5), (5,8), (8,13) and so on. Geometrically,such structures can be described as lattices on the sur-face of a cylinder. The appearance of Fibonacci numbersin a cylindrical structure is called cylindrical phyllotaxis .Another type of phyllotactic patterns occurs whenunits of a plant are packed on a disk. For example, in asunflower or a daisy, the lines connecting neighboring flo-rets define two families of parastichy spirals [3]. Again, ineach family the number of spirals is Fibonacci, and in thetwo families the numbers are consecutive in the sequence(Fig. 1(b)). For such structures, it is common to use thename spiral phyllotaxis . A detailed discussion of the twogeometric models of phyllotaxis, cylindrical and spiral,
FIG. 1: (a) Geometric model of a cylindrical plant, such asa pineapple or a pine cone. The lines connecting neighboringscales form two families of parastichy helices with the rightand left helicity. Each family contains a Fibonacci number ofhelices. The figure illustrates a model of such a structure by alattice on a cylinder of circumference a , and its mapping ontoa periodic lattice in the plane. The two lattice vectors areshown that connect the nearest and next-nearest neighbors,and define parastichy helices, The divergence angle d and thehelix rise r are defined in Eq.(1); the parameters x = − d / π and y = r /a are defined in Eq.(4). (b) Geometric model of aplant with a disk geometry, such as a sunflower. The lines con-necting neighboring florets generate two families of parastichyspirals of the right and left helicity. The figure illustrates amodel of the floret arrangement by a spiral lattice. The pointson the spirals are numbered according to the distance fromdisk center. can be found in the article by Rothen and Koch [20].The first studies of phyllotaxis probably go back to theobservation made by Leonardo da Vinci (notebook, 1503)who noted thatNature has arranged the leaves of the latestbranches of many plants so that the sixth isalways above the first , ...(Here 6 − B. Why Fibonacci numbers?
In the XX century the trends and interests in this fieldfocused on the question of the origin of phyllotaxis. Hereis how the challenge was stated by a mathematical biol-ogist: The fascinating question: “Why does the Fi-bonacci sequence arise in the secondary rightor left spirals seen on plants?” seems to beat the heart of problems of plant morphol-ogy. In atomic physics, Balmer’s series hasopened the way to Bohr’s theory of the atomand then to quantum mechanics and to quan-tum electrodynamics. The great hope of bio-mathematicians is that some day they maybe able to do for biology what has been doneby mathematical physicists in physics. (R. V.Jean, Ref. [4], p. 641)Indeed, despite the long history, no general agreementwas reached with regard to the origin of phyllotaxis.Many explanations proposed in the past simply state thatthe Fibonacci packings are optimal in some sense. Al-ready Leonardo, after noticing the regularity in the num-bers of leaves, remarks that the reason they are arrangedin such a way probably has something to do with a bet-ter exposure to the sun light through minimizing mutualshadowing. Later, the logic of the discussion of the ori-gin of phyllotaxis basically followed that of Leonardo’s,by adding more optimization factors, such as air circula-tion in between the leaves, or the density of packing of theseeds or florets. The review of such theories can be found,for example, in the D’Arcy Thompson’s book [1]. TheDarwinian evolutionary theory, by asserting that naturalselection gives rise to the optimal structures that are bet-ter designed for a local environment, apparently endowssuch theories with an air of confidence.However, such views are known to be at odds with cer-tain experimental observations. Indeed, as emphasized inthe evolutionary theory, the evolutionary pressure arisesfrom variability within the species. Namely, if the speciesare optimized with respect to a certain factor by evolu-tionary pressure, there must also occur, however rarely,the species that are not optimal but near-optimal. Ifthis was the case in the phyllotactic growth, a sunflowerthat normally exhibits 55 spirals could sometimes have54 or 56 spirals. In other words, the most frequent kindof an exception from the Fibonacci rule would be asso-ciated with the numbers that are close but not equal tothe Fibonacci numbers. And yet such numbers never oc-cur; instead a very different kind of non-Fibonacci pat-terns is observed. The most common exception knownto occur is described by the numbers from the so-called
Lucas sequence [2]: 1, 3, 4, 7, 11, 18, 29 ... . The se-quence is constructed according to the Fibonacci additionrule, however it starts with a different pair of numbers.Another aspect of phyllotaxis that seems to be hard toreconcile with the evolutionary optimization theories, isthat the non-Fibonacci numbers patterns, even those ofthe most common Lucas type, occur extremely rarely,with the probability of few percent or lower.The high stability of Fibonacci numbers could makeone suspect that perhaps they are in some way “pro-grammed.” However, the conjecture that these numbersare merely encoded genetically seems improbable, be-cause individuals within the same species often exhibitdifferent Fibonacci numbers, but hardly ever the non-Fibonacci (Lucas) numbers. Moreover, as in the case ofspiral phyllotaxis, different Fibonacci numbers can occurwithin one plant. In addition, by simply making referenceto genes one does not come closer to the understandingexactly what is special about these numbers. Thus oneis led to seek an explanation elsewhere.In this article we review a mechanical theory [15, 16]that establishes a general physical mechanism of phyl-lotactic growth. Namely, we consider the role of mechan-ical stresses in a densely packed arrangement of units ina growing plant. The stresses, which build up due to thegrowth anisotropy, lead to a shear strain in the structureincreasing gradually and then relaxing by a sequence ofabrupt structural rearrangements. As demonstrated be-low, this simple mechanical process gives rise, exclusivelyand deterministically, to Fibonacci structures.The mechanical explanation of phyllotaxis, which linksFibonacci numbers to transitions in the system upon stress build-up in the process of growth, is completelygeneral and robust. In particular, these transitions areinsensitive to the specific form of interactions and, oc-curring one after another, generate the entire sequenceof Fibonacci structures. Further, besides establishingthe prevalent character of Fibonacci numbers in phyl-lotactic patterns the mechanical theory explains why theexceptions are predominantly of a Lucas type. This isachieved by allowing irregularities (e.g. due to fluctu-ations or noise) to trigger a single mistake at an earlygrowth stage. All predictions of the mechanical theoryare therefore in agreement with the observations.Besides the mechanical scenario of phyllotaxis reviewedin this article, several other approaches focused on var-ious physical growth-related effects that may lead tophyllotaxis. Those include theories of growth medi-ated by a diffusion of inhibitor (Mitchison [5]), by areaction-diffusion process (Meinhardt [6]), and by me-chanical interactions that control the largest availablespace (Couder and Douady [17]). Some of these theo-ries are discussed in other chapters of this volume.Are the mechanisms emphasizing different physicalprocesses mutually exclusive or complementary to oneanother? On a first thought it may appear that, giventhe clear differences between these approaches, identify-ing the right mechanism will invalidate other explana-tions. Yet, below we argue that the situation is con-siderably more interesting. Our analysis of mechanicalstresses and their impact on growth establishes the prop-erty of robustness . The robustness basically means a widestability of phyllotactic patterns with respect to param-eter variation in the model. The stability property, byan extension, suggests that there is a degree of truth inall theories which invoke some form of repulsion/stressduring growth — mechanical, chemical, or else. In otherwords, all theories of phyllotaxis which use lattices of ob-jects with some kind of repulsive interactions, no matterthe origin and specific geometry (e.g. cylinder, disk, orcone), are equivalent in a “coarse-grained” sense. Thesituation here is similar to that in the theory of patternformation, where many different microscopic models areknown to lead to identical types of patterns on a largerscale. In physicist’s language, such models, while differ-ent in details, belong to the same universality class.At the same time, there is an open question of compar-ing the proposed scenarios with the underlying physiolog-ical processes in biological systems. This is a fascinatingexperimental problem that will have to be addressed byfuture work. However, whatever the outcome of thesestudies might be, it is worth noting that phyllotacticgrowth occurs in a large variety of biological systems. Itis therefore possible that no unique physiological processapplicable to all systems will emerge, and instead sev-eral different microscopic scenarios must be considered.At the same time, the robustness of phyllotactic growth,as established in the mechanical theory, will ensure thatthe resulting patterns are insensitive to the microscopicdetails of the growth.This article is organized as follows. In Sec. II we sum-marize the mechanical theory of phyllotaxis and providea brief review of its history. In Sec. III we introduce a ge-ometric model that involves cylindrical lattices and fam-ilies of parastichy helices, defined in terms of the shortestlattice vectors. In this section the notion of a phase spaceof all cylindrical lattices is defined, which will be centralin our subsequent discussion. Next, in Sec. IV we in-troduce the energy model, and consider its symmetries.These symmetries are found to form a large family, de-scribed as an infinite group of modular transformations.Then, in Sec. V the main result of this article is derived.We use the modular symmetries to relate different growthstages and thereby show that Fibonacci numbers are uni-versal in phyllotaxis. This result, established within theenergy model, is cast in a rigorous form of a theorem.From the discussion in Sec. V the robustness property ofphyllotactic growth becomes evident. In Sec. VI we dis-cuss possible modifications of the model which can bringit closer to other growth geometries, such as those of spi-ral phyllotaxis, and argue that the stability of Fibonaccipatterns and their universality remain unaffected.
II. MECHANICAL THEORYA. Growth under stress
The mechanical model discussed here involves, in itssimplest form, a regular cylindrical lattice of repulsivelyinteracting objects (see Fig. 1a). Cylindrical lattices pro-vide a convenient representation for a wide variety ofphyllotactic patterns. Crucially, in this model the latticesare taken to be deformable, with the lattice geometry notfixed rigidly but instead controlled by the balance be-tween external forces and repulsive interactions betweendifferent lattice points. The repulsion can be chosen sothat it mimics the contact rigidity of the structural unitsof a plant. However, rather than restricting the repul-sive interaction to be a short-range type, it is beneficialto consider a generic repulsive interaction which includesboth the short-range and the long-range parts.Further, the effects of growth can be naturally incorpo-rated in the cylindrical lattice model through an externalforce applied along the cylinder axis, which gradually in-creases as the growth progresses. The uniaxial stress dueto such a force mimics the stresses arising in the growth ofan elongated object in the presence of the closed-volumeor confinement constraints. As we will see, under theuniaxial stress the system deforms in such a way that,starting from a simplest quasi-one-dimensional chain-likestructure, it goes through the sequence of Fibonacci phyl-lotactic patterns in a completely deterministic way.In the mechanical model of phyllotaxis a cylindricallattice deforms upon a gradual increase in stress; this de-formation is described in terms of a gradually developingshear that tends to minimize the total repulsion energyat each given stress value. It is crucial, however, that in such a process the structure does not track the globalenergy minimum. The reason is that in optimizing itsenergy the system can explore only the nearby states bydeveloping small deformations in the lattice. Therefore,the Fibonacci structures, appearing upon a gradual in-crease of the stress, characterize the progression, or timesequence, of deformation, rather than the result of globalenergy minimization. The Fibonacci structures emergein this framework as the stages of rearrangement in a de-formable system, starting from a certain simple structureand appearing one by one upon an increase in the stressduring the growth.In treating the problem quantitatively, the first stepis to define the mechanical energy of a lattice using amodel interaction between lattice points, and use this in-teraction to study the system evolution under stress (seeSecs. III,IV). This analysis, which is straightforward tocarry out numerically [15], reveals a very robust behav-ior: for a variety of repulsive interactions, increasing thestress drives the system through all Fibonacci structures,exclusively and without exceptions. Further, the modelhas a distinct advantage in that it can be treated ana-lytically for a large family of interactions (see Sec. V).In this case the robustness and universality of Fibonaccistructures can be established rigorously [16]. We subse-quently argue (Sec. VI), using the robustness property,that the results obtained for cylindrical geometry can begeneralized to the disk geometry.Anisotropic mechanical stress accompanying thegrowth is a key ingredient of the theory of phyllotaxisadvocated in this article. What can be the origin ofsuch a stress? In that regard we note that, while therequired stress cannot originate from isotropic hydro-static pressure, any growth anisotropy can in generallead to anisotropic stresses of the form that generatephyllotaxis. The details of the relationship between thegrowth anisotropy and stress depend on the specifics ofthe system at hand and, in particular, the system geome-try. However, as argued below, these details are inessen-tial for understanding the general relation between thestress buildup during growth and the formation of phyl-lotactic patterns.The relation between anisotropic growth and mechan-ical stresses is most clear for the growth of a cylindricalstructure. Consider a cylindrical lattice, such as thatpictured in Fig. 1a, which is growing while being encap-sulated in a fixed volume. Such growth can describe, forinstance, the early developmental stages of objects likepine tree cones, which have a growth center at the apex.At the growth center, new structural units of the would-be cone are produced at a constant rate, and as a resultthe cone extends forward. In free space, such a growthwould yield one-dimensional chain-like structures. How-ever, if the growth is taking place within a closed volume,the growing cone, or a similar cylindrical object, soonmeets a constraint (namely, a boundary) which does notallow it to freely extend forward. As a result, as moreand more units are generated at the growth center, thecone will deform in order to fit inside the enclosure. It isinstructive to characterize the result of such deformationby rescaling to fixed density on the cylindrical surface(see Fig. 1a). Upon such rescaling, the growth can be de-scribed as a gradual expansion in the direction transverseto the cylinder axis accompanied by compression alongthe axis, such that the two-dimensional areal density re-mains fixed. As the our analysis below predicts, sucha process leads to, exclusively, the Fibonacci structures.The actual numbers achieved will depend on the growthduration: the longer the growth continues, the higher arethe Fibonacci numbers that can be accessed.What are the implications of the behavior found incylindrical lattices for other geometries of interest? Auseful example that helps to answer this question is thegrowth of a circular object in a disk geometry, as illus-trated in Fig. 1b. Such a growth, e.g. describing thedevelopment of a sunflower or a daisy, has its center atthe center of the disk. During the growth, as more units(e.g. florets) are added at the center, they push the previ-ously added units, forcing them to move outwards alongthe disk radius. One can argue that this process leads totransformations of the structure equivalent to those of adeformable cylindrical lattice.This can be done, e.g., by focusing on the evolutionof a small rectangular patch of the structure, a distance r away from the center. Namely, one can choose therectangle sides to be aligned with the cylindrical coor-dinates, radial and azimuthal, equal to δr and rδθ , re-spectively, and consider how they change upon growth.As the growth progresses, this patch moves radially to anew location r ′ > r , where the new sides of the rectan-gle become ( r/r ′ ) δr and r ′ δθ , since the rectangle angulardimension δθ and its area rδrδθ remain constant. Therectangle aspect ratio rδθ/δr increases thereupon by afactor ( r ′ /r ) . As a result, the part of the lattice withinthe rectangle, after moving from r to r ′ ≫ r , is stronglydeformed. The deformation degree ( r ′ /r ) gradually in-creases along the radius, whereas the areal density of theflorets remains roughly constant. Which implies that,since the main forces in the system are due to nearestneighbors exerting pressure on each other, locally the me-chanics of the deformed spiral lattice is similar to that ofa cylindrical lattice. To map one problem to the other,we can simply replace the spirals of the disk by the he-lices of the cylinder, which gives a one-to-one relationbetween the parastichy lines connecting the adjacent lat-tice points. One unique aspect of the disk growth ge-ometry that distinguishes it from cylindrical geometry isthat the deformation varying as a function of radius, asdiscussed above, gives rise to structural transitions anddifferent phyllotactic domains occurring within the samedisk (see discussion in Sec. VI). The concentric annulus-shaped phyllotactic domains have the numbers of spirals(parasticies) which are smaller for the inner domains andgreater for the outer domains, increasing with radius (seeSec. VI B). The numbers of spirals, which are Fibonaccinumbers, are therefore largest near the disk edge. Another geometry of interest corresponds to thegrowth of a plant at an apex of its shoot. In this case,the growth center where new units of a plant (e.g. leaves,scales, or spines) are generated is located at the centerof the rounded top part of a shoot. The geometry of acurved cone describing such a growth can be viewed asintermediate between the disk and cylinder geometriesdiscussed above. Near the growth center, the growthprocess can be described by the disk geometry, whereassome distance away from it, as the shape of the shootis curved towards a cylinder, it resembles the cylindricalgrowth. As in the disk growth case, a thicker shoot givesrise to bigger Fibonacci numbers.These three basic growth geometries represent differentvarieties of the general problem of a deformable latticegrowing under stress. All three exhibit a very similar be-havior. While in each case there are specifics that lead tosome modifications in the description of the growth, theydo not affect the stability of Fibonacci numbers. The de-velopmental stability follows from the robustness of thesequence of structural changes induced as the mechani-cal deformation increases under stress (see Sec. VI). Thisproperty is central to understanding the universality ofFibonacci patterns. B. Early mechanical scenarios
The first suggestion that mechanical forces play a keyrole in phyllotaxis seems to have been made by Hu-bert Airy. In his work “On Leaf-arrangements” (Ref. [9],p.177) he wrote:Take a number of spheres (say oak-galls) torepresent leaves, and attach them in two rowsin alternate order (1/2) along opposite sidesof a stretched india-rubber band. Give theband a slight twist to determine the directionof the twist in the subsequent contraction,and then relax the tension. The two rows ofspheres will roll up with a strong twist into atight complex order, which, if the spheres areattached in a close contact with the axis, theorder becomes condensed into (nearly) 2/5,with great precision and stability. And it ap-pears that further contraction, with increaseddistance of the spheres from the axis, will nec-essarily produce the order (nearly) 3/8, 5/13,8/21, etc. in succession, and that these suc-cessive orders represent successive maxima ofstability in the process of change from thesimple to the complex.”This description contains the key element of the mechani-cal theory: a spiral structure of rigid objects that deformsunder a stress. However, the way it is described does notmake it clear why such a process leads exclusively to Fi-bonacci numbers. Moreover, Airy’s statement seems tobe more of a conjecture rather than a conclusion, sincethe analysis verifying this suggestion has not been doneat the time. To put these ideas on a firm ground a theoryof the development of a spiral structure under mechanicalstress is needed.An important step in this direction was made by vanIterson [11]. As a starting point, he considered cylindricallattices, a geometric model that has been used to describephyllotactic patterns since the work of Bravais [10], towhich van Iterson added a new ingredient. As in the ear-lier work, a plant is represented as a cylinder with a spiralarrangement of lattice points on the surface of a cylinder.In that, the lattice points represent individual units of aplant such as leaves, scales, or spines (see Fig. 1). etc.The lattice points, if ordered according to their relativeheights along the cylinder axis, h m , m = 0 , ± , ± ... , canbe described geometrically as a set of points on a single generating helix : h m = r m, θ m = d m, (1)where r is the rise of the helix, and d is the angular stepknown as the divergence angle in the literature. Thequantities h and θ , and the lattice points labeled by m ,are illustrated in Fig. 1a; here all lattice points are as-sumed to have different height values h m . The new in-gredient of van Iterson’s model is identical disks centeredat the lattice points, with the arrangement of the disksconstrained by the dense packing requirement. Namely,the disks are packed in a lattice such that each disk is ina direct contact with at least four nearest-neighbor diskswhich touch but do not overlap.The densely packed disks resemble the densely packedstructural units of a plant, defining the families of rightand left parastichy helices which connect the nearest-neighbor lattice points as illustrated in Fig. 1a. The num-bers of the helices in the two families define the cylindri-cal lattice “parastichy type”. (This model, and in par-ticular the relation of the pairs ( r , d ) to the hyperbolicplane, will be discussed in greater detail in Sec. III.) VanIterson noticed that the dense packing requirement con-strains the disk arrangements in a very special way. Hedemonstrated that the lattices of densely packed disks,when described in terms of the generating helix, Eq.(1),define possible pairs ( r , d ) that form a branching struc-ture known as the Cayley tree (see Fig. 2). The Cay-ley tree consists of the arcs of circles in the ( r , d ) planeconnected with each other at the triple branching pointsoccurring at the ends of the arcs. The triple branchingpoints correspond to the maximum-density disk packing— the symmetric triangular lattices in which each diskmakes contact with exactly six neighbors. Further, theCayley tree helps to understand the properties of thedensely packed disk lattices encoded in their parastichytypes — the numbers of the right and left helices whichare uniquely determined by the values of r and d .These numbers show an interesting behavior whenplaced on the corresponding branches of the Cayley tree.First, the numbers do not change within one branch.Each branch can therefore be labeled by a single pair of FIG. 2: Van Iterson diagram of pairs ( r , d ) corresponding todense packing of disks [11]. The diagram has a structure ofa Cayley tree with the branching number z = 3. Parastichypairs (the numbers of the right and left helices) are shown nextto each arc. Inset: The contact pressure theory by Adler [14]predicts lifting of the 3-fold symmetry at the branching pointsby opening gaps in some of the arcs. The dotted line repre-sents unstable disk packings with negative contact pressure. integers ( n , n ). This is illustrated in Fig. 2 where thesenumbers are shown next to the corresponding Cayley treebranches. Second, as noted above, the triple branchingpoints correspond to perfect triangular lattices of maxi-mum density. Since in a perfect triangular lattice thereare three shortest lattice vectors of equal length, and thusthree distinct families of helices associated with these vec-tors, each branching point of the Cayley tree is naturallylabeled by three integers ( n , n , n ). As discussed be-low, these integers obey the identities n + n = n , upto a permutation, which might suggests a link to the Fi-bonacci series.Is there a relation between this construction and Fi-bonacci numbers? A quick inspection of Fig. 2 indicatesthat the numbers labeling the arcs are in general non-Fibonacci, with the Fibonacci numbers found only in asmall subset of the tree. At the same time, the arcs thatare labeled by the pairs of Fibonacci numbers form twocontinuous paths going from the Cayley tree top all theway down to the bottom. Along these paths the pairsof consecutive Fibonacci numbers are found (1 , , , ... , and (2 , , , ... ; the numbers grow asthe rise parameter r of the generating helix decreases.Taking the decrease in r to be associated with the plantgrowth (e.g., see Ref. [23]), it is tempting to conjecturethat the arcs of the Fibonacci paths are somehow linkedto different stages of the growth.These observations were first connected to the effectsof mechanical stress by Adler [14]. He considered a forcethat is being applied to a cylindrical lattice along thecylinder axis and increases gradually. Under growingstress the lattice deforms so that the rise of the helix r gradually diminishes. In the disk model the effect ofgrowing stress is accounted for by a gradual downwardmovement along the Cayley tree of the point ( r , d ) thatrepresents the lattice. As r decreases, the evolution of thepoint ( r , d ) tracks one of the tree branches until a branch-ing point is reached. At which point the system mustchoose between the two arcs going further down. If weassume that the system, represented by the point ( r , d ),would always choose, by some mechanism, to evolve fromthe branching point along the arc marked by a pair ofFibonacci numbers, the model would predict phyllotac-tic growth. However, in the original formulation of thedisk model it was unclear what mechanism could guidethe growth in the right direction at the branching points.Adler proposed to replace the purely geometric con-dition of a dense packing of the disks by a more phys-ical notion of an external pressure which is applied tothe system of rigid disks and is balanced by the contactpressure between the disks (see also [12]). The contactpressure theory “resolves” the triple points on the Cay-ley tree diagram; this happens because the condition formechanical equilibrium predicts that one of the arcs thatmeet at each branching point represents a system withnegative contact pressure. The negative pressure regionscorrespond to unstable structures and must be eliminatedfrom the diagram. This leads to gaps opening in the arcswhich divide the tree into separate paths of arcs that rundownward (see inset of Fig. 2). The evolution along eachof these paths, owing to the absence of branching, is de-terministic. Further, the two Fibonacci paths discussedabove remain unaffected by the negative pressure gaps.Now, as r diminishes, the system can evolve through thesequence of Fibonacci states in a perfectly determinis-tic fashion. (Further details of this construction are dis-cussed in the chapter by Adler in this volume.)To summarize this discussion, by introducing contactpressure between the disks Adler completed and solvedvan Iterson’s disk model. His solution provided a proofof phyllotaxis for the contact-pressure interaction. Whilerealistic, this interaction is not found in any biologicalsystem. Moreover, since the phenomenon of phyllotaxisis so widespread, and the form of interaction presum-ably varies strongly among different species, a solutionfor a particular interaction only partially addresses theproblem, leaving the question of robustness unanswered.Proving the universality of phyllotaxis requires estab-lishing the stability of Fibonacci growth in a sufficientlygeneric class of models.This is precisely what is accomplished by the energymodel [15, 16]. The energy model is essentially a gen-eralized contact pressure model, with the contact inter-actions between disks replaced by a generic repulsive in-teraction between points of a cylindrical lattice. A bigadvantage of the energy model is that it cleanly delin-eates the physics ingredients (interactions) from the ge-ometry ingredients (cylindrical lattices). This makes itpossible to treat the problem of phillotaxis in full rigorand generality. These results will be summarized and reviewed below. Quite unexpectedly, the separation ofinteractions and geometry reveals that the problem pos-sesses a hidden symmetry, namely the modular symmetrygroup GL (2 , Z ). With the help of the GL (2 , Z ) symme-try the occurrence of Fibonacci numbers in phyllotacticpatterns can be established for generic repulsive inter-actions. This analysis demonstrates the universality ofFibonacci phyllotaxis, providing an explanation of theFibonacci growth.Further insight into the mechanical origin of phyl-lotaxis was provided by an experimental work publishedrecently by Couder and Douady [17] (see their chapter inthis volume). They devised a hydrodynamic system thatmodels spiral phyllotaxis with the help of magneticallypolarized droplets of a ferrofluid on an oil surface. In thisexperiment [17], the ferrofluid droplets appear in a reg-ular sequence at the middle of an oil disk, representingunits of a phyllotactic pattern generated at the growthcenter. The droplets repel each other by the dipole r − interaction; this gives rise to an effective pressure fromthe newer droplets on the older ones that organizes themin Fibonacci patterns similar to those seen in the spi-ral phyllotaxis. This work proves that, as envisioned byAiry, mechanical forces are indeed sufficient to producephyllotactic patterns.In addition to the non-biological realization in the oildrop experiment, it was proposed that phyllotaxis canbe observed in other physical systems, such as cells ofBenard convection [22] and flux lattices in supercon-ductors [15]. III. THE GEOMETRIC MODELA. Cylindrical lattices and parastichy helices
In this section we review in greater detail the cylindri-cal lattice model of phyllotaxis. This model provides, inparticular, a definition of the lattice parastichy type givenby the numbers of the right and left helices. In our discus-sion we will focus on the dependence of these quantitieson lattice geometry. Further, we will establish a relationbetween the phase space of cylindrical lattices and thehyperbolic plane. This relation will prove quite usefullater in the analysis of the energy model. Our discussionof cylindrical lattices overlaps in part with Refs. [18–20].Cylindrical lattices can be described by mapping themto lattices in a plane. For a given cylindrical lattice (1)the planar lattice is obtained by unrolling the cylinder asa wallpaper roll. Upon unrolling a two-dimensional lat-tice is obtained that has periodicity given by the cylindercircumference. Namely, working in Cartesian coordinates r = u i + v j , (2)we align the cylinder axis parallel to the j axis and rollit along the i axis direction. Accordingly, each point ofthe generating helix (1) is mapped on a one-dimensionalarray of points parallel to the i axis, r p = r + pa i , with a the cylinder circumference. In this way we obtain atwo-dimensional periodic lattice: r pm = u pm i + v pm j = a (( p − mx ) i + my j )= √ A (cid:18) p − mx √ y i + m √ y j (cid:19) , (3)where p and m are integers, A = a r is the unit cell area(see Fig. 1a). The first line of Eq.(3) gives the planar lat-tice obtained by unrolling the cylinder; in the second linenew parameters are introduced that will be used through-out our discussion below. These are the relative row-to-row displacement x of the generating helix, Eq.(1), andthe height-to-circumference ratio y , related to the param-eters r , d as x = − d / (2 π ) , y = r /a. (4)The helical lattice (1) is completely specified by threeparameters d , r , and a . Accordingly, the lattice (3) ischaracterized by the parameter values x , y , and A .The representation involving the quantities x and y ,which is used in our discussion below, has a number ofadvantages, some obvious and some less obvious. First,this representation is scale-independent (i.e. is invariantupon rescaling), since x and y depend only on the anglesbetween the basis vectors of the lattice and their relativesizes, but neither on the cylinder circumference a nor thelattice unit cell area A .Further, as demonstrated below in Secs. IV,V, theseparameters present an intrinsic advantage from a geo-metric viewpoint, helping to link cylindrical lattices tohyperbolic geometry. Namely, the lattice geometry canbe described by a single complex parameter z = x + iy (5)that can be treated as a variable in the hyperbolic plane.The energy of the lattice will be shown to be invariantunder modular transformations z → z ′ = ( az + b ) / ( cz + d ) . (6)This result will be crucial for our analysis of hidden sym-metries of the lattice energy in Sec. V, and for establish-ing stability of the Fibonacci phyllotaxis.In a lattice (3) the parastichy helices are introducedwith the help of shortest lattice vectors. In a generic lat-tice, each point has two nearest neighbors, where “near-est” refers to the metric in the Cartesian ( u, v ) plane (2),∆ = p ( u − u ′ ) + ( v − v ′ ) . Alternatively, the distancecan be measured on the curved cylinder surface usinggeodesics. This defines, up to a sign, a lattice vector r pm connecting each lattice point to its nearest neighbor.Similarly, there is a pair of next-nearest neighbors thatdefines a vector r qn . Generally, | r qn | > | r pm | , however,for the lattices with a rhombic unit cell the lengths oftwo vectors are equal. The vectors r pm and r qn will becalled the pair of shortest vectors . Given r pm , for any lattice point r we can define thecorresponding parastichy line by drawing a straight linethrough r and the nearest lattice points r ± r pm . Thisgives a line r + s r pm with s a real-valued parameter, whichcontains an array of lattice points { r + k r pm } , where k is an integer. On the cylinder, this line corresponds to ahelix drawn through r and connecting it to nearest neigh-bors. Different helices obtained in this way for different r ’s form a parastichy family . Likewise, the second paras-tichy family is defined by connecting lattice points in thenext-shortest vector r qn direction. The parastichy num-bers P , P which define the lattice parastichy type aredefined as the numbers of the helices in these two families(see also Refs. [18–20]).The definitions above are of course all but natural asthe helices on a cylindrical plant picked by a human eyedo connect nearest neighbors. In the remaining part ofthis section we review several useful properties of theshortest lattice vectors and their relation with the num-bers P , P . Namely, for the shortest and next-shortestvectors r pm and r qn given by (3) the numbers P and P equal, up to permutation, | m | and | n | .First, we show that the vectors r pm and r qn are prim-itive vectors, i.e. they provide a basis for the lattice (3).Although this property is nearly obvious, we sketch aquick proof both for completeness and as a reminder ofthe fundamental property of the lattice unit cell area thatwill be useful below.To prove that the shortest vectors r pm and r qn forma basis, it is sufficient to show that they define a paral-lelogram of a minimum area corresponding, in terms ofthe areal density, to exactly one lattice point, A = r a .Suppose the latter were not true, then the parallelogramspanned by r pm and r qn , besides the lattice points atthe vertices, would have contained an inner lattice point.We arrive at a contradiction by noting that the short-est distance from any inner point of a parallelogram toits vertices is smaller than one of the parallelogram sides | r pm | or | r qn | .As a side remark, the shortest vectors provide whatmay be viewed as a natural basis. Because of the shortestvector property the parallelogram spanned by r pm or r qn has a minimum size compared to those for other pairs ofprimitive vectors. As a quick reminder, the choice of abasis in a lattice is not unique, and any pair of primitivevectors can serve as a basis. For example, the lattice (3)has a basis defined by primitive vectors e = a i , e = − r x i + r y j . While this is a legitimate basis for this lattice, the lengthsof the vectors e , can greatly exceed the distance be-tween the nearest lattice points in the ( u, v ) metric. Tothe contrary, the basis formed by the shortest vectors r pm and r qn would be optimal (and thus natural) in thesense of this metric, in loose analogy with the role of theWigner-Seitz cell in studying crystal lattices.Next, we establish the relation between the vectors r pm and r qn coordinates and the parastichy numbers, givenby P = | m | and P = | n | (up to a permutation). We firstanalyze the family of parastichy lines obtained with r pm .Consider the domain in the ( u, v ) plane (2) defined by 0 ≤ u < a and 0 ≤ v < v pm = r m . This is a rectangle of area mA which contains exactly | m | lattice points including . Because these points are not nearest neighbors of eachother in the sense of r pm , any parastichy line passingthrough one of them does not pass through any otherof these points. On the other hand, two opposite sidesof the parallelogram are related by a translation by r pm ,which means that each parastichy must pass through oneof these points. Therefore, the number of parastichiesgenerated by r pm is | m | . By the same argument, thenumber of parastichies in the other family, generated by r qn , equals | n | .Since the shortest vectors r pm and r qn are defined upto a sign, we can assume, without loss of generality, thatboth m and n are positive integers. With this convention,used thoughout the article, the parastichy numbers arejust P = m and P = n . B. The space of cylindrical lattices
Here we consider the relation between the parastichynumbers P , P and the parameters A , x , y of cylindri-cal lattices. This relation, as will shortly become clear,is central to understanding the geometry of phyllotacticpatterns. Our analysis will demonstrate a nontrivial rela-tion between the quantities P , and the hyperbolic planeparameterized by the complex variable z = x + iy . Thenumbers P , will be shown to be encoded in a “hyper-bolic wallpaper”, consisting of domains in the hyperbolicplane that are mapped on each other by modular trans-formations. This connection will provide a tool that willhelp us to explain Fibonacci phyllotaxis and understandits stability.The cylindrical lattices (3) are completely defined, upto a rescaling factor √ A , by specifying the parameters −∞ < x < ∞ and 0 < y < ∞ . Hence the upper half-plane y > x, y ) plane can serve as the phasespace of all such lattices. The rescaling factor will beirrelevant for most of our discussion. Thus, unless statedotherwise, we focus on the lattices of unit density, A = 1.The pairs ( x, y ) are in one-to-one correspondence withthe planar lattices (3), since any such pair defines a lat-tice and vice versa. For the lattices on a cylinder, how-ever, all pairs ( x + n, y ) with integer n correspond to thesame lattice, because changing x by an integer amountsto changing d in (1) by a multiple of 2 π , which obvi-ously does not change the generating spiral. Hence x can be chosen in the interval 0 ≤ x <
1. Furthermore,the transformation x → − x maps a cylindrical lat-tice with one helicity to an equivalent lattice with theopposite helicity, right-hand to left-hand or vice versa.Therefore, all non-equivalent lattices can be parameter-ized by 0 . ≤ x < ≤ x < .
5) upto interchanging the right-hand and the left-hand helic- ity. Below we adopt this convention in all drawings andfigures. Later, however, when discussing the energy of alattice in Secs. IV, V, VI, it will be more convenient toconsider the entire halfplane y > x, y ) plane,because the transformations x → n ± x are merely a sub-group of a symmetry group of the energy function.We now proceed to discuss how the ( x, y ) plane is orga-nized by the parastichy pairs P , P . Given a particularparastichy pair, P = m and P = n , what is the set of x and y for which it is realized? We will call this seta parastichy domain in the ( x, y ) plane corresponding tothe parastichy pair ( m, n ). (The scale-invariant defini-tion of the parastichies makes the rescaling factor √ A anirrelevant parameter.) To find the shape of the paras-tichy domains, let us assume that the pair of shortestvectors r pm and r qn is fixed, and study how this restricts x and y .One requirement for r pm and r qn comes from the factthat the parallelogram spanned by these vectors musthave area A , since it is a unit cell of the lattice (seeabove). Therefore we must have | r pm × r qn | = A . Thiscondition is equivalent to pn − qm = ± , (7)which restricts possible combinations of p , m , q , and n .First, it follows from (7) that the integers m and n aremutually prime. Also, if m and n are positive and notequal to 1 simultaneously, the integers p and q are eitherboth positive, or both negative.Given the m and n values, what can be said about p and q ? Different pairs ( p, q ) corresponding to the same( m, n ), satisfying condition (7), are related by the trans-formation p ′ = p + km , q ′ = q + kn , where k is an integer.According to Eq.(3), the change ( p, q ) → ( p ′ , q ′ ) can alsobe realized by a transformation of the lattice parame-ter x ′ = x + k . Combining it with the transformation x ′ = 1 − x , changing helicity, we see that the pair ( m, n )uniquely determines p and q for x within the interval0 . ≤ x < x and y are constrained bythe condition that r pm and r qn are the shortest vectors.Since each of the diagonals of the parallelogram spannedby r pm and r qn must in this case be longer than its sides,we can write four inequalities:(a) | r pm ± r qn | > | r pm | ; (b) | r pm ± r qn | > | r qn | . (8)To see what these conditions mean for x and y let us writethem explicitly using the form (3) of the lattice vectors.The condition (a), taken with both the plus and minussigns, restricts the pair ( x, y ) to be in the region( n + 2 m ) (cid:0) ( x − x )( x − x ) + y (cid:1) > , ( n − m ) (cid:0) ( x − x )( x − x ) + y (cid:1) > , (9)where x = q + 2 pn + 2 m , x = qn , x = q − pn − m The condition (b) takes a similar form. Taken together,(a) and (b) define a curvilinear domain in the ( x, y ) plane,0
FIG. 3: The x, y plane partitioned by parastichy domains(thick lines). The domains are marked by parastichy pairs( P , P ). Only the 0 . < x < x → x ± , x → − x . Fareytriangles are shown by thin lines (see Fig. 5). Inset:
Aparastichy domain bounded by four semicircles (9). as shown in the inset of Fig. 3. The domain boundariesare arcs of semicircles with the diameters on the x axis.For any x and y inside the domain the lattice (3) is de-scribed by the parastichy pair P = m , P = n .Similar domains can be obtained for all combinationsof p , m , q , and n , which satisfy the condition (7). Thesedomains partition the ( x, y ) plane as shown in Fig. 3. Inthis figure we display the parastichy domains only in theregion 0 . < x < x → x + k and reflections x → k − x . This invariance property accounts for the factthat the pairs ( x, y ), ( x + 1 , y ), and ( − x, y ) represent thesame lattice up to helicity change.There is a useful relation of the plane partition by theparastichy domains P = m , P = n and van Iterson dia-gram (Fig. 2). The boundaries of the parastichy domainscorrespond to rhombic lattices, because at the border-line where the parastichy pair changes the lattice musthave next-shortest vectors of equal length. Therefore,the points where the domains join in three correspond toperfect triangular lattices (this property wil be discussedfurther in Sec. V B). At the same time, these points arethe triple branching points of van Iterson tree. Thus eachbranch of the tree resides within a single parastichy do-main, connecting its opposite corners (compare Fig. 2 andFig. 3). After associating each van Iterson branch with aparastichy domain we see that the parastichy numbers invan Iterson diagram are in fact identical to those intro-duced above using the plane partition. Therefore, eachbranch of van Iterson diagram is characterized by paras-tichy numbers P , P which are constant within it andcan change only at the branching of the tree. It is interesting to understand the organization of theparastichy numbers P , P throughout the x, y plane inFig. 3. Inspecting the parastichy pairs written in the ad-jacent domains shown Fig. 3 inset gives a simple addi-tion rule: For a domain with P = m and P = n , thetwo domains adjacent to it from below (closer to the x axis) are lebeled by the pairs P = m + n, P = n and P = m, P = m + n . As illustrated in Fig. 3, this re-lation between parastichy pairs of neighboring domains,so far conjectured empirically, describes organization ofpairs P , P in the entire x, y plane. The addition rulefor the numbers P , P will be proven in Sec. IV B.For now taking the addition property for granted, wecan infer that for any pair of mutually prime integers m > n > IV. THE ENERGY MODELA. Phyllotactic growth under stress: the sequenceof deformations and adjustments
To describe interaction between different units of aphyllotactic pattern, we employ a repulsive potential U ( r ) of a generic form [16]. In our model, this poten-tial defines forces by which the points constituing thecylindrical lattice (1) repel each other. From a micro-scopic point of view, the interaction U models the effectsof rigidity of structural units and of contact pressure be-tween them, as well as other similar effects arising atshort distances (e.g. the volume constraint in dense pack-ing of scales of a pine-tree cone). While the interactionsof a short range or hard core type perhaps would be themost relevant for systems of interest, it is beneficial toallow for a generic repulsive interaction, since this helpsto assess the robustness of phyllotactic growth.Next, we introduce the interaction potential describ-ing forces between the lattice points and define the en-ergy functional. To that end, we will focus on the centralforce potential model, in which the force F = −∇ U ( r )points along the line connecting two interacting points.The potential U ( r ) in this case is a function of the dis-tance measured using the Euclidean metric in the planeobtained by unrolling the cylinder and the cylindricallattice into 2D, as discussed above. To define the en-ergy functional for a lattice we first note that the total1energy obtained from all pairwise interactions, given bythe expression E total = 12 X pmp ′ m ′ U ( | r pm − r p ′ m ′ | ) , (10)is formally divergent for an infinite lattice, either or acylinder or in 2D. It is therefore more natural to considerthe ‘energy density’ defined as energy per one lattice site.This quantity is given by[27] E ( x, y ) = X pm U ( r pm ) , r pm = | r pm | , (11)where the sum runs over all vectors of the lattice (3).The quantity (11) is equal to the energy density per unitarea times the unit cell area A . Hereafter we suppressthe dependence on A , focusing mainly on the constantdensity case, A = 1. We note parenthetically that the r = 0 term can be taken out from the sum (11) with noimpact on the discussion. Indeed, this term correspondsto self-interaction, and thus eliminating it would merelychange the energy by a constant which is independent ofthe lattice geometry and thus inessential for our analysis.Below we will focus on the case of repulsive inter-actions, dU ( r ) /dr <
0. In addition, for mathematicalconvenience, we assume that that U ( r ) decays rapidlyenough to assure convergence of the sum (11). The func-tional form of the interaction U ( r ) can be taken, for ex-ample, as an expeonential e − r/r , a gaussian e − r /r , ora power law r − γ . It will be clear from our discussion thatthe qualitative behavior is independent of the particularform of interaction.One caveat associated with our definition of energy,Eq.(11), is that this quantity is defined as a sum over all lattice points (3), of which the points r = u i + v j with equal v coordinates correspond to the same pointof the spiral structure on the cylinder (1). Despite this,we treat these points as distinct in the sum over i and j .Such a choice of the expression for the energy, Eq.(11), isdeliberate: as we demonstrate below, the quantity (11)has hidden symmetries which underpin the robustnessand universality of phyllotactic growth.Further, one can argue that the approximation madein replacing the energy of a cylindrical lattice by that ofa lattice in the two-dimensional plane is better than itsounds. Indeed, comparing the energy (11) to the en-ergy of the cylindrical structure (1), with the interac-tion defined using the shortest distance on a cylinder, wenote that the difference of the two energies correspondsto the part of the sum (11) with the u component ofthe radius vector exceeding half of the cylinder circum-ference, | u | > a/
2. However, as we will see below, thedevelopment under uniaxial stress makes the circumfer-ence a increase so that already at a relatively early stageof growth it can exceed the range of interaction set by U ( r ). As soon as this happened, the inaccuracy in theexpression (11) becomes insignificant. This observationjustifies using the expression (11) instead of a marginally more accurate but less symmetric expression for cylindri-cal lattices (1).Next, we proceed to analyze the transformations ofcylindrical lattice structures under stress. Increasing thepressure acting along the cylinder axis has two distincteffects. One is a compression of the system in the z -axisdirection. Another is a buildup of pressure accompaniedby a density change in the system. For simplicity, belowwe will treat these two effects as decoupled, assumingthat the variation of stress does not affect pressure inthe system, and hence the lattice density ρ = A − isconstant. This approximation is made on the groundsthat the system compressibility depends mostly on thelattice average density and not as much on the detailssuch as the angles between nearest-neighbor bonds andother geometric details [24]. We will therefore treat thedistinction between the conditions of constant pressureand constant density as inessential, ignoring it in ourdiscussion.We implement the constant density approximation bydescribing the lattice defomration through changing thelattice parameter y and, at the same time, maintainingconstant A . As y varies gradually from higher to smallervalues, for each given y value we have to adjust x for thelattice energy to attain a local minimum, ∂E∂x = 0 , ∂ E∂x > , A, y = const (12)and then to examine the evolution of optimal x as y de-creases from infinity to zero. The process (12) makes x an implicit function of y . One can say that the func-tion x ( y ) describes the progression (or, “history”) of thedeformation. That is, we assume that y is controlled ex-ternally, and x is a free parameter in which the systemis trying to reach local equilibrium. The reason the rolesof x and y are different becomes obvious from their geo-metric meaning: y controls the spacing of lattice pointsalong the cylinder axis, while x corresponds to a latticeshear, or to a twist of the cylindrical structure (since A is constant, neither x nor y affect the lattice density).The energy functional E ( x, y ) governs the deformationand shear developing in the compressed lattice. The re-sulting trajectory in the lattice phase space ( x, y ) canbe linked to the behavior of the energy minima, Eq.(12).Namely, as the parameter y is becoming smaller underuniaxial compression, the energy E ( x, y ), taken at a fixed y value, acquires more and more extremal points in x . Inthe phyllotaxis problem we are interested in tracing outthe (local) minimum towards which the structure evolveswithout jumping to other minima. As we will see, it isinstructive to analyse all local minima on equal footing.Such analysis will provide, in particular, an insight intothe structure of the patterns resulting from the minimaevolution under variation of y . This is exemplified inFig. 4 where the trajectories for different energy minimaobtained for the interaction U ( r ) = e − r are shown. Theresulting pattern of trajectories has a number of interest-ing properties:2 FIG. 4: The dependence x vs. y (thick lines) for the localminima of the energy (10), describing cylindrical lattice (3)evolution under increasing uniaxial stress. The interaction U ( r ) = e − r and unit density, A = 1, was used in calcula-tion. Farey triangles (Fig. 5) are shown to help exhibit thestructure of the pattern (thin lines). Inset:
The trajectoryof the global minimum of the energy (10) switches erraticallybetween different branches of the local minima shown in themain panel.
1. For y large enough there is only one energyminimum, located at x = 0 .
5, which correspondsto the angle α in the generating spiral (1) equalto 180 o . In terms of cylindrical lattices, x = 0 . y ≈ .
25 this structure becomes unstable, andacquires a twist with a right-hand or left-handhelicity. The two opposite-helicity states arerelated by mirror symmetry, and thus have equalenergies. At this y value a single energy minimumtransforms into a pair of adjacent minima. This ismanifest in the x = 0 . y ≈ .
25, and can be associated with a bifurcationof the evolving structure.3. As y is lowered further, trajectories in Fig. 4display no other branching or bifurcation belowthe first bifurcation encountered at y ≈ . y decreases, however each of these minima appears some distance away fromthe existing minima.4. There are two principal trajectories that emerge atthe bifurcation. Because of what has been said, theprincipal trajectories display smooth evolution be-low the bifurcation, displaying no secondary bifur-cations. Hence, provided that the system evolutiongoes through the bifurcation point, at subsequentcompression it will follow one of the two principaltrajectories.Of course, while the detailed behavior of the trajecto-ries depends somewhat on the choice of the interaction U ( r ), the general behavior of the minima appears to berobust. In particular, the properties 1,2,3,4 highlightedabove are insensitive to the particular form of the inter-action. This is exemplified, for instance, by the analysisin Refs. [15] and [16], which focusses on the interactions U ( r ) = exp (cid:0) − r (cid:1) and U ( r ) = − ln r , respectively.Next we proceed to discuss what numbers of spirals, i.e.the parasticy numbers, can be realized for the structuresobtained for different trajectories. This can be analyzedmost easily by superimposing the trajectories with theparastichy domains constructed above in Sec.III B Aftersuperimposing the trajectories with the domains in Fig.3,one makes a striking observation that the principal trajec-tory, following the bifurcation, passes exclusively throughthe Fibonacci domains. This result means that, startingfrom the simplest structure, the system does not haveany other choice but to evolve into a Fibonacci struc-ture. The system, responding to y gradually varying bymaintaining local equilibrium in x , is driven through theFibonacci sequence of states. This is so because in orderto switch to a non-Fibonacci state it has to overcome afinite energy barrier.We note that in our approach to modeling phyllotaxisby a requirement that the system traces local energy min-ima the word “local” is absolutely crucial. If instead onewould have chosen to analyze the global energy minimum,for each y seeking the x value yielding the lowest energystate, the result would have been quite different. Thisis illustrated in the inset of Fig. 4 which shows the po-sition of the global minimum which jumps in a perfectlyerratic way between the principal and non-principal tra-jectories. Such jumping trajectory shows no regularitywhatsoever, in particular it passes through many non-Fibonacci parastichy domains and generates structures ofa mostly non-Fibonacci type. This obviously means thatthe global energy criterion, while being useful in a varietyof other problems, does not provide a good guidance inthe phyllotaxis problem. One might argue that this is ina sense natural since the phyllotactic systems of interest,which develop into Fibonacci structures, are macroscopiceven at an early stage of growth. Indeed, even an em-bryo at an early developmental stage is a macroscopicobject which is unlikely to transcend structural barri-ers due to thermal or environmental fluctuations, even if3that would take it to a lower energy state. This is obvi-ously related to the fact that for a macroscopic systemthe time required for reaching the lowest energy state byjumping over barriers would be very large, presumablymuch larger than the growth time. We therefore concludethat, while local energy minimization has a clear phys-ical significance in describing growth, the models basedon a global energy minimization, which are blind to thepresence of barriers, are of limited utility.Returning to the discussion of diffrent trajectories inFig. 4, it is intresting to note that the next principal tra-jectory in Fig. 4 generates structures from the Lucas se-quence 1, 3, 4, 7, 11, 18, etc. This is in good agreementwith (and provides an explanation for) the well knownfact that the Lucas numbers are the most common excep-tion in phyllotaxis. To obtain these numbers through ourmechanical development model one simply has to assumethat the system makes one mistake at the very beginningby jumping over the energy barrier to the lattice in the(3 ,
1) domain, after which it strictly follows the rules ofthe game. One can crudely estimate the barrier height,and see that it increases inversely with y . This impliesthat a mistake, if happened at all, would be most likelyto occur at higher y values, i.e. at the beginning of thegrowth. In other words, the Lucas sequence is associatedwith the mistake in the growth which is the most likelyone to occur.We verified, by performing numerical simulations andotherwise, that the results described above show con-siderable robustness and are not interaction-specific.Namely, we find that all ‘reasonable’ repulsive interac-tions, dU/dr <
0, fit the bill. One might argue thata particular form of the repulsive interaction would notmatter as long as it renders the lattice stable. This con-jecture is indeed true, as will be discussed in the next sec-tion where we show that the energy model can be treatedanalytically and rigorously. After describing the rigorousresults, we will return to the robustness property andformulate more precisely the conditions on the potentialunder which the energy model leads to Fibonacci struc-tures. The robustness property provides a lot of freedomin varying the form of interaction. Furthermore, one cangeneralize the results to the interactions that vary duringthe growth (see Sec. VI A), which provides insight intophyllotactic growth in the geometries other than cyclin-drical (see Sec. VI B).Finally, we make a cautionary remark that using thedegree of compression y as a control parameter may notbe entirely physical. In a real system, the external forcesproducing stresses in a growing system correspond topressure and strain, with the resulting deformation gov-erned by the conditions of mechanical stability. However,it can be shown [24] that in our problem the stress andthe deformation are in a one-to-one relation, and so, toavoid unnecessary complications, in what follows we willreplace the actual external forces by the quantity y de-fined above, which will act as a control parameter . B. Symmetries of the energy
In this section we will discuss the symmetry proper-ties of the energy function, Eq.(11), focusing on the to-pography defined by E ( x, y, A ) in the x − y plane. Aswe will see, the behavior of E ( x, y, A ) is rather peculiar:there are infinitely many energy minima which are alldegenerate, i.e. correspond to identical energy values.Furthermore, the E ( x, y, A ) topography is such that theminima are organized in an intricate network resemblinga mountain range with a system of valleys surroundedby peaks and passes. To understand the resulting struc-ture, we will introduce a group of modular symmetries,comprised of the transformations of the x − y plane thatleave the function E ( x, y, A ) invariant. We will choosethe fundamental domain of the symmetry group and de-fine the corresponding partition of the x − y plane suchthat it elucidates the network induced by the E ( x, y, A )topography. As we will see, this construction greatly fa-cilitates the analysis and leads to a simple explanationof the pattern of the growth trajectories, such as thatin Fig. 4. In subsequent sections the modular symmetrywill be used to treat the problem analytically and rigor-ously, and to prove the stability of Fibonacci structuresin a rather general way.To make our discussion of modular symmetry moretransparent, it will be convenient to introduce the com-plex variable z = x + iy and view the ( x, y ) plane asa complex z plane. The energy function E ( x, y, A ) willthen be defined on the so-called modular space, allowingour analysis to benefit from a high symmetry revealed bysuch a construction. Accordingly, we will use the nota-tion E ( z, A ) instead of E ( x, y, A ).Furthermore, it is also convenient to use complexparameterization for the ‘physical space’ (the unfoldedcylinder) where the lattices (3), obtained by unrollingcylindrical lattices, are defined. In passing to the com-plex notation we identify i → j → i , after which Eq.(3)reads r pm ( z ) = u pm ( z ) + iv pm ( z ) = s A Im( z ) ( p − mz ) , (13)where z = x − iy .General modular transformations of a complex planeare defined in a standard way as fractional linear trans-formations of a complex variable. These trasformationsand the associated mappings of the complex plane areoften encountered in the complex-variable calculus andits applications [25]. In particular, these transformationsplay an important role in the hyperbolic geometry. Thisrelation, as will become clear shortly, will be pivotal forour discussion.Integer modular transformations of a complex vari-able [25] are defined by using a unimodular 2 × A = (cid:18) a bc d (cid:19) , | det A | = 1 . (14)4Then the corresponding modular transformations of z aredefined as z → z ′ = ( az + b ) / ( cz + d ) , if det A = 1 ,z → z ′ = ( az + b ) / ( cz + d ) , if det A = − . (15)Separate definitions for the two signs of the determinantdet A are required to assure that the halfplane Im z ′ > A = 1, and an anti-analytic function for det A = − z → z ′ → z ′′ is a modulartransformation associated with a matrix A ′ A , where A and A describe the transformations z → z ′ and z ′ → z ′′ ,respectively.In group theory, the 2 × SL (2 , Z ).Here we are interested in a bigger group known as GL (2 , Z ) comprised of integer matrices with the deter-minant +1 or − SL (2 , Z ) asa subgroup, GL (2 , Z ) /SL (2 , Z ) = Z .The significance of these transformations is elucidatedby the following Theorem:
The energy E ( z, A ) is invariant under anymodular transformation (15). Proof:
As a first step, we show that under the trans-formations of z given in Eq.(15) the lattices (13) changein avery simple way. Namely, these transformations de-fine Euclidean rotations of the lattice, when det A =1,and a rotation combined with a mirror reflection, whendet A = − r pm ( z ) = e iφ r p ′ m ′ ( z ′ ) for det A = 1 ,r pm ( z ) = e iφ r p ′ m ′ ( z ′ ) for det A = − , (16)where p ′ , m ′ and the rotation angle φ are defined by p ′ = ap + bm, m ′ = cp + dm, exp(2 iφ ) = ( cz + d ) / ( cz + d ) . (17)Eq.(16) means that the lattices (13), under the transfor-mations (15), are mapped to isometric lattices.Next, to complete the proof, we note that the latticeenergy E ( z, A ), Eq.(11), is given in terms of an isotropic,angle-independent interaction U ( r ). As a result, thequantity E ( z, A ) is invariant under distance-preservingtransformations of the lattice, such as the rotations andmirror-reflections in Eq.(16). QED
Symmetry can be used to gain insight into the prop-erties of energy E ( z, A ) in a much the same way as,for example, periodic functions are analyzed by studying them within the fundamental period of the translationsymmetry group and then using the periodicity propertyto extend the function to the entire space. To applythis strategy we need to identify a suitable tesselation ofthe complex z plane induced by the GL (2 , Z ) symmetry.This tesselation will play a role in our analysis analogousto periodic tesselations of space used for describing thestructure of periodic functions.An appropriate group-theoretic vehicle used to charac-terize functions invariant under some symmetry group isthat of a fundamental domain . The fundamental domainis defined a geometric shape in the space of the action ofthe group such that its images obtained by applying allgroup elements cover the whole space with no gaps andno overlaps (except, possibly, over boundaries). Looselyspeaking, the fundamental domain is a minimal spaceregion that, after having been replicated by the grouptransformations, partitions the entire space.For the symmetry group (15), we need a domain in the z plane chosen so that the mappings (15) of the domaincover the z plane, overlapping only along boundaries.The fundamental domain construction for the modulargroup, introduced by Gauss, is widely discussed in math-ematical literature. (For example, see Ref. [25], Chap. 2,and Ref. [26], Chap. 5.) However, the standard Gaussianfundamental domain known as the “modular figure” willnot be that useful, because for our purpose it is too small.Instead, we will use a 3 times larger domain, which is nota truly fundamental domain in the sense of the standarddefinition. We will see below that the larger domain ismore natural from the point of view of the energy topog-raphy in the z plane. The relation of our domain to theGaussian domain will be discussed in Sec. V B.Let us begin with introducing in the halfplane Im z ≥
0a family of semicircles: y = r(cid:16) x − pm (cid:17) (cid:16) qn − x (cid:17) , pm ≤ x ≤ qn , (18)where p , m , q , and n are integers with pn − mq = ± m > n >
0. (In Figs. 3, 4 the semicircles are shownby thin lines.) We will denote each semicircle by its endpoints: [ p/m, q/n ]. It is convenient to allow formally thecombinations in which either m or n is zero, by addingvertical straight lines x = q . (It follows from pn − mq = ± m = 0, then n = 1.) Our notation for such“generalized semicircles” that connect x = q with x = ∞ ,will be [1 / , q/ x − y plane into curvilinear triangles with thevertices on the x axis at the points x = p/m, q/n, r/s ,where r = p + q , s = m + n . We will denote suchtriangles by specifying the x -coordinate of their ver-tices: [ p/m, q/n, ( p + q ) / ( m + n )]. These triangles arecalled Farey triangles , and are closely related to the socalled Farey numbers. Farey triangles and Farey num-bers constitute a nice subject on the borderline between5
FIG. 5: Relation between the Farey triangles (Fig. 3, thicklines) and parastichy domains (Figs. 3,4, thin lines), whichpartition the complex z plane in two different ways, is illus-trated. Each Farey triangle is centered at a corner of threeadjacent parastichy domains. Each parastichy domain ovelapswith two adjacent Farey triangles. arithmetic and elementary geometry, and are well ac-counted in mathematical literature (see, e.g., the text-books [25, 26].)In the number theory, Farey numbers is a name fora construction that organizes all rational numbers 0 ≤ p/m ≤ / /
1, one applies the Farey sum rule, p/m ⊕ q/n =( p + q ) / ( m + n ), and successively generates more and morerational numbers (see [25], Chap.5). The order in whichthe numbers are generated coincides with the hierarchyof our Farey triangles: one vertex of each triangle is Fareysum of two other vertices (see Fig. 5).Let us comment on the relation between the Farey tri-angles and the parastichy domains. Both partitions ofthe z plane are invariant under modular transformations.Comparing Eq. 18 with the parastichy domains construc-tion of Sec. III B, it is obvious that for a given Fareytriangle [ p/m, q/n, ( p + q ) / ( m + n )], the sides [ p/m, q/n ],[ p/m, ( p + q ) / ( m + n )], and [ q/n, ( p + q ) / ( m + n )] belong tothe parastichy domains with the parastichy pairs ( m, n ),( m, m + n ), and ( n, m + n ), respectively (see Eq. 9 andFig. 5). Hence, each Farey triangle overlaps with threeparastichy domains, one side per one domain, and con-versely, each parastichy domain overlaps with two Fareytriangles. The relation of the two ways of partitioningthe z plane will be studied in more detail in Sec. V, andthen used for the discussion of the stability of Fibonaccinumbers.To see the role of Farey triangles, in Fig. 6 we draw acontour plot of the energy for the potential U ( r ) = e − r ,together with the Farey triangles partitioning the z plane.It is evident from the figure that the energy contours dis-play similar behavior in each triangle. Qualitatively, fromthe topography point of view, there is a “valley” inside FIG. 6: A contour plot (thick lines) of the lattice energy forthe interaction U ( r ) = exp( − r ). Farey triangles (thin lines)are shown for guidance. The contours behavior in all Fareytriangles is similar. each Farey triangle, with a minimum at the center. Foreach valley there are three neighboring valleys in the ad-jacent triangles, “connected” with it by “passes” throughthe saddles located at the middle of each side of the Fareytriangles. At the corners of the triangle there are peaksof infinite height.Schematically, the valleys in Fig. 6 are connected withthe neighboring valleys in a structure that can be charac-terized as a Cayley tree with a branching number three.This tree can be compared to the pattern of trajecto-ries in Fig. 4 obtained by looking at the energy minimain x at fixed y . Obviously, no matter what y , the min-ima will reside somewhere within the valleys, away fromthe peaks. Therefore, when y varies continuously, theminima trace out the valleys and their interconnections.Comparing Fig. 4 to Fig. 6, it is evident that the patternof the minima approximately repeats the pattern of thevalleys. From a topological point of view, the branchingpoints of the Cayley tree of valleys correspond to quasi-branching of the minima trajectories that one defines bybridging different trajectories across the gaps [16]. As y varies, each trajectory explores a sequence of neighbor-ing valleys, and quasibranchings manifest that there arethree neighboring valleys around each valley.However, let us emphasize that the trajectories of min-ima in Fig. 4 are not invariant under the modular trans-formations. Technically speaking, the reason is that weare looking for minima of a modular symmetric function E ( z, A ) under a non-symmetric constraint: y = const.Obviously, if the pattern of minima were symmetric,all parastichy pairs would be equally likely to appear.Hence, the absence of symmetry is an ultimate cause ofthe appearance of Fibonacci numbers.6It is important to realize that the problem has in-complete symmetry: the lattice energy is modular-symmetric, but the deformation process is not. Thusa concept of asymmetry emerges, which implies the ab-sence of symmetry, but nevertheless certain closeness toan exact symmetry. We will see that the asymmetry isa much more powerful notion than just total absence ofsymmetry. The modular symmetry transformations en-able one to compare the energy minima trajectories indifferent Farey triangles. It is even correct to say thatthe problem of stability of Fibonacci numbers in phyl-lotaxis is now reduced to the analysis of the asymmetryof the pattern of energy minima trajectories. This ideawill be the bottom line of the following discussion. V. UNIVERSALITY OF FIBONACCI NUMBERSA. Plan of the discussion
In this section we will prove a rather general statementabout any energy minima trajectory. Although, strictlyspeaking, we are interested only in the principal trajec-tory (which is Fibonacci), there is a lot of advantage inconsidering all trajectories together.To formulate our main result, let us recall that a gen-eralized Fibonacci sequence is a sequence of integers Φ n obeying the Fibonacci recursion relation,Φ n +1 = Φ n + Φ n − , (19)where the starting numbers Φ , Φ can be any integers.For example, the standard Fibonacci sequence and theLucas sequence are given by Φ = 1 , Φ = 1 and Φ =1 , Φ = 3, respectively. Main theorem:
For any interaction U ( r ) from thefamily defined below, a continuous trajectory of the en-ergy minima in the z plane goes through the sequenceof parastichy domains with the parastichy pairs from ageneralized Fibonacci sequence. The first two numbersof the sequence are determined by the beginning of thetrajectory.Once the theorem is proven, the result about havingonly Fibonacci numbers on the principal trajectory fol-lows from the fact that it begins in the (1 ,
1) parastichydomain.It is worth remarking that in the theorem the word“continuous” is crucial! Given that the deformation ofthe lattice is continuous, the theorem states that the de-formation stages all are (generalized) Fibonacci. On theother hand, if the deformation makes the lattice unsta-ble, and it abruptly transforms to some other structureor lattice, with a jump on the z plane, the theorem is notclaiming anything.The class of potentials U ( r ) for which the theoremholds is defined by the following two conditions. a) The deformation is a continuous process, and does notmake the lattice unstable. b) The topography of the energy E ( z ) is the simplestrequired by modular symmetry: no critical points otherthan the △ and (cid:3) points (see Fig. 6 and Sec. V B).Condition a) is stated in the theorem, and the role andmeaning of condition b) will become evident in the proof(see Sec. V C). The conditions a) and b) define implicitlya large family of interactions, apparently including all re-pulsive potentials. At the moment, however, an explicitcharacterization of this family is lacking. (Maybe it isworthwhile to mention that we failed to find exceptionsamong repulsive interactions.)The proof of the theorem relies on the modular sym-metry of the energy E ( z ) (see Sec. IV B), which givesthe topography of E ( z ) in the entire z plane from thatwithin a single Farey triangle. Also, we will use the re-lation between Farey triangles and parastichy domains,already mentioned in Sec. III B. We will discuss it againand summarize in Sec. V B. The plan of the proof in-cludes the following two steps. First, by making use ofthe modular symmetry, we replace the global statementabout the trajectories behavior in the entire z plane by anequivalent local statement about the behavior in a singleFarey triangle (see Lemmas 1, 2, 3 in Sec. V C). Then,by using modular symmetries, we treat a minima trajec-tory within a single Farey triangle, and find the relationwith the parastichy domains. For that, we map the lines y = const in an arbitrary Farey triangle onto a certainfamily of curves in a “reference” Farey triangle, and then(locally) minimize the energy on these curves. At thisstep, we establish a Fibonacci type relation between theparastichy domains traced by a trajectory.Loosely speaking, the Farey triangles play an “organiz-ing role” in the z plane. They represent something likestandard blocks, or units for the trajectories. In somesense, the behavior of trajectories inside all triangles issimilar, and the topology of the trajectories pattern canbe constructed by replicating one triangle with the mod-ular symmetries. B. Modular symmetry and Farey triangles
Here we review basic facts about Farey triangles onthe hyperbolic plane, and study the relation of the Fareytriangles and parastichy domains. Necessarily, our dis-cussion be brief. (We refer interested reader to very goodtexts by Apostol [25] and Iversen [26].)The geometric role of modular transformations is thatthey preserve the hyperbolic metric dl = | Im z | − d ¯ zdz .The Im z ≥ hyperbolic plane . The transformations thatpreserve the hyperbolic metric (they are called isome-tries ) form a group known as P SL ( R ) or SL ( R ) ⊗ Z ,that plays a role in the hyperbolic geometry similar tothat of space translations and rotations in the Euclideangeometry. All hyperbolic isometries have the form (15)with arbitrary real a , b , c , and d , such that det A = ± a , b , c ,7 d , and det A = ±
1, form an infinite discrete subgroupof the group of all isometries, as, for instance, symme-tries of a crystal in the Euclidean space do. The grouptheory notation for the modular group is
P SL ( Z ) or SL ( Z ) ⊗ Z .The geodesics [26] of the hyperbolic metric are semi-circles with the diameters on the real axis: y = p ( x − x ) ( x − x ) , x < x < x . (20)Obviously, the isometries (15) map any geodesic to ageodesic. There are two different classes of hyperbolicisometries: a) The mappings with det A = +1 that have just twofixed points on the real axis; b) The mappings with det A = − z ′ − x z ′ − x = − ¯ z − x ¯ z − x . (21)Such a transformation is analogous to a Euclidean reflec-tion, with the geodesic (20) corresponding to a mirror.A composition of an even number of transformations ofthe type b) is a transformation of the type a) .The property of Farey triangles that makes them usefulis that the modular transformations (15) map any Fareytriangle either onto another Farey triangle, or onto itself.To see why is that, one first notes that by a transforma-tion (15) with unrestricted a , b , c , and d , we can mapany three points on the real axis onto any other threepoints (because it is a fractional-linear function). Thenone checks that, since vertices of Farey triangles are ra-tionals, the numbers a , b , c , and d can be chosen to beinteger, with det A = ±
1. This follows from the explicitform of the transformation (15) that maps an (arbitrary)Farey triangle [ p/m, q/n, ( p + q ) / ( m + n )] onto the trian-gle [0 / , / , / A = (cid:18) m − pm − n − p + q (cid:19) . (22)Applied to the lattice energy E ( z ), modular symme-tries relate the values of E at different points of the z plane. Given E ( z ) inside one Farey triangle, the modu-lar transformations extend it throughout the whole plane.However, let us emphasize that the function E ( z ) withinone Farey triangle has certain symmetry properties, andhence it is not completely arbitrary. The reason is that,as it was mentioned above, the Farey triangles are largerthan the actual fundamental domains of the modulargroup. For each triangle there are six modular trans-formations (including the indentity) that map it ontoitself. The energy E ( z ) is invariant under these trans-formations, and hence there is a usual relation betweenthe symmetries of a triangle and the behavior of the func-tion E ( z ) inside it: a) There are extremal points of E ( z ) at all symmetrypoints of the triangle transformations; b) The contours of E ( z ) are normal to the invariant linesof the transformations.Since the topography of the energy E ( z ) within a Fareytriangle will be crucial for our discussion, let us describehere the triangle symmetries in some detail.There is a total of six modular transformations whichmap a Farey triangle onto itself. From the group theorypoint of view, the symmetry group of a Farey triangle isidentical to that of an equilateral triangle in a Euclideanplane that includes three reflections, two rotations, andthe identity. Each permutation of the triangle verticesdefines a modular transformation (since it is a fractional-linear transformation). For example, the transformationof the triangle [ p/m, q/n, ( p + q ) / ( m + n )] that inter-changes p/m with q/n , and preserves ( p + q ) / ( m + n ),is a reflection (21) with respect to the geodesic [( p + q ) / ( m + n ) , ( p − q ) / ( m − n )]. The two other reflectionscorrespond to the geodesics [ p/m, ( p + 2 q ) / ( m + 2 n )] and[ q/n, (2 p + q ) / (2 m + n )]. Like in the Euclidean plane,where a composition of two reflections is a rotation, forthe Farey triangle the composition of two reflections is atransformation giving rise to a cyclic permutation of thevertices. The intersection of the three geodesics z = pm + qn + pn/ mq/ m + mn + n (23)+ i √
32 1 m + mn + n . is a symmetry point of the triangle invariant under allsix symmetries of the triangle.The symmetry lines of the Farey triangle are thethree geodesics listed above, and also the three sidesof the triangle, [ p/m, q/n ], [ p/m, ( p + q ) / ( m + n )], and[ q/n, ( p + q ) / ( m + n )], since the reflections (21) aboutthem are modular transformations that map the triangle[ p/m, q/n, ( p + q ) / ( m + n )] to the three adjacent trian-gles. There are three symmetry points on the sides ofthe triangle given by the intersections of the symmetrylines. For example, the side [ p/m, q/n ] intersects with[( p + q ) / ( m + n ) , ( p − q ) / ( m − n )] at z = pm + qnm + n + i m + n . (24)The symmetry lines and the points (23) and (24) havea simple meaning in terms of the geometry of the lat-tice (3). The Farey triangle sides correspond to thelattices with a rectangular unit cell. The geodesics[( p + q ) / ( m + n ) , ( p − q ) / ( m − n )], [ p/m, ( p +2 q ) / ( m +2 n )],and [ q/n, (2 p + q ) / (2 m + n )] correspond to the latticeswith a rhombic unit cell. One verifies this immedi-ately by checking that the conditions r pm · r qn = 0,and | r pm | = | r qn | define the geodesics [ p/m, q/n ], and[( p − q ) / ( m − n ) , ( p + q ) / ( m + n )], respectively. Conse-quently, the points of the form (23) correspond to per-fect triangular lattices, and the points (24) correspond to8square lattices. We will call them △ points and (cid:3) points,respectively. Each Farey triangle contains one △ pointin its interior and three (cid:3) points, one per each side.The parastichy domains, as we found in Sec. III B, aredefined by the lines corresponding to rhombic lattices,i.e., by the symmetry lines. The corners of the paras-tichy domains are the △ points. The boundaries are theportions of the symmetry lines extending from the △ point to the Farey triangle vertices. At each △ pointthree parastichy domains are adjacent. The parastichypairs of the three parastichy domains overlapping witha Farey triangle [ p/m, q/n, ( p + q ) / ( m + n )] are ( m, n ),( m, m + n ), and ( n, m + n ).The symmetry points (23) and (24) are extremal pointsof the energy E ( z ). At the △ points the energy musthave local minimum, since a perfect triangular lattice isa lowest energy configuration for an isotropic repulsiveinteraction. The minima at the centers of the valleys in-side the Farey triangles in Fig. 6 are at the △ points. Atthe corners of the Farey triangles the energy is maximal,since at these points one of the lattice periods vanishes,and the energy diverges. Given that the maxima are atthe corners, and the minima at the centers, the (cid:3) pointshave to be the saddle points. (The square lattice hashigher energy than the triangular lattice, and it is unsta-ble with respect to a deformation towards triangular lat-tice through a continuous sequence of rhombic lattices.)Such a configuration of extremal points, located only atthe symmetry points, is a minimal combination requiredby symmetry. For an arbitrary interaction, in princi-ple, one could have other extremal points of the energy,not associated with the symmetry points of the Fareytriangles. However, that would make the discussion un-necessarily complicated, and moreover, for typical simplerepulsive interactions (exponential, gaussian, power law,etc) there are no additional extremal points. Therefore,in the following discussion we will assume that there areno other extremal points besides those required by sym-metry.Another useful property is that, by symmetry, the con-tours of the energy E ( z ) are normal to the Farey trianglesides and the symmetry lines, everywhere except the (cid:3) and △ points. This result will be crucial in the followingdiscussion. C. Farey partition of the hyperbolic plane and thetrajectories of energy minima
To prove the main theorem, we have to see why thewhole infinite sequence of the parastichy domains tracedby an energy minima trajectory has the parastichy pairsthat obey the Fibonacci addition rule. For that we willstudy the relation between the Farey triangles and theminima trajectories. The theorem statement is global ,since it states the trajectories behavior in the entire z plane. Our first step will be to replace it by local state-ments which are equivalent to the theorem, but talk only about one Farey triangle. One can say that we reducethe mystery all over the z plane to the mystery inside oneFarey triangle. The following three Lemmas together areequivalent to the statement of the main theorem. Lemma 1
The minima trajectories enter or leave Fareytriangles only through the (cid:3) points.
Lemma 2
The trajectories enter or leave the parastichydomains only through the △ points. Any △ point belongsto one of the trajectories. Lemma 3
At a △ point, by going from one to anotherparastichy domain the trajectory obeys the Fibonaccirule: after exiting the domain ( m, n ), m < n , it entersthe domain ( n, m + n ).Given the starting point of a continuous trajectoryof minima, the three Lemmas completely determinethrough which sequence of Farey triangles and parastichydomains the trajectory proceeds. The Lemmas 1 , statethat, no matter what interaction, there are (cid:3) and △ points on each trajectory. Although the trajectory maybe sensitive to the interaction, the sequence of the (cid:3) and △ points through which it goes is decided purely geomet-rically, and does not depend on the interaction.As for the △ points, since perfect triangular latticehas lowest energy, it is understandable that the △ pointsalways lie on trajectories. However, the fact that (some) (cid:3) points also lie on the trajectories comes as a surprise.We will see that both Lemma 1 and follow from themodular symmetry. Lemma 3 states that the behavior of trajectories nearthe △ points determines the evolution of parastichy pairs.Each △ is a corner of three parastichy domains (say, withthe pairs ( m, n ), ( m, m + n ), and ( n, m + n )). The rulethat the change of the pair is always ( m, n ) → ( n, m + n ),where m < n , applied iteratively along the trajectory,evidently leads to a generalized Fibonacci sequence. Ac-cording to Lemma 3 , there is no trajectory branching:new trajectory in the parastichy domain ( m, m + n ), m < n , emerges at a finite distance from the △ point.Note that Lemma 3 is not applicable to the triangle[0 / , / , / m = n = 1), where branching occurs atthe bifurcation point. Proof:
To prove the lemmas for a given Farey triangle(say [ p/m, q/n, ( p + q ) / ( m + n )], m < n ) we use modularsymmetry (22), and map this triangle onto [0 / , / , / y transformunder this mapping, we have to find out how the con-straint y = const transforms. Under the transformation(22), the family of lines y = λ is mapped to a family ofcircles tangent to the real axis, | z − ( x + iy ) | = y , (25)where x = m/ ( m − n ), and y = 1 / (2( m − n ) λ ). Theimages of the minima can be found by (locally) minimiz-ing the energy on the circles (25).This task is greatly facilitated by the following variablechange, z → w = ( − i √ z + ǫ ) / ( i √ z + ǫ ) , (26)9 FIG. 7: (a)
A contour plot (thick lines) of the lattice energyfor the interaction U ( r ) = exp( − r ), mapped to the w planeby (26). The triangle [ C, B, A ] is a mapping of the Fareytriangle [0 / , / , / (b) A trajectory of energy minimumenters the triangle [
A, B, C ] through the point D , and passesthe point O to exit through the point F . A new trajectory iscreated near the point E , and exits through the point E . where ǫ = e iπ/ = 1 / i √ /
2. It maps the upper z planeonto a unit circle | w | = 1. Fig. 7a shows the contours ofthe energy taken as a function of w .Under the variable change (26), the Farey trian-gle [0 / , / , /
2] maps onto the triangle [ ǫ , ǫ ,
1] =[
C, B, A ]. The △ point maps to the point O and thethree (cid:3) points map to the points D , E , and F . Theparastichy boundaries map onto the lines OA , OB , and OC . The curved triangles AOB , BOC , COA belong toparastichy domains ( n, m + n ), ( m, n ), and ( m, m + n ),respectively.The family of circles in Eq. 25 maps to a new family of circles | w − w (1 − r ) | = r,w = ( m − ( m + n ) ǫ ) / ( m − ( m + n ) ǫ ) , (27) r = √ / (cid:0) √ λ (cid:0) m + mn + n (cid:1)(cid:1) . which are tangent to the unit circle at w . We computearg( w ) = 2 arctan ( n + m ) / ( n − m )) , (28)and find 2 π/ < arg( w ) ≤ π , since 0 ≤ m ≤ n .The behavior of the energy contours within the trian-gle [ A, B, C ] follows from the results of Sec. V. The con-tours are smooth curves everywhere except the extremalpoints: A , B , C , D , E , F , and O . Also, according to thelast statement of Sec. V, the contours are everywhere per-pendicular to the arcs AB , BC , CA , and to the segments AD , BE , CF .Now we are ready to study the trajectories. Geomet-rically, the trajectories points are those where the con-tours are tangent to the circles (27). It is simple to seethat the circles (27) cannot be perpendicular to AB , BC , CA , AD , BE , CF , provided that arg( w ) = π . Butarg( w ) = π can occur only for the special case m = 1, n = 1, excluded by the condition m < n . Thus, weproved Lemma 1 and
Lemma 2 .Modular symmetry implies that in the vicinity of the △ point O the contours are circles. From that, it fol-lows that at the point O the trajectory is tangent tothe straight line joining O and w (Fig. 7(b)), and alsothat the trajectory branching does not happen. Since2 π/ < arg( w ) < π , the trajectory in BOC passes O toenter AOB . Recalling that parastichy pairs in
BOC and
AOB are ( m, n ) and ( n, m + n ), respectively, one gets theFibonacci rule, thus proving Lemma 3 . QED
Finally, let us discuss the case m = 1, n = 1, in which Lemma 3 does not hold. However, if the method of theproof is applied to this case, it yields arg( w ) = π , thatis a family of circles (27) perpendicular to AD . Then,obviously Lemma 1 is still valid, since the circles cannotbe perpendicular to AB , BC , CA . The trajectory goesalong the segment AD up to some point within the seg-ment AO , where the bifurcation occurs. At this point,the trajectory gives rise to a pair of symmetric princi-pal trajectories, which reflects the x → − x symmetry.Thereafter, the principal trajectories are described by thelemmas, and hence correspond to Fibonacci structures. VI. ROBUSTNESS OF THE PROBLEMA. Deformation versus anisotropic growth
So far, we studied cylindrical lattices under the fixeddensity condition (or, equivalently, under fixed pressure),and found a relation between phyllotaxis and the de-formation caused by a uniaxial stress applied along thecylinder axis. However, in order to understand the ubiq-uity of phyllotaxis, as discussed in Sec. II A, one has to0apply these results to the growth in different geometries:cylinder, disk, or cone. In this section we will discussimplications of the theorem proved in Sec. IV B for thesegrowth problems. We will see that, because of the ro-bustness of the problem, phyllotaxis can be understoodin each case by an appropriate modification of the defor-mation picture.Let us start with the simplest case of a cylindricalstructure that grows anisotropically: faster in the cylin-der crossection, and slower along the axis. Intuitively, upto a dimension rescaling by a factor √ ρ , where ρ is thedensity, we return to the deformation problem at fixeddensity. However, since the rescaling changes the inter-action U ( r ), and this might affect the growth, let us lookat the problem closer.The anisotropic growth of a cylindrical structure (1)can be accounted for by assuming that the density andthe interaction are some functions of y : ρ ( y ), U y ( r ), sincedecreasing y in this problem plays a role of the time.By density rescaling, the energy E can be brought tothe form in which A = ρ − = 1, and the interactionparametrically depends on y : E ( x, y ; A, U ( r )) = E ( x, y ; 1 , U eff ( r )) ,U eff ( r ) = U y ( √ Ar ) . (29)This formula replaces the density change by an effectiveinteraction evolution. To find out the effect of the latteron the energy minima trajectories, one has to go backto Lemma 1, 2, and . Obviously, as long as the effec-tive interaction in Eq. 29 satisfy the requirements listedat the end of Sec. V A, the lemmas hold: the minimatrajectories enter and exit Farey triangles through the (cid:3) points; they enter and exit the parastichy domainsthrough the △ points; and at the △ points the Fibonaccirule is obeyed. The reason is that each of the lemmas is a“local” statement with respect to y , as it can be verifiedby considering the situation at one particular y .From that, we conclude that the deformation and thegrowth problems are equivalent. By a similar argument,one can extend the result to the case when the inter-action varies during the growth in some arbitrary way,consistent with the conditions required by the theorem.Finally, there are interactions for which the deforma-tion and the growth problem are equivalent exactly, notjust topologically, or in the sense of parastichy pairs. Anypower law potential U ( r ) = U r − γ has such a property,which follows from its scale invariance. B. The disk and cone geometries
Plants with disk geometry exhibit spiral phyllotaxis(see Fig. 1(b)), which apparently is quite similar to thecylindrical phyllotaxis. Naively, one can say that a com-plex mapping z → ln z transforms one problem to theother, as it maps punched z plane to a cylinder, andtransforms log–spirals in the z plane to helices on the cylinder. However, there is one interesting feature thatdistinguishes the spiral phyllotaxis from the cylindricalone. For example, a sunflower is divided into circularconcentric domains of ring shape which have differentparastichy numbers. From one domain to another, theparastichy numbers change so that they grow outwardsthe disk, and the transition of the numbers across thedomain boundary (called “parastichy transition”) followsthe Fibonacci rule.To apply the energy model, we note that sufficiently faraway from the center, outside a core region, the structurelocally looks like a lattice. So, one can write a latticeenergy, as we did for the cylindrical problem, and studythe evolution under a stress. The origin of the stressesin this case, as we discussed in Sec. II A, is at the diskcenter, where during the growth new structure units arebeing generated. Older units are pushed by newer ones,and move outwards as the structure grows.Let us model the structure by the points of a spirallattice (Fig. 1(b)), and order the lattice points, given by( ρ i , θ i ) in cylindrical coordinates, in the order of increas-ing radius ρ i . To characterize the structure locally as aperiodic lattice, in analogy with the discussion in Sec.II B(cf. Eq.(1), we introduce two sets of parameters: r i = ρ i +1 − ρ i , a i = 2 πρ i , d i = θ i +1 − θ i , (30)and x i = − d i / π, y i = r i /a i , A i = r i a i . (31)We assume that outside a core region, i > i min , the pa-rameters are slowly varying as function of i , so that thenotion of local periodic lattice can be used. The param-eters (30) and (31) have the same geometrical meaningas the corresponding parameters (1), (3) of a cylindricallattice.The density is approximately constant throughout thedisk, outside the core region. In terms of the parameters(30) and (31), A i = A = const for i > i min . From that,outside the core region, y i decreases with the radius: y i = A/ (2 πρ i ) . Then, since A i and y i are completely fixedby the growth process, x i is the only free parameter left,and we again have a problem of optimizing x i so that thelattice energy attains a local minimum.For a short range interaction of radius ρ , the size ofthe core region ρ c ti is comparable to ρ . Inside this regionone cannot identify periodic lattice, even approximately,and the growth has to be described in a different way.However, the study of the core region lies beyond thescope of this paper, and we refer the reader to otherchapters of this volume as well as Refs. [5, 6, 17]. Belowwe assume that, due to the short range of interaction,the core has essentially no effect on the lattice structurefar away, where parastichy transitions occur. The role ofthe core is merely in setting right initial conditions forthe deformation process, which can be compared withthe selection of initial lattice in the cylindrical problem.Similarly to that problem, as we will see, once initial1structure is given, the final state of its deformation isdetermined anambiguously.Now, let us discuss the growth process. Total energyof the structure is given by E tot = X i>i c E i , E i = 12 X j = i U ( | r j − r i | ) , (32)where E i is the energy of the point i interacting withother points. For a short range interaction, E i can bereplaced by the energy E ( x i , y i ; A i ) of a periodic latticewhose structure is locally identical to that of the spirallattice, E tot = X i>i c E ( x i , y i ; A i ) . (33)For slowly varying parameters, the summation over i canbe replaced by an integration over ρ , E tot = Z ρ c dρ πρ A − ( ρ ) E ( x ( ρ ) , y ( ρ ); A ( ρ )) , (34)where 2 πρ A − ( ρ ) is the radial density. We note that E tot depends only on x ( ρ ), since y ( ρ ) and A ( ρ ) are fixed.To find a stable structure of a spiral lattice, one has tominimize E tot with respect to x ( ρ ), with a constraint that x ( ρ ) is continuous. Evidently, this minimization problemleads to the same minima trajectories as in the cylindri-cal problem, since ρ and y are in a one-to-one relation.This equivalence fixes the solution of the disk problem x ( ρ ), making the sequence of parastichy pairs in the twoproblems identical.One can say that the time of the cylindrical problembecame the radial dimension of the disk problem. Thereason is that the degree of compression y in the diskproblem decreases as ρ − at ρ ≥ ρ . Hence, the changeof a parastichy pair is represented in the disk problemby a circular boundary at which a parastichy transitionoccurs. Finally, we comment on the universality of the paras-tichy transitions. Evidently, the particular form of A ( ρ )and y ( ρ ) assumed above can be replaced basically by anyother reasonable dependence, by an argument similar tothat of Sec. VI A. The parastichy transitions are thereforea universal property of spiral lattices, independent of thedetails of the density profile. In this way our discussionof parastichy transitions can be extended to other geome-tries, e.g. the cone geometry intermediate between thoseof a disk and cylinder, supplying a general argument forthe stability of Fibonacci numbers. VII. SUMMARY
The mechanical theory explains phyllotaxis by stip-ulating that the growth of a plant at a critical stage(perhaps embryonal) is anisotropic. It argues that stressbuildup and relaxation in a deformable lattice induced byits growth is such that it leads to Fibonacci phyllotacticpatterns, exclusively and deterministically. In the me-chanical theory reviewed above this result is derived rig-orously by analyzing the mechanics of deformable cylin-drical lattices. This theory explains the predominanceof Fibonacci numbers in phyllotaxis, as well as why themost frequent exceptions in phyllotaxis are described byLucas numbers.The robustness of phyllotactic growth in other geome-tries, such as a disk or a cone, can be understood bycombining the stability if Fibonacci numbers in cylin-drical lattices with the general properties of deformablelattices growing under stress. Indeed, stress buildup andrelaxation in non-cylindrical lattices, if viewed locally,follow the same rules as in the cylindrical case. As aresult, the outcome of the anisotropic deformation is in-sensitive to the specifics of the interaction energy of thephyllotactic pattern or its geometry, be it a cylinder, adisk, or a cone. This provides further insight into theubiquity of Fibonacci numbers in natural growth. [1] D’Arcy W. Thompson,
On Growth and Form (CambridgeUniversity Press, 1942)[2] H. Weyl,
Symmetry , (Princeton University Press, Prince-ton, N.J., 1952)[3] Parastichy (from Wikipedia, the free encyclopedia)[4] R. V. Jean,
J. Theor. Biol. , , 639-660 (1978)[5] G. H. Mitchison, Science , , 270 (1977)[6] H. Meinhardt, Ber. Deutsch. Bot. Ges. , , 101 (1974)[7] R. V. Jean, Mathematical approach to patterns and formin plant growth , (Wiley and Sons, 1984);
Phyllotaxis,a systematic study in plant morphogenesis , (CambridgeUniversity Press, 1994)[8] R. F. Lyndon,
Plant development, the cellular basis , (Un-win Hyman, 1990)[9] H. Airy, Proc. R. Soc., , 176-179 (1873); Proc. R. Soc., , 299-307 (1873) [10] L. Bravais & A. Bravais, Ann. Sci. Nat. (2), , 42-110,193-221, and 291-348 (1837)[11] G. Van Iterson, Mathematische und Mikroskopisch-Anat-omische Studien ¨uber Blattstellungen (Gustav Fischer,Jena, 1907);[12] S. Schwendener, Mechanische Theorie der Blattstellun-gen, (Leipzig, Engelman, 1878)[13] I. Adler, J. Theor. Biol., , 1 (1974)[14] I. Adler, J. Theor. Biol., , 29 (1977)[15] L. S. Levitov, Phys. Rev. Lett. , 224 (1991).[16] L. S. Levitov, Europhys. Lett. , 533 (1991); and Pis’maZh. Eksp. Teor. Fiz. , 542 (1991) [JETP Lett. , 546(1991)].[17] S. Douady and Y. Couder, Phys. Rev. Lett. , , 2098(1992)[18] R. O. Erickson, Science , , 705 (1973) [19] N. Rivier, Mod. Phys. Lett. B, , 953 (1988)[20] F. Rothen and A. -J. Koch, J. Phys. France , 633-657and 1603-1621 (1989);[21] P. Prusinkievics and A. Lindenmayer, The AlgorithmicBeauty of Plants (Springer-Verlag, New York, 1990).[22] N. Rivier, R. Occelli, J. Pantaloni, and A. Lissowski, J.Phys. (Paris), , 49 (1984).[23] R. D. Meicenheimer, Am. J. Botany , 557 (1979). [24] H.-W. Lee and L. S. Levitov, unpublished[25] T. M. Apostol, Modular Functions and Dirichlet Seriesin Number Theory , Chaps. 2, 5, (Springer-Verlag, 1976)[26] B. Iversen,