Universal Hyperuniform Organization in Looped Leaf Vein Networks

Recently, it was discovered that apparently disordered organizations in biological systems across scales, such distributions of photoreceptor cones in avian retina Jiao et al. (2014), self-organization of immune cells Mayer et al. (2015), and large-scale patterns of vegetation Ge (2023), can possess an exotic hidden order called “hyperuniformity”. Disorder hyperuniform (DHU) systems lack conventional long-range order, yet they completely suppress normalized infinite-wavelength density fluctuations like crystals Torquato and Stillinger (2003); Torquato (2018). In this sense, DHU systems can be considered to possess a “hidden order” in between that of a perfect crystal and a completely random system (e.g., an ideal gas). DHU is characterized by a local number variance $\sigma_{N}^{2}(R)$ associated with a spherical window of radius $R$ that grows more slowly than the window volume (e.g., $\sim R^{d}$ in $d$-dimensional Euclidean space) in the large-$R$ limit Torquato and Stillinger (2003); Torquato (2018), i.e., $\lim_{r\rightarrow\infty}\sigma_{N}^{2}(R)/R^{d}=0$. This is equivalently manifested as the vanishing static structure factor in the infinite-wavelength (or zero-wavenumber) limit, i.e., $\lim_{|{\bf k}|\rightarrow 0}S({\bf k})=0$, where ${\bf k}$ is the wavenumber and $S({\bf k})$ is related to the pair-correlation function $g_{2}({\bf r})$ via $S({\bf k})=1+\rho\int e^{-i{\bf k}\cdot{\bf r}}[g_{2}({\bf r})-1]d{\bf r}$ and $\rho=N/V$ is the number density of the system. For statistically isotropic systems, the structure factor only depends on the wavenumber $k=|{\bf k}|$. The small-$k$ scaling behavior of $S(k)$, i.e., $S({k})\sim k^{\alpha}$, determines the large-$R$ asymptotic behavior of $\sigma_{N}^{2}(R)$, based on which all DHU systems can be categorized into three classes: $\sigma_{N}^{2}(R)\sim R^{d-1}$ for $\alpha>1$ (class I); $\sigma_{N}^{2}(R)\sim R^{d-1}\ln(R)$ for $\alpha=1$ (class II); and $\sigma_{N}^{2}(R)\sim R^{d-\alpha}$ for $0<\alpha<1$ (class III), where $\alpha$ is the hyperuniformity exponent Torquato (2018).

Besides the aforementioned DHU examples in biological systems, a wide spectrum of equilibrium Torquato et al. (2015); Batten et al. (2009); Jiao (2022) and non-equilibrium Gabrielli et al. (2002); Donev et al. (2005); Zachary et al. (2011a, b); Yuan et al. (2021) many-body and material systems, in both classical Kurita and Weeks (2011); Hunter and Weeks (2012); Dreyfus et al. (2015); Chen et al. (2023a); Zhang et al. (2023); Chremos and Douglas (2018) and quantum mechanical Torquato et al. (2008); Feynman and Cohen (1956); Gerasimenko et al. (2019); Sakai et al. (2022); Rumi et al. (2019); Sánchez et al. (2019); Zheng et al. (2020); Chen et al. (2021, 2022, 2023b) varieties, have been identified to possess the property of hyperuniformity. Notably, certain driven non-equilibrium systems Hexner and Levine (2015); Jack et al. (2015); Weijs et al. (2015); Salvalaglio et al. (2020); Nizam et al. (2021); Zheng et al. (2023), active-particle fluids Lei et al. (2019); Lei and Ni (2019); Huang et al. (2021); Zhang and Snezhko (2022); Oppenheimer et al. (2022), and dynamic random organizing systems Hexner and Levine (2017); Hexner et al. (2017); Weijs and Bartolo (2017); Wilken et al. (2022), which have been used as models for, e.g., self-organization of biological cells, school of fish or flock of birds, are also found to exhibit emergent hyperuniform behaviors. Novel DHU materials have been engineered that can possess superior properties compared to their crystalline counterpart, such as high-degree of isotropy and robustness against defects for transport applications Florescu et al. (2009); Klatt et al. (2022); Zhang et al. (2016); Chen and Torquato (2018); Torquato (2021); Xu et al. (2017); Leseur et al. (2016); Yu (2023); Shi et al. (2023); Torquato and Chen (2018).

