Statistical hydraulic model for the Leonardo’s rule

Leonardo da Vinci was a man with innate curiosity and exceptional intellectual capacity that allowed him to understand nature [1]. Leonardo is best known for his iconic paintings the Mona Lisa and the Last Supper [2] as well as his remarkable drawing the Vitruvian Man [3]. However, his legacy goes beyond arts with seminal contributions to science. Based on observation and experience he deciphered different patterns in nature. In plants he unveiled the so-called sixth leaf, rings of growth and branching patterns in trees nowadays known as Leonardo’s rule [1, 3]. Regarding the latter, he found that the branches at different heights are proportional to the stem. In specific, the sum of the transverse areas of the branches at different bifurcation levels remains constant and equal to the transverse area of the stem. Leonardo, based on his knowledge of the transportation of water in plants, attributed this relationship to the conservation of sap flux [1]. The modern version of the Leonardo’s rule is given in Pythagorean fashion $d^{\Delta}=\sum^{n}_{i=1}d_{i}^{\Delta}$, where $d$ is the diameter of the stem, $d_{i}$ is the diameter of the $i$-th branch at a particular bifurcation level, $n$ the number of branches at the mentioned level and $\Delta$ the so-called Leonardo’s exponent [4]. The nominal value of the Leonardo’s exponent is 2, however it is well known that it depends on the tree species [4, 5, 6]. The Leonardo’s rules is also obeyed by the radical system of trees [7]. Furthermore, the Leonardo’s rule in conjunction with the external fractal architecture of trees is involved in the optimization of different biological processes [8, 9]. For instance, the modeling of trees shows that the mechanical stress induced by wind loads is minimized by the self-similar organization of the branches as well as the area-preserving condition [8]. Likewise, the competition for light is optimized due to the complex-fractal growth of trees in which the Leonardo’s rule is implicated [9].

There are also several works assessing the impact of hydraulics in plant vascular systems [10, 5, 11, 12, 13]. From the rather simple pipe model to the most sophisticated proposals based on internal and external fractal architectures. The common factor in all these studies is the derivation of well-known allometric scaling laws. In the case of the pipe model [10, 5, 11], trees are modeled as a uniform tube from the base to the top, with the tube diameter equivalent to the cross-section of the stem. The pipe model entails the Leonardo’s rule as well as the non-linear scaling relationships for the height and mass of a plant with respect to the basal stem diameter. The pipe model has also been useful to determine the biomass of plants. A historic review of the pipe model can be consulted in [11]. The WBE model is a more refined proposal in which a general theory of resource distribution through hierarchical branching networks is considered [12]. In particular, the model takes into account external space filling (self-similar branching), hydraulic optimization (tube tapering) and scale free canopy (invariant leaf-petiole size). In this model despite tubes taper the flow rate is size independent due to the increase of the flow velocity in small radius branches. The WBE model also shows that the well-known allometric scaling laws for animals [14] apply for plants [15], supporting the idea that they are universal in biology and can be modeled by simple geometric and hydrodynamic principles [16]. A decade later of the WBE model a theory that also considers internal space filling as well as trade-offs between hydraulic safety and efficiency was proposed [13]. The interrelationship between the vascular and branching networks give rise to more general allometric scaling laws that can respond to the particularities of tree species. In addition, the model gives predictions for the sap (solutes and water) flow, the taper and frequency of xylem conduits. In contrast to the WBE model the flow velocity remains constant regardless of the tapering of the xylem conduits, and the flow rate size independent as a result of the internal fractal architecture. The area preserving condition is also fulfilled in this model. As we have documented, it is well accepted that the external branching network minimizes the mechanical stress in trees, optimizes the capture of light for the photosynthesis process and shapes the fundamental allometric scaling laws. However, for the internal vascular network there is no consensus and depending on the model we can have different outcomes, as in the case of the fluid velocity. Finally, it is important to remark that in all these models the area preserving condition (Leonardo’s rule) and the size-independent flow rate (flow conservation) are obeyed. So, it is interesting to see how these conditions as fundamental determine by and large the size distribution of xylem conduits.

In this work we explore the seminal idea of Leonardo that the conservation of the sap flux is directly implicated in the Leonardo’s rule. We propose a statistical hydraulic model based on a power law distribution for the xylem elements. We corroborate the Leonardo’s rule by measuring the radius of the stem and branches of different tree species. We also obtain the size distribution of the xylem elements after subjected one of the species to a maceration process. A fairly good agreement is obtained between the theoretical prediction and the experimental results for the size distribution of the xylem elements.

As we have documented the hydraulic model has been modified and refined throughout years. From the rather simple pipe model to the more elaborated proposals of the last decades. In the present work, we adopt an statistical hydraulic version based on a power law distribution for the xylem elements that accounts of the space filling internal fractal architecture of trees as well as memory-correlation effects of the conduction process. We also consider that the flux through each xylem element is non ideal due to the intricacies related to the internal architecture of the conduction vessels [17, 18]. In addition, we assume the Leonardo’s seminal hypothesis as valid, that is, the Leonardo’s rule is a consequence of the sap flux conservation [1]. In fact, by attending these requirements we will show that the power law exponent is three and that the sap flux is optimized (ideal fluid behavior) as a consequence of the collective behavior of the xylem elements.

