The Hidden Diversity of Vascular Patterns in Flower Heads

Vascular systems are intimately related to the shape and spatial arrangement of the plant organs they support. We investigate the largely unexplored association between spiral phyllotaxis and the vascular system in Asteraceae flowers heads. We imaged heads of eight species using synchrotron-based X-ray micro-computed tomography and applied original virtual reality and haptic software to explore head vasculature in three dimensions. We then constructed a computational model to infer a plausible patterning mechanism. The vascular system in the head of the model plant Gerbera hybrida is qualitatively different from those of Bellis perennis and Helianthus annuus, characterized previously. Cirsium vulgare, Craspedia globosa, Echinacea purpurea, Echinops bannaticus, and Tanacetum vulgare represent variants of the Bellis and Helianthus systems. In each species the layout of the main strands is stereotypical, but details vary. The observed vascular patterns can be generated by a common computational model with different parameter values. In spite of the observed differences of vascular systems in heads, they may be produced by a conserved mechanism. The diversity and irregularities of vasculature stand in contrast with the relative uniformity and regularity of phyllotactic patterns, confirming that phyllotaxis in heads is not driven by the vasculature.

The simulation of vascular system development is preceded by preprocessing steps.First, receptacle contours and active ring positions measured at discrete time points during head development are interpolated by the C-model to describe the receptacle contour and the position of the active ring as continuous functions of time (Sec.2).Both of these functions are employed by the model of phyllotaxis (P-model) to generate time-stamped initial primordia positions on the growing receptacle (Sec.3).In addition, the H-model uses a sequence of receptacle contours output by the C-model to create a corresponding sequence of the entire receptacle sections, including their interior (Sec.4.2).The development of the vascular system is then simulated iteratively, in a loop consisting of two operations taking place over incrementally advancing time: • Insertion of new vascular strands (V-model) (Sec.5).Typically, these strands are induced by the addition of new primordia and connect them to the existing vasculature (initially, the procambial ring at the base of the receptacle).Consistent with our observations, the gerbera model also incorporates centripetally extending adaxial strands, which support floret veins formed later.
• Receptacle growth (G-model) (Sec. 4.3).This operation entails incremental deformation of all strands embedded in the growing receptacle, as well as the repositioning of primordia on the receptacle surface due to the growth of the receptacle over a simulation step.
The simulation process is illustrated in Fig. SMD2.

The contour expansion model (C-model)
The C-model produces a continuous-time description of the expanding receptacle contour by interpolating digitized profiles of receptacles scanned at different stages of development (Fig. S9a,d,g).For each plant species, it is constructed by following the steps described by Zhang et al. (2021): 1. Volumetric scanning of several (7-10) heads capturing head development from an early stage to the stage at which all primordia and vascular strands have been patterned; 2. Extraction of a virtual longitudinal section of each scan; 3. Manual tracing of the contour of each section, including positions of selected primordia serving as landmarks, and the boundary of the phyllotactic pattern defining the morphogenetically active zone; 4. Interpolation of these data.
3 The phyllotactic patterning model (P-model) The P-model operates on the receptacle surface obtained by revolving the receptacle contour around the head axis (Fig. S9b,e,h).To generate the phyllotactic pattern in a gerbera head, we used the model described by Zhang et al. (2021), augmented with the information on the fate of floret primordia (trans vs. disc) (Fig. S19c).To generate the phyllotactic patterns in sunflower and Bellis heads (Fig. S9f,i) we modified the gerbera model by using C-models of receptacle contour expansion based on the sunflower and Bellis data, and sizes of primordia estimated from measurements of sunflower and Bellis heads.
4 The receptacle model (H and G-models)

Overview
Our data consist of 3D scans of heads at different developmental stages, as opposed to time-lapse images of the same developing head.With these data we were able to identify primordia that served as landmarks for determining local growth rates of the receptacle surface, but we were unable to identify internal landmarks that would provide information on the local growth within the head.Consequently, we approximated the growth of the receptacle interior on the basis of the expansion of its surface.We assume that the receptacle grows in a radially symmetric manner.It is thus sufficient to simulate the growth of one half of a single longitudinal cross-section only, since all sections grow the same way and can be obtained by rotating a single (half-)cross-section around the head axis (Fig. S9a,b,d,e,g,h).We conceptualize this cross-section as a (visco)elastic membrane attached to the section contour, although in our model it is defined in geometric rather than physical terms.In each simulation step the contour expands as described by the C-model, slightly deforming the membrane attached to it.This deformation changes the positions of points embedded in the membrane -and thus in the interior of the receptacle -including those representing the developing vascular system.
Mathematically, our model of the receptacle interior growth (G-model) is based on caging: the idea of smoothly deforming complex objects by embedding them into a geometric structure -the cage -manipulated with a small number of control points (Joshi et al., 2007).A historical predecessor to caging is the transformation method proposed by Thompson (1942) to describe how organisms change their shape during evolution or development.With that method, the shape of interest was positioned in a Cartesian coordinate system and deformed through linear or non-linear transformations of this system.Instead, to model the growth of the receptacle interior we employ harmonic coordinates: a flexible generalization of barycentric coordinates allowing for an arbitrary number of control points to be placed on the receptacle contour (Joshi et al., 2007).As the control points move due to contour expansion, points with given coordinates follow along, simulating the growth of the receptacle interior.