Although it is not completely understood how DHU emerges in many of the aforementioned systems, the attempt to achieve certain functional optimality under constraints seems to be a crucial underlying mechanism Ge (2023). In the case of avian photoreceptors, the cones of the same type favor the ordered triangular-lattice arrangement that maximizes light sampling efficiency. However, this ideal arrangement is frustrated by the size polydispersity, shape irregularity and incompatible number ratios of five cone types. Therefore, the resulting DHU arrangement is one that optimizes light sampling under these constraints Jiao et al. (2014). Similarly, it was found the DHU pattern of arid and semi-arid vegetation ecosystem is likely resulted from optimal utilization of water resources under constraints Ge (2023). These exciting discoveries naturally lead to the question: Are there any other biological systems that have developed hyperuniform organizations for optimal functionality under constraints?

In this letter, we report the discovery of universal hyperuniform organization in looped vein networks. Leaf veins are present throughout the leaf lamina, providing mechanical support of the leaf structure Scarpella et al. (2010); Biedroń and Banasiak (2018). Another key functionality of leaf veins is to transport water, minerals and nutrients to the leaf cells. These veins form a hierarchical network: On the macroscopic scale, the thick primary veins organize into a branched network, which is very efficient to mechanically support the leaf structure, see Fig. 1(a). On the other hand, in between the main branches, the secondary thinner veins can form looped networks Katifori and Magnasco (2012) for efficient transport, see Fig. 1(b).

Quantitative characterizations of the venation patterns of angiosperm leafs have revealed strong local regularity ref1 ; ref2 and topological hierarchy Katifori and Magnasco (2012) of the vein networks, which are intimately related to their functionality. For example, it is found local vein density is strongly correlated with the transpiration rate ref1 , and the number of stomata within each areole was linearly related to the length of the looping vein contour ref2 . Katifori and Magnasco Katifori and Magnasco (2012) introduced a novel hierarchical loop decomposition framework to quantify the hierarchical topological organization of the looped venation patterns. Ronellenfitsch and Katifori ref3 showed that the development of vein networks spanning from looped to tree-like patterns corresponds to an optimal trade-off between transport efficiency, cost and robustness against defects, a property that has been observed in many disordered hyperuniform systems Torquato (2018).

Our detailed analysis below indicates the distribution of the geometrical centers associated with the cellular structures in the looped vein network, for all leaf species analyzed here, possesses a high degree of hyperuniformity, characterized by a vanishing structure factor in the zero-wavenumber limit, i.e., $S(k)\sim k^{\alpha}$ with $\alpha\approx 0.64\pm 0.021$, providing a novel example of class III hyperuniformity in biology. Importantly, we show the DHU organization of the vein networks leads to superior transport behaviors, compared to other uncorrelated or correlated disordered organization models, which is crucial to the functionality of the veins.

Sample preparation and image processing. We collected and analyzed a large number of air-dried leaf samples of 6 distinct species including ficus religiosa (FR), ficus caulocarpa (FC), ficus microcarpa (FM), smilax indica (SI), populus rotundifolia (PR), and yulania denudate (YD), commonly found in middle and northern part of Asia, Europe and North America. The selected species are also representative of a larger group of plant species with looped leaf venation ref1 . We took high-resolution optical images of the samples (see Fig. 1(a) for an example of ficus religiosa leaf), and used Matlab for image correction and enhancement to highlight the vein network (Fig. 1(b)). We employed the OpenCV package in Python to identify the “cells” formed by the veins based on two criteria: (i) a cellular region is enclosed by a closed loop of veins; and (ii) a cellular region contains an open-end branched vein in the center (see SI for details). Such open-end branched structures play a similar role of capillaries, releasing the minerals and nutrients to individual leaf cells Scarpella et al. (2010). An example of the resulting cellular structure is shown in Fig. 1(c), based on which the geometrical center of each cell is obtained (with number $N$ varying from $\sim 100$ to $\sim 400$) for subsequent analysis, see Fig. 1(d). The average number density of derived point configurations is $\rho\approx 125$ cm^-2.

Hyperuniform organization of cellular structures of vein networks. We obtain statistics of the point configurations derived from the cellular structures of the secondary vein networks, including the static structure factor $S(k)$, the number variance $\sigma_{N}^{2}(R)$, as well as the pair-correlation function $g_{2}(r)$, by averaging over 20 independent realizations randomly selected from the leaf samples (see SI for details). To better understand the unique hyperuniform organizations of these structures, we also consider two non-hyperuniform benchmark systems at the same number density: (i) a Poisson distribution of points and (ii) a hard-disk packing generated via random sequential addition (RSA). The former represents a totally uncorrelated disordered system, while the latter represents a correlated system through the mutual exclusion effects.