Our starting point is the mentioned Leonardo’s hypothesis [1]. In specific, if the stem is bifurcated in two daughter branches we can established a relationship between their corresponding fluxes as follows

The first task that we carried out it was the corroboration of the Leonardo’s rule by directly measuring the stem-branches radii of four different trees species. We measured the stem-branches of the first bifurcation of Quercus rugosa, Eucalyptus camaldulensis, Schinus molle and Prunus serotina. We have followed the measuring protocol of Aratsu [5]. In particular, we carefully choose healthy trees, with no damage in the branches and with an evident bifurcation. In addition, we have measured the stem and the daughter branches at a reasonable distance from the bifurcation to have reliable data, particularly, to avoid the possible data dispersion associated to the natural branch thickening around the bifurcation. We have considered 91, 50, 48 and 58 individuals for Quercus rugosa, Eucalyptus camaldulensis, Schinus molle and Prunus serotina, respectively.

In order to corroborate the Leonardo’s rule we will assume that the perimeter of the stem and daughter branches corresponds to the one of a circle. So, the Pythagorean relationship between radii Eq. (7) is also valid for the perimeters

The second task that we performed was to unravel the trees xylem elements through the so-called maceration process [19]. This process is recommended for young trees, thus, we implemented it only in Prunus serotina. In specific, we macerated three individual according to the following protocol:

1. 1.

   For each individual we prepare a 300 mL maceration solution, of which 120 mL are distilled water, 30 mL are hydrogen peroxide and the remaining 150 mL correspond to glacial acetic acid.

2. 2.

   With a handsaw we take pieces of the tree right above and below the bifurcations by cutting the branches transversely. In total we collect five pieces, three at the first bifurcation stage and two more at the second stage. From each piece we extract fine slices with a scalpel.

3. 3.

   The fine slices of each piece are placed in a sterile flask with 60 mL of maceration solution. The five flasks of each individual are sealed and introduced in an incubation oven for a five day period at a temperature of 60^∘C.

4. 4.

   After the time in the incubation oven each flask sample is subjected to three rinses in order to drain the glacial acetic acid. The remaining material is let to rest a whole night in distilled water.

5. 5.

   Finally, the distilled water is drained and the remaining solution centrifuge to obtain (separate) the xylem elements.

After the maceration process we proceed to observe the xylem elements. In order to do so, we put the resulting macerated solution in a 0.100 mm depth Neubauer chamber. Specifically, we put a drop of the solution in each of the squares of the chamber, protecting them with a cover glass. Once the sample is ready it is placed in an optical microscope under 100X magnification. We then proceed to scan from left to right and from up to down the observing area by taking pictures of it in the process, see Fig. 1. The pictures are processed with the software ImageJ. We consider two methods to obtain the size distribution of the xylem elements. One in which we measure directly the diameter of the xylem elements and other in which the images are black and white contrasted to measure the shaded areas of the xylem elements. More details of these methods can be found in the supplementary material.

In this section we present the main results of the field and experimental work. Most importantly, we contrast the resulting size distribution of the xylem elements after the maceration process with the statistical hydraulic model outcomes.

Regarding the field work, in table 1 we show the Leonardo’s exponent of the tree species studied. The exponents are obtained by computing the zeros of Eq. (8) and averaging the exponents of all individuals of each tree specie. The standard error of the exponents is also obtained and shown in table 1. As we can see the Leonardo’s exponent for all tree species lies in the well-known range $1.8<\Delta<2.3$ [8]. It is worth noting that the exponent for adult trees is above 2, while the corresponding one to young trees (Prunus serotina) is below 2. As the Leonardo’s exponent depends on the tree specie, we also expect that the size distribution of xylem elements does. In fact, in table 1 we also show the expected exponents for the xylem elements size distribution based on our statistical hydraulic model outcomes. As the only tree specie that was subjected to the maceration process was Prunus serotina we expect a size distribution exponent $x\approx 3.11$ in the experimental work.

Regarding the experimental work, we firstly show the results of measuring directly the diameters of the xylem elements after the maceration process. The size distributions for the three Prunus serotina individuals are shown in Fig. 2. As we can see the size of the xylem elements ranges from 8 $\mu$m to 30 $\mu$m. Furthermore, small size xylem elements are the vast majority, resulting in a power law distribution for all individuals. In fact, the best fitting function for the experimental data is of the form $f(r)=ar^{-b}$, see the solid-blue curves in Fig. 2. In all cases the correlation coefficient is over 90 %, indicating a fairly good fitting for the power law function. We also determine the relative error of the resulting exponent with respect to the expected one. The resulting exponents for the three individuals as well as its corresponding relative errors are shown in table 2. As we can notice the exponent of the third individual ($b\approx 3.13$) is closest to the expected one ($x\approx 3.11$), with a relative error of 0.64 %. Here, it important to mention that measuring directly the diameter of the xylem elements is time demanding and consequently impractical to process all images of the three Prunus serotina individuals. So, the exponents that we obtained are the result of processing a reduced number of xylem elements images.