Calculation of harmonic coordinates (H-model)
Let Ω denote one half of the receptacle cross-section, and ∂Ω be the boundary of Ω.We partition this boundary into segment Γ representing the contour, and segment Γ m that closes ∂Ω along the head axis (Fig. SMD3): Assume that Γ is defined by a set of control points c i ∈ Γ, i = 1, ..., N .Some of these points match landmarks such as bracts, and others may be inserted inbetween to ensure sufficiently precise contour representation.Following Joshi et al. (2007), the position of an arbitrary point p ∈ Ω can be defined as the affine combination (Goldman, 2002) of the control point positions, in which weight fields (functions Each field w i , i = 1, . . ., N, is calculated assuming boundary conditions (BC): where δ i,j is the Kronecker delta function: Taking Equation 5into account, boundary condition (3) specifies that weight field w i evaluates to 1 at control point c i and to 0 at all remaining control points.It also implies that weights w 1 (p), . . ., w N (p) of any point p ∈ Ω sum to 1, as required by the definition of the affine combination of points.Boundary condition ( 4) ensures that all points on the head axis Γ m remain on the axis, as a move across or away from the axis would violate the assumption of head symmetry.The N -tuple of weights w 1 (p), . . ., w N (p) constitutes harmonic coordinates of point p.They are sample values of the harmonic coordinate fields w 1 , . . ., w N at point p.
We solve Laplace's equation numerically for each contour in the developmental sequence produced by the C-model.To this end, we first (Delaunay) triangulate the Ω domain using the Triangle library (Shewchuk, 1996).Laplace's equation is then discretized for this triangulation using the cotangent formula (Crane et al. 2013, see also Ringham et al. 2021) implemented in the libigl library (Jacobson et al., 2018).This discretization yields a system of linear equations associated with the triangle vertices (Jacobson, 2013), which we solve using the simplicial Cholesky factorization algorithm in the eigen library (Guennebaud et al., 2010).This computation is repeated for all N different boundary conditions (choices of control point c i at which the weight field is equal to 1), yielding N weight fields w 1 , . . ., w N .The results are stored in the coordinate fields file, on which the receptacle growth model (G-model) subsequently relies (Fig. SMD1).
Three examples of weight fields w i are shown in Fig. SMD4.Weight values for points p inside the triangles were calculated by linearly interpolating the weights at the triangle vertices using barycentric coordinates (Anisimov et al., 2016).Consistent with the membrane analogy, each control point c i affects weights w i (p) of points p positioned close to it more strongly than weights of points p positioned further away.

The receptacle growth model (G-model)
The input to the receptacle growth model is a sequence of M + 1 harmonic coordinate fields w t 1 , . . ., w t N representing snapshots of the receptacle development at time steps t = 0, . . ., M , where M is the number of simulation steps.We decrease the simulated time interval between consecutive steps t as the receptacle grows, because the geometry of larger receptacles changes faster than the geometry of receptacles at the onset of their development.With these harmonic coordinate fields in place, the propagation (change of position) of a point p ∈ Ω t resulting from the receptacle growth from stage t to stage t + 1 is modelled by expressing its position at time t, p t , using harmonic coordinate values (weights) w t i (p t ) associated with section Ω t , and using the same values for the subsequent section Ω t+1 (Joshi et al., 2007): As the control points move, point p will move in a "shape-aware" manner, i.e., its position will be most strongly affected by the contour points nearby.Harmonic coordinates are well suited for both convex and concave control cages, but large deformations from convex to concave configurations may reposition some interior points p ∈ Ω t outside receptacle domain Ω t+1 .With small time steps, the changes in the receptacle shape between stages t and t + 1 are also small, and we have rarely encountered this problem in practice.Nevertheless, we test for its occurrence, and if it happens, we project the misplaced point to its nearest location in contour ∂Ω t+1 .
5 The vascular patterning model (V-model)

Background
In constructing the vascular model we assumed that the basic processes of vascular patterning in heads are similar to those analyzed in inflorescences of plants not forming heads, such as Arabidopsis, tomato and Brachypodium (Reinhardt et al., 2003;Bayer et al., 2009;O'Connor et al., 2014).In brief, as the meristem grows, maxima of auxin concentration emerge in the active (peripheral) zone of the epidermis through feedback between auxin transport and the positions of PIN1 auxin efflux carriers.These maxima are sources of auxin flowing towards previously formed vascular strands, acting as sinks.Following the canalization hypothesis (Sachs 1969(Sachs , 1991)), the flow of auxin is patterned into narrow streams or canals, along which the vascular tissues differentiate.Experimental evidence for the above processes has been supported by cell-level computational models, which simulate the formation of auxin maxima in the epidermis (Smith et al., 2006;Jönsson et al., 2006) and the formation of canals of auxin flow between sources and sinks (Mitchison, 1980;Rolland-Lagan and Prusinkiewicz, 2005), and integrate both processes (Bayer et al., 2009;O'Connor et al., 2014); see (Cieslak et al., 2022) for a recent review).
Modeling large multicellular structures at the cellular level is a computationally daunting task.Consequently, models abstracting cell-level processes into higher-level geometric constructs have been proposed to simulate vascular patterning of entire plant organs (Cieslak et al., 2022).In the particle-based model introduced for leaves by Rodkaew et al. (2002Rodkaew et al. ( , 2003) ) and extended to developing inflorescences by Owens et al. (2016), vascular strands originate at the organ periphery, then extend towards the organ base by the addition of short segments ("particles") at the vein tips.With the space colonization algorithm proposed by Runions et al. (2005), leaf veins also extend by the addition of short segments, but patterning proceeds in the opposite direction: from the leaf base towards its margin.These dynamics are simplified further in the resistance minimization model proposed by Runions et al. (2017), in which entire vein segments connecting auxin sources to the previously formed vasculature are inserted in single simulation steps.Our model of vascular patterning in flower heads is based primarily on the resistance minimization algorithm, which is well suited to simulated experiments and refinements (Runions et al., 2017).