Figure 2(a) and (b) respectively shows $S(k)$ in the linear and log scale for the three systems. As expected, the uncorrelated Poisson system possesses a constant $S(k)=1$, while the RSA system shows a strong suppression of scattering for small $k$, but is non-hyperuniform. On the other hand, the point configurations derived from the cellular structures of the leaf vein network exhibit a high degree of hyperuniformity, characterized by a vanishing $S(k)$ in the zero-$k$ limit with a hyperuniformity index Torquato (2018) $H=S(0)/S(k_{p})\approx 10^{-3}$ (where $k_{p}$ is the wavevector associated the first and highest peak and $S(0)=\lim_{k\rightarrow 0}S(k)$) and a hyperuniformity exponent $\alpha\approx 0.64\pm 0.021$ (illustrated as the dashed line), corresponding to class III hyperuniformity. To understand the boundary effects, we also carry out a system-size study and confirm that the same small-$k$ behavior and hyperuniformity exponent can be robustly obtained in systems with $N\sim 150$ to $N\sim 400$ (see Fig. S2 in SI), with the latter corresponds to the largest patches of secondary veins in between the two main branches.

Figure 2(c) shows $\sigma_{N}^{2}(R)$ of the three systems in the log scale, where the Poisson and RSA systems follow the scaling $\sigma_{N}^{2}(R)\sim R^{2}$ and the leaf vein system possesses $\sigma_{N}^{2}(R)\sim R^{\beta}$, with $\beta\approx 1.35\pm 0.017$, consistent with the estimate from $S(k)$ analysis, i.e., $\beta=2-\alpha\approx 1.36$. Figure 2(d) shows $g_{2}(r)$, characterizing the short-range correlations. Again, the Poisson system possesses a $g_{2}(r)=1$, indicating its totally random nature. $g_{2}(r)$ of the RSA system possesses a clear exclusion region, corresponding to the diameter of the congruent circular disks, which is immediately followed by a sharp peak, resulted from the nearly contacting particle configurations. $g_{2}(r)$ of the vein networks is reminiscent of that of poly-disperse disk packing, characterized by a well-defined exclusion region, followed by a smooth increase and a broad peak. Similar $g_{2}(r)$’s were also observed for the hyperuniform distribution of the photoreceptors in avian retina Jiao et al. (2014).

Superior transport properties induced by hyperuniformity. As noted above, one main function of the secondary vein network is to transport water and minerals to the leaf cells. As shown in Fig. 1b, there is an open-end branched structure in the middle of each “cellular region” of the secondary vein network, from which the minerals can be released and diffused to the individual cells. Here we employ the time-dependent diffusion spreadability $\mathcal{S}(t)$ Torquato (2021) to quantify the efficiency of diffusive transport supported by the DHU structures. Specifically, a fixed amount of solutes are initially confined within a subset of regions $\mathcal{V}_{2}$ (e.g., within the veins), which then start diffusing into the remaining regions $\mathcal{V}_{1}$ at $t=0$, assuming the diffusivity $D$ is identical for both phases (see inset of Fig. 3b). $\mathcal{S}(t)$ is defined as the fraction of total solutes present in $\mathcal{V}_{1}$ as a function of time $t$, which is a measure of the “spreadability” of the solutes as a function of time. A related quantity is the excess spreadability $\mathcal{S}(\infty)-\mathcal{S}(t)$, i.e.,

where $\phi_{2}$ is the volume fraction of $\mathcal{V}_{2}$, and $\tilde{\chi}_{V}({\bf k})$ is the spectral density function associated with $\mathcal{V}_{2}$, which is the Fourier transform of the corresponding autocovariance function $\chi_{V}({\bf r})=S^{(2)}_{2}({\bf r})-\phi_{2}^{2}$ and $S^{(2)}_{2}({\bf r})$ is the two-point correlation function of $\mathcal{V}_{2}$ Shi et al. (2023). The large-$t$ behavior of $\mathcal{S}(\infty)-\mathcal{S}(t)$ quantifies how fast the solutes can achieves the asymptotic uniform distribution in the system, reflecting the efficiency of diffusive transport.