In order to overcome this obstacle, we have used an alternative method based on black and white contrast the xylem elements images. This method allows us to measure in an automatic way the shaded areas of the xylem elements and to determine indirectly the corresponding diameters. More details can be found in the supplementary information. In Fig. 3 we show the outcomes of this alternative method for the first Prunus serotina individual. In this case we have measured the area of about $10,000$ xylem elements. As we can notice the size of xylem elements ranges from 20 $\mu$m to 180 $\mu$m, being the vast majority small size xylem elements. The experimental data fit quite well to a power law function (solid-blue curve), with a size distribution exponent $b\approx 3.08$ and a 99 % correlation coefficient. In addition, the relative error of the size distribution exponent is 0.96 %.

The first aspect that we want to discuss is related to the characteristics of the xylem elements. It is well known that as trees grow some xylem elements lose its conductive ability, as they fill with resin, giving place to the so-called heartwood [17]. At the same time new xylem elements arise, keeping the tree conducting capability. The elements that compose the xylem have in its walls lignin [17, 21], a substance that makes them rigid and woody, providing mechanical stability and support to the tree. Both conducting and non-conducting elements contribute to the thickness of the stem and branches. With the maceration process what we obtain is a disperse sample with conducting and non-conducting elements. Our first image processing method allows us to choose xylem structures with conducting capabilities such as tracheids and vessel elements. With the black and white contrast method we measure the areas of conducting and non-conducting elements. Though, heartwood structures are not participating in the conduction process, they are relevant in the determination of the size distribution because they contribute to the thickness of the stem and branches. Under this context, we expect results in better agreement with respect to the predicted ones with the black and white contrast method. In fact, this is what happened with the first individual. After considering conducting and non-conducting elements the exponent of the size distribution ($b\approx 3.08$) results in close agreement with respect to the predicted one ($x\approx 3.11$) with a relative error of 0.96 %.

The second aspect that we want to address it is about the type of distribution function that was obtained. In specific, a power law distribution function alike to the so-called Lévy-type distribution functions [22, 23, 24]. In fact, these distribution functions are characterized by memory and correlation effects. Taking into account that xylem structures form a complex conducting network [17], in which memory and correlation effect are plausible, we consider that it is totally reasonable to obtain a power law distribution for the xylem elements. In fact, a power law distribution, in which small size xylem elements dominate, optimized the conduction process by changing the non-ideal fluid behavior of individual xylem elements to an ideal one after summing over all xylem elements. Actually, at first sight, it is counterintuitive that a network of non-idea fluid xylem elements can optimized the conduction process. However, what really matters it is the collective effect of all conducting structures rather than its individual behavior. This intricate organization of the xylem elements it is likely to be related to the space-filling maximization as in the case of other biological complex networks [25, 26, 27]. Our results are also indicative that the generation of the xylem elements, from the meristematic tissue called Cambium [28], is taking place in such a way that small radius xylem structures dominate over large ones, giving rise to a power law distribution.

Finally, we want to discuss the results obtained with the black and white image-processing method. In principle, it is possible that with the maceration process a lot of conducting structures be broken. If so, what we are obtaining with the black and white image-processing method would be the distribution function of a rupture process [29, 30], as we are measuring automatically all structures in our sample. However, the exponent of the distribution function associated to a rupture process it is well known [29, 30] and it is not the one that we obtained experimentally. Furthermore, taking into account the similarities of the distribution functions obtained with both image-processing methods we truly believe that the results obtained with the black and white contrast method correspond to the conducting elements rather than a rupture process.

In summary, we have assessed the Leonardo’s rule of tree branching from the hydraulic standpoint. We have proposed a model based on the flux conservation, the Hagen-Poseuille law and a xylem elements size-distribution of the Lévy-type. The model tells us that when the Leonardo’s rule is fulfilled the size-distribution presents a power law exponent of three. Furthermore, the total flux is optimized, going from a non-ideal fluid for each xylem element to an idea fluid after summing over all xylem elements. We consider that this optimization process is reasonable because the minerals transported by xylem have to reach all parts of the three to make the photosynthesis process possible. We have corroborated the Leonardo’s rule by directly measuring the stem-branches diameter of four different tree species: Quercus rugosa, Eucalyptus camaldulensis, Schinus molle and Prunus serotina. For the latter we carried out a maceration process in order to unveil the size and number of xylem elements. We found a power law distribution for the xylem elements with an exponent close to three. This value is in close agreement with our theoretical prediction, implying that the seminal hypothesis of Leonardo, the flux conservation is directly related in the trees branching, is valid.