Vascular strand representation
We model vasculature as a directed graph V, E with vertices v ∈ V connected by edges e ∈ E. The edges represent segments of vascular strands, and the vertices represent points at which these segments meet.Some vertices may represent branching or joint points of several strands, while others may represent connecting points between short segments of a longer strand.The ground tissue is not represented explicitly and is defined by the absence of vasculature.The impact of receptacle growth on the shape of vascular strands is approximated by expressing positions of all vertices v using harmonic coordinates, which deform as the receptacle expands.Intuitively, we may thus think of the harmonic coordinate fields as playing the role of the ground tissue, in which the vasculature is embedded.Note, however, that the harmonic coordinates of all vertices are recomputed after each growth step, and thus the mathematical description of this embedding changes from one step to the next.

Vein initials
The patterning of each vascular strand begins with the placement of a primordium on the receptacle surface by the phyllotaxis model (P-model).This placement is followed by the formation of a short vein initial that penetrates the receptacle interior in a direction approximately normal to the receptacle surface.Biologically, it may originate in the supported bract or flower, and extend into the subepidermal tissue, where it connects to other strands of the emerging vascular network.In the following description we refer to this connecting point as the primordium, even though it is at some distance from the actual primordium returned by the P-model.

Minimum resistance path
Once a primorium p has been established, a strand segment that connects it to the already existing vasculature is formed.The connecting vertex v is chosen such that the resistance between the primordium and the procambial ring located at the receptacle base is minimized: Here R G (p, v ) is the resistance of the path from primordium p to a "candidate" vertex v passing through the ground tissue, and R V (v , ring) is the resistance of the path from v to the procambial ring passing through the already differentiated vascular strands.While modeling leaf vasculature, Runions et al. (2017) assumed that resistance R G (p, v ) was proportional to the (Euclidean) distance d(p, v ) between points p and v : where r G is the resistivity (resistance per unit distance) of the ground tissue.
Similarly, the resistance R V (v , ring) between point v and the ring through the existing vasculature V was calculated as where r V is the resistivity of vascular strands, and l e denotes the length of the strand segment e in the path E v connecting vertex v to the procambial ring.Equations 8 and 9 imply that resistivities of the ground tissue (r G ) and vascular strands (r V ) are constant scalar values.With this assumption it is possible to simulate vascular development in diverse leaves (Runions et al., 2017), but to reproduce vascular patterns in heads we needed to vary resistivities according to time and direction.