Figure 3a shows the excess spreadability $\mathcal{S}(\infty)-\mathcal{S}(t)$ of the DHU vein networks by directly computing $\tilde{\chi}_{V}({\bf k})$ of the “vein phase” from the imaging data. Also shown is $\mathcal{S}(\infty)-\mathcal{S}(t)$ of a Voronoi network derived from the point configuration associated with the vein networks, including congruent disks centered at those points. The size of the disks is chosen such that their area is equivalent to the average area of the open-end branched structures in the areola. This simple model captures very well both short-time and long-time behavior of $\mathcal{S}(\infty)-\mathcal{S}(t)$ for the actual vein networks. Notably, the DHU vein network possesses a long-$t$ asymptotic scaling $\mathcal{S}(\infty)-\mathcal{S}(t)\sim 1/t^{\gamma}$, where $\gamma\approx 1.38$. This is to contrast the scaling $\mathcal{S}(\infty)-\mathcal{S}(t)\sim 1/t$ for a general nonhyperuniform disordered system Torquato (2021).

To further illustrate the diffusive transport efficiency of the DHU structures, we compare in Fig. 3(b) $\mathcal{S}(\infty)-\mathcal{S}(t)$ of three model structures derived from the Poisson, RSA and vein-network point configurations, see Fig. 3(c). For the DHU system, $\tilde{\chi}_{V}({k})=\rho\tilde{m}^{2}(k)S(k)$ for the disk regions, where $\tilde{m}(k)$ is the Fourier transform of the shape function of the disks, and $\tilde{m}(0)$ is a constant. Thus, the small-$k$ scaling $\tilde{\chi}_{V}({k})\sim S(k)\sim k^{\alpha}$ leads to $\mathcal{S}(\infty)-\mathcal{S}(t)\sim 1/t^{(1+\alpha/2)}$. The scaling analysis from the numerical data yields $\sim 1/t^{1.33}$, which is consistent with the calculated exponent $\alpha\approx 0.64$. This indicates uniform solute distribution can be achieved much faster in the DHU system than in other two disordered systems, characterized by $\mathcal{S}(\infty)-\mathcal{S}(t)\sim 1/t$. These analyses clearly indicate the superior diffusive transport efficiency of the DHU networks, compared to other disordered organizations.

Universal hyperuniformity in looped vein networks across different leaf species. Last but not least, we present the analyses on the leaf samples of other species including ficus caulocarpa (FC), ficus microcarpa (FM), smilax indica (SI), populus rotundifolia (PR), and yulania denudate (YD), as shown in Fig. 4(a). The statistics including $S(k)$, $\sigma_{N}^{2}(R)$ and $g_{2}(r)$ for all six leaf species, when rescaled to the same number density, all collapse on to the corresponding universal curves, as shown in Fig. 4(b)-(e). Although the small-$k$ behaviors of $S(k)$ exhibit some small discrepancy, possibly due to the finite-size effect, the estimated hyperuniformity exponent $\alpha$ for all systems agree reasonably well with one another, with the numerical values falling into the narrow interval [0.626, 0.648]. These results indicate the secondary vein networks for all six leaf species exhibit a universal hyperuniform organization on large scales, characterized by $S(k)\sim k^{\alpha}$ with $\alpha\approx 0.64$ which belongs to Class III hyperuniformity. In addition, the different leaf species also possess very similar local structures when properly rescaled, manifested as the collapsed pair-correlation functions.

How such universal DHU organizations emerge during the development and regulation of the vein network? It is well established that the leaf vascular system develops in response to the flow and localization of the growth regulator auxin as continuous strands of conducting tissues arranged in repetitive spatial patterns Scarpella et al. (2010); Biedroń and Banasiak (2018); Fujita and Mochizuki (2006a, b). Although the auxin diffusion implies a regular length scale above which a new vein would be created, we show in SI that such local regularity and topological hierarchy do not spontaneously lead to the observed hyperuniformity, using both an ideal geometrical model that aligns points with a narrow distance distribution along the main vein network and a more sophisticated growth model ref4 explicitly incorporating the diffusion of auxin and considers the hierarchy of vein growth. It has been shown ref3 that hierarchically organized reticulation in venation patterns can be constructed and maintained through spatially correlated load fluctuations across length scales, offering a potential mechanism for regulating the vein networks to induce hyperuniformity.

In summary, we have discovered a universal hyperuniform organization in looped leaf vein networks across a wide spectrum of representative tree species, which does not seem to be revealed by previously reported characterizations of the venation Katifori and Magnasco (2012); ref1 ; ref2 , nor captured by previous models of leaf vein development ref3 ; Fujita and Mochizuki (2006a, b); ref4 . Such DHU organization leads to superior efficiency of diffusive transport in the system. Although globally optimal transport is only achieved by periodic structures, the geometrical constraints and frustrations (e.g., curvature of primary branches and relatively short inter-branch distances) rule out the possibility to develop perfectly ordered venation patterns. In this sense, the DHU venation pattern may represent an sub-optimal organization of secondary vein networks for transport under constraints induced by the primary vein branches.