Time dependence
Vascular strands supporting the first few bracts in all three modeled species run approximately parallel to each other (Figs. 5 and S6b).To promote direct connections of these bracts to the procambial ring without forming a branching structure, we assumed that the ground tissue resistivity r G is low in very young heads.It increases to its target value used in the remainder of the simulation after the first few strands have been patterned.
An inherent feature of spiral phyllotaxis is that primordia are arranged into parastichies: conspicuous spirals running in opposite directions.In some flower heads, such as Bellis, developing vascular strands consistently follow a specific parastichy family (Fig. 5i).In contrast, strands simulated using the minimum resistance algorithm with constant resistivities switch their course between different parastichies in an irregular manner (Fig. S10a).To simulate the development of vascular strands that follow a specific family of parastichies, we extend the minimum resistance model by assuming that ground tissue resistivity r G and vascular strand resistivity r V depend locally on time in a manner promoting connections to older strands.Specifically, ground tissue resistivity r G (p, v ) adjusted for the age of nearby veins is calculated using the formula where S p is a set of vertices u close to vertex p. Resistivity r age G (p, v ) is thus increased with respect to the reference value r G for vertices v younger than the oldest vertex u in proximity of p. Correspondingly, Equation 8 is modified to Resistivity of an existing strand e is assumed to change from the ground tissue resistivity r G to its final value r V over maturation time τ m .This change is effected using a sigmoidal function r age V (e) of the segment age measured from the time of the segment creation.Equation 9 is thus modified to: Taking the time-dependent resistivities into account, the path of least resistance is found by substituting Equations 11 and 12 into Equation 7.An example of a vascular system following a single family of parastichies obtained in this manner is shown in Fig S10b.

Reticulation
To complete the reticulate vasculature exemplified by Bellis, the vascular connections that correspond to the opposite family of parastichies must also be produced.Data (Fig. 5c) indicate that these connections are formed after strands from the first family of parastichies, as opposed to appearing concurrently.We find them by searching for the minimum-resistance connection using Equation 7again, after a delay τ d , while excluding vertices v that already belong to the previously formed vascular strand originating at p. Once a candidate vertex v has been found, the total resistance R(p, ring) of the path passing through v is compared to the resistance of the original path.Unless the resistance of the new path is disproportionally high (the ratio of the resistances of the new and old paths is above a predefined threshold κ), the connection is made.The effect of this process is illustrated in Fig. S10c.

Sectoriality
In sunflower heads, vascular strands supporting floret primordia penetrate the receptacle deeply, forming highly branching tree-like structures attached to the radially spreading abaxial strands.These tree structures extend sectorially, which gives them an almost two-dimensional character.We account for this sectoriality by assuming the presence of a mechanism that polarizes the ground tissue radially, producing an anisotropic resistivity field.To model it, we extend the minimum resistance framework by introducing a non-linear function that increases ground tissue resistivity r G (p, v ) as point v deviates from the longitudinal section p occupies.Some degree of sectoriality is also observed in other heads; the impact of its magnitude is illustrated in Fig. S11.
We use the perpendicular distance between v and the radial longitudinal sector including p to modulate ground tissue resistivity r G (p, v ) where n is the normal of the sector including p. Parameter µ sec thus controls the increase of ground tissue resistivity in directions not aligned with the longitudinal sector; µ sec = 0 disables this feature.Sectoriality can be combined with time dependence by replacing r G in Eq. 10 with the function r sec G (p, v ) (Eq. 13).

Adaxial vascular strands of gerbera
A distinctive feature of gerbera heads is the presence of adaxial strands that originate near the head rim and extend radially towards the head center (Figs. 2b,c and S1b,c).These strands have a different character from both the Bellis strands, which follow parastichies, and the sunflower strands, which form branching structures deeply penetrating the receptacle.Remarkably, the adaxial strands of gerbera develop ahead of the florets they support, which indicates that their development is not guided by the floret primordia (Figs.5b,c and  S5d).To model these strands, we assume the presence of a diffuse auxin source, represented geometrically as a set of "attraction points", towards which the adaxial strands gradually extend.The attraction points are distributed subepidermally ahead of the active ring, such that they maintain constant density.The strand extension process is simulated using the space colonization algorithm (Runions et al., 2005(Runions et al., , 2007)).The veins connecting florets to the already formed vasculature are then patterned using the resistance minimization algorithm, taking the already formed adaxial strands into account.

Randomization
Single-shot vascular connections generated by the resistance minimization algorithm are straight lines connecting primordia (more precisely, the ends of vein initials originating at the primordia) to the existing vasculature.In contrast, the strands observed in heads meander, deviating from straight lines to various degrees.One source of these deviations is the irregular lattice of ground tissue cells within which the strand patterning and differentiation occur.As our geometric model is an abstraction that does not explicitly simulate ground tissue cells or other possible factors in path deviation (except for the random distribution of attraction points inherent in the space colonization algorithm), we introduce deviations by randomly perturbing the vascular paths.To this end, we build a strand from p to v as a polygonal line with segments of approximately equal length δ step joined at vertices p 0 = p, p 1 , p 2 , . . ., p n = ṽ.The endpoint ṽ is an existing vascular vertex in the proximity of, but not necessarily coinciding with, v.The process is controlled by two parameters: the probability η that a perturbation will take place, and the maximum magnitude of this perturbation, γ max .The polygonal trajectory is constructed according to the formula: where ŝt is the unit vector from vertex p t towards the target vertex v, dt is a randomly chosen direction of deviation perpendicular to ŝt (in 3D), and γ t is the magnitude of this deviation, chosen randomly (with uniform distribution) from the interval [0, γ max ].The strand construction ends when the most recently constructed vertex p t is within a predefined threshold distance from a vertex ṽ of the existing vasculature, at which time the connection is completed by adjusting position of p t to coincide with ṽ.

Calculation of strand diameter
The model elements described so far are focused on generating the skeletons of vascular structures.To facilitate visual comparisons between scanned and simulated patterns, we also model the widths (diameters) of the vascular strands.We adapt for this purpose the method proposed for leaf venation by Runions et al. (2005).The computation proceeds basipetally, from the organ primordia to the receptacle base.Veins directly supporting organ primordia, as well as the tips of adaxial strands in gerbera, are assigned predefined widths ρ specified as model parameters.Vascular segments that extend a strand outside of branching points have the same diameter as the preceding strands.At the points where strands with diameters d 1 and d 2 meet in their course toward the receptacle base, diameter d 12 of the resulting, more basal strand is calculated using Murray's (1927) law, where n is a model parameter (MacDonald, 1983).These rules suffice to determine strand widths in open vascular patterns (without loops).
In the reticulated vasculature of Bellis (Fig. 7) and its less regular variants (Fig. S11), we also consider configurations in which vascular strands split.While modeling leaf venation, Runions et al. (2005) assumed that veinlets reaching a reticulation point contribute evenly toward the width of the outgoing strands.In the model of Bellis, we modify this heuristic by assuming that the strand or strands reaching a reticulation point contribute to the width of outgoing strands unevenly, favoring the outgoing strand that was formed first.This asymmetry results in different diameters of the strands following opposite parastichies (Fig. 7), as observed in nature (Fig. 3g).

Figure SMD1 :
Figure SMD1: Data flow in the dynamic model of vascular patterning in heads.Dashed lines indicate model components implemented as separate programs.Bold font indicates data exchanged as files.Red arrows highlight the main iteration loop.

Figure SMD2 :
Figure SMD2: Simulation process in a section of a head.(a) Initialization.A receptacle (green) and procambial ring (blue dots) are established.(b) Receptacle growth.Initially no primordia or vascular strands are present.(c) Patterning.New primordia (red) are positioned on the receptacle surface according to the P-model, and strands (brown) connecting these primordia to the procambial ring are established by the V-model.(d) Next growth step.Vascular strands embedded in the receptacle are deformed, and primordium positions on the receptacle surface are adjusted by the G-model as the receptacle grows.(e) Result of iterating steps (c) and (d) several times.As the simulation progresses, new primordia emerge, giving rise to strands that connect these primordia to the procambial ring or to the strands formed previously.Arrows indicate additional adaxial strands specific to gerbera.

Figure SMD3 :
Figure SMD3: Elements of a receptacle cross-section representation.

Figure SMD4 :
Figure SMD4: Three examples of harmonic weight fields associated with receptacle sections.The influence of the selected point c i (red) on field w i is the highest (red) close to c i and decreases (white to blue) further away.

Table SMD1 :
Numerical values of key model parameters used in the simulations of the Bellis, sunflower and gerbera heads.Blue rows show quantities derived from other parameters or emerging in the simulations.