Statistical hydraulic model for the Leonardo's rule
O. Sotolongo-Costa, P. Villasana-Mercado, L. Sánchez-Calderón, I. Rodríguez-Vargas
SStatistical hydraulic model for the Leonardo’s rule
O. Sotolongo-Costa ∗ Centro de Investigaci´on en Ciencias,Instituto de Investigaciones en Ciencias B´asicas y Aplicadas,Universidad Aut´onoma del Estado de Morelos,Av. Universidad 1001, Col. Chamilpa,62209 Cuernavaca, Morelos, Mexico.
P. Villasana-Mercado
Unidad Acad´emica de F´ısica, Universidad Aut´onoma de Zacatecas,Calzada Solidaridad Esquina con Paseo la Bufa S/N, 98060 Zacatecas, Zac., Mexico.
L. S´anchez-Calder´on
Unidad Acad´emica de Ciencias Biol´ogicas,Universidad Aut´onoma de Zacatecas,Calzada Solidaridad Esquina con Paseo la Bufa S/N, 98060 Zacatecas, Zac., Mexico.
I. Rodr´ıguez-Vargas † Unidad Acad´emica de Ciencia y Tecnolog´ıa de la Luz y la Materia,Universidad Aut´onoma de Zacatecas,Carretera Zacatecas-Guadalajara Km. 6,Ejido La Escondida, 98160, Zacatecas, Zac., Mexico (Dated: September 23, 2020) a r X i v : . [ phy s i c s . b i o - ph ] S e p bstract More than five hundred years ago Leonardo Da Vinci found a pattern in the growth of treesnowadays known as the Leonardo’s rule. This rule relates the thickness of the stem with thethickness of the branches at different bifurcation stages in a Pythagorean fashion. He argued thathis rule was the result of the conservation of sap flux. In the present work, we explore this idea byassuming that the sap flux through each xylem element behaves as a non-ideal fluid and the size-distribution of the xylem elements obeys a power law distribution. We find that the simultaneousfulfillment of Leonardo’s rule and the conservation of the sap flux, lead to a global behavior of thesap like that of an ideal fluid after summing over all xylem elements. These results are supportedby field and experimental work. In particular, we corroborated the Leonardo’s rule in differenttree species by measuring the stem-branches thickness at different bifurcations stages. We alsodetermined the statistical size-distribution of the xylem elements through a maceration process,finding that it corresponds to a power law distribution with an exponent close to three, which isthe exponent required for the Leonardo’s rule. As far as we know this is the first time that astatistical hydraulic model supported by experimental data is presented for the Leonardo’s rule.
INTRODUCTION
Leonardo da Vinci was a man with innate curiosity and exceptional intellectual capacitythat allowed him to understand nature [1]. Leonardo is best known for his iconic paintingsthe 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 onobservation and experience he deciphered different patterns in nature. In plants he unveiledthe so-called sixth leaf, rings of growth and branching patterns in trees nowadays known asLeonardo’s rule [1, 3]. Regarding the latter, he found that the branches at different heightsare proportional to the stem. In specific, the sum of the transverse areas of the branches atdifferent 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 thisrelationship to the conservation of sap flux [1]. The modern version of the Leonardo’s ruleis given in Pythagorean fashion d ∆ = (cid:80) ni =1 d ∆ i , where d is the diameter of the stem, d i is thediameter of the i -th branch at a particular bifurcation level, n the number of branches at2he mentioned level and ∆ the so-called Leonardo’s exponent [4]. The nominal value of theLeonardo’s exponent is 2, however it is well known that it depends on the tree species [4–6].The Leonardo’s rules is also obeyed by the radical system of trees [7]. Furthermore, theLeonardo’s rule in conjunction with the external fractal architecture of trees is involved inthe optimization of different biological processes [8, 9]. For instance, the modeling of treesshows that the mechanical stress induced by wind loads is minimized by the self-similarorganization of the branches as well as the area-preserving condition [8]. Likewise, thecompetition for light is optimized due to the complex-fractal growth of trees in which theLeonardo’s rule is implicated [9].There are also several works assessing the impact of hydraulics in plant vascular systems[5, 10–13]. From the rather simple pipe model to the most sophisticated proposals basedon internal and external fractal architectures. The common factor in all these studies is thederivation of well-known allometric scaling laws. In the case of the pipe model [5, 10, 11],trees are modeled as a uniform tube from the base to the top, with the tube diameterequivalent to the cross-section of the stem. The pipe model entails the Leonardo’s rule aswell as the non-linear scaling relationships for the height and mass of a plant with respect tothe basal stem diameter. The pipe model has also been useful to determine the biomass ofplants. A historic review of the pipe model can be consulted in [11]. The WBE model is amore refined proposal in which a general theory of resource distribution through hierarchicalbranching networks is considered [12]. In particular, the model takes into account externalspace filling (self-similar branching), hydraulic optimization (tube tapering) and scale freecanopy (invariant leaf-petiole size). In this model despite tubes taper the flow rate is sizeindependent due to the increase of the flow velocity in small radius branches. The WBEmodel also shows that the well-known allometric scaling laws for animals [14] apply forplants [15], supporting the idea that they are universal in biology and can be modeled bysimple geometric and hydrodynamic principles [16]. A decade later of the WBE modela theory that also considers internal space filling as well as trade-offs between hydraulicsafety and efficiency was proposed [13]. The interrelationship between the vascular andbranching networks give rise to more general allometric scaling laws that can respond to theparticularities of tree species. In addition, the model gives predictions for the sap (solutesand water) flow, the taper and frequency of xylem conduits. In contrast to the WBE modelthe flow velocity remains constant regardless of the tapering of the xylem conduits, and the3ow rate size independent as a result of the internal fractal architecture. The area preservingcondition is also fulfilled in this model. As we have documented, it is well accepted that theexternal branching network minimizes the mechanical stress in trees, optimizes the captureof 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 modelwe can have different outcomes, as in the case of the fluid velocity. Finally, it is important toremark 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 theseconditions 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 sapflux is directly implicated in the Leonardo’s rule. We propose a statistical hydraulic modelbased on a power law distribution for the xylem elements. We corroborate the Leonardo’srule by measuring the radius of the stem and branches of different tree species. We alsoobtain the size distribution of the xylem elements after subjected one of the species to amaceration process. A fairly good agreement is obtained between the theoretical predictionand the experimental results for the size distribution of the xylem elements. STATISTICAL HYDRAULIC MODEL
As we have documented the hydraulic model has been modified and refined throughoutyears. From the rather simple pipe model to the more elaborated proposals of the lastdecades. In the present work, we adopt an statistical hydraulic version based on a powerlaw distribution for the xylem elements that accounts of the space filling internal fractalarchitecture of trees as well as memory-correlation effects of the conduction process. Wealso consider that the flux through each xylem element is non ideal due to the intricaciesrelated to the internal architecture of the conduction vessels [17, 18]. In addition, we assumethe Leonardo’s seminal hypothesis as valid, that is, the Leonardo’s rule is a consequence ofthe sap flux conservation [1]. In fact, by attending these requirements we will show that thepower law exponent is three and that the sap flux is optimized (ideal fluid behavior) as aconsequence of the collective behavior of the xylem elements.Our starting point is the mentioned Leonardo’s hypothesis [1]. In specific, if the stem isbifurcated in two daughter branches we can established a relationship between their corre-4ponding fluxes as follows Q = Q + Q , (1)where Q is the flux of the main branch and Q and Q are the fluxes of the daughter branches.Now, by assuming that the flow through individual xylem elements is of non ideal nature,that is, it obeys the Hagen-Poseuille law [17, 18], it is possible to establish a relationshipbetween the flux ( q ) and radius ( r ) of the xylem elements as, q ( r ) ∝ r , (2)where the proportionality constant (not shown) comes in terms of intrinsic characteristics ofthe fluid. For the moment, we consider that the constant is not essential for our derivation.The effective flux through branches can be computed by averaging over all xylem elements Q ( R ) = (cid:90) R f ( r ) q ( r ) dr, (3)where f ( r ) represents the size distribution of xylem elements. By considering a power lawdistribution f ( r ) = r − x , (4)we can obtain a relationship between the flux and radius of a branch as Q ( R ) ∝ R − x , (5)with x an exponent to be determined. By substituting Eq. (5) in Eq. (1) we obtain R − x = R − x + R − x , (6)where R is the radius of the main branch, and R and R are the radii of the daughterbranches.Finally, according to the Leonardo’s rule the radii of tree branches are related in Pythagoreanfashion R ∆ = R ∆1 + R ∆2 . (7)5onsidering that the Leonardo’s exponent is approximately ∆ ≈
2, and by comparing Eqs.(6) and (7) we obtain a power law exponent x ≈ FIELD WORK: MEASURING STEAM AND BRANCHES
The first task that we carried out it was the corroboration of the Leonardo’s rule bydirectly measuring the stem-branches radii of four different trees species. We measured thestem-branches of the first bifurcation of
Quercus rugosa , Eucalyptus camaldulensis , Schinusmolle and
Prunus serotina . We have followed the measuring protocol of Aratsu [5]. Inparticular, we carefully choose healthy trees, with no damage in the branches and with anevident bifurcation. In addition, we have measured the stem and the daughter branchesat a reasonable distance from the bifurcation to have reliable data, particularly, to avoidthe possible data dispersion associated to the natural branch thickening around the bifur-cation. We have considered 91, 50, 48 and 58 individuals for
Quercus rugosa , Eucalyptuscamaldulensis , Schinus molle and
Prunus serotina , respectively.In order to corroborate the Leonardo’s rule we will assume that the perimeter of the stemand daughter branches corresponds to the one of a circle. So, the Pythagorean relationshipbetween radii Eq. (7) is also valid for the perimeters P ∆1 P ∆ + P ∆2 P ∆ = 1 . (8)Once the perimeters are measured we can compute the Leonardo’s exponent ∆. This processis performed for each individual, to finally obtain a net Leonardo’s exponent by averagingover all individuals of each tree specie. In the case of Quercus rugosa , Eucalyptus camal-dulensis and
Schinus molle the perimeters are determined with a measuring tape of onemeter length and an experimental error of ± . Prunus serotina we have used a digital Vernier with a measuring range of 150 mm and aprecision of 0.1 mm. In this case the individuals are young trees. For more details aboutthe characteristics of the trees see the supplementary information.6
XPERIMENTAL WORK: MACERATION PROCESS
The second task that we performed was to unravel the trees xylem elements through theso-called maceration process [19]. This process is recommended for young trees, thus, weimplemented it only in
Prunus serotina . In specific, we macerated three individual accordingto the following protocol:1. For each individual we prepare a 300 mL maceration solution, of which 120 mL aredistilled water, 30 mL are hydrogen peroxide and the remaining 150 mL correspondto glacial acetic acid.2. With a handsaw we take pieces of the tree right above and below the bifurcationsby cutting the branches transversely. In total we collect five pieces, three at the firstbifurcation stage and two more at the second stage. From each piece we extract fineslices with a scalpel.3. The fine slices of each piece are placed in a sterile flask with 60 mL of macerationsolution. The five flasks of each individual are sealed and introduced in an incubationoven for a five day period at a temperature of 60 ◦ C.4. After the time in the incubation oven each flask sample is subjected to three rinses inorder to drain the glacial acetic acid. The remaining material is let to rest a wholenight in distilled water.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 todo 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, protectingthem with a cover glass. Once the sample is ready it is placed in an optical microscope under100X magnification. We then proceed to scan from left to right and from up to down theobserving area by taking pictures of it in the process, see Fig. 1. The pictures are processedwith the software ImageJ. We consider two methods to obtain the size distribution of thexylem elements. One in which we measure directly the diameter of the xylem elements and7
IG. 1. Images of the xylem elements after the maceration process. (a) 10X magnificationof a region of the Neubauer chamber, in which we can appreciate the xylem elements dispersedthroughout the Neubauer chamber grid. (b)The same as (a), but under a 40X magnification. Atthe center we can see xylem structures known as xylem vessels [20]. other in which the images are black and white contrasted to measure the shaded areas of thexylem elements. More details of these methods can be found in the supplementary material.
RESULTS
In this section we present the main results of the field and experimental work. Mostimportantly, we contrast the resulting size distribution of the xylem elements after themaceration process with the statistical hydraulic model outcomes.Regarding the field work, in table I we show the Leonardo’s exponent of the tree speciesstudied. The exponents are obtained by computing the zeros of Eq. (8) and averaging theexponents of all individuals of each tree specie. The standard error of the exponents is alsoobtained and shown in table I. As we can see the Leonardo’s exponent for all tree species liesin the well-known range 1 . < ∆ < . Prunus serotina ) is below2. As the Leonardo’s exponent depends on the tree specie, we also expect that the sizedistribution of xylem elements does. In fact, in table I we also show the expected exponentsfor 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 ≈ .
11 in the experimental work.8
ABLE I. Leonardo’s exponent of the tree species analyzed in the field work. The expectedexponent for the size distribution of the xylem elements is also included for each specie.Regular name Scientific name Leonardo’s exponent Size distributionexpected exponentPirul
Schinus molle . ± .
08 2 . ± . Eucalyptus camaldulensis . ± .
10 2 . ± . Quercus rugosa . ± .
06 2 . ± . Prunus serotina . ± .
07 3 . ± . Regarding the experimental work, we firstly show the results of measuring directly thediameters of the xylem elements after the maceration process. The size distributions for thethree
Prunus serotina individuals are shown in Fig. 2. As we can see the size of the xylemelements ranges from 8 µ m to 30 µ m. Furthermore, small size xylem elements are the vastmajority, resulting in a power law distribution for all individuals. In fact, the best fittingfunction for the experimental data is of the form f ( r ) = ar − b , see the solid-blue curves inFig. 2. In all cases the correlation coefficient is over 90 %, indicating a fairly good fitting forthe power law function. We also determine the relative error of the resulting exponent withrespect to the expected one. The resulting exponents for the three individuals as well as itscorresponding relative errors are shown in table II. As we can notice the exponent of thethird individual ( b ≈ .
13) is closest to the expected one ( x ≈ . Prunus serotina individuals. So, the exponents that we obtained are the result of processinga reduced number of xylem elements images.In order to overcome this obstacle, we have used an alternative method based on blackand white contrast the xylem elements images. This method allows us to measure in anautomatic way the shaded areas of the xylem elements and to determine indirectly thecorresponding 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. Aswe can notice the size of xylem elements ranges from 20 µ m to 180 µ m, being the vast9 IG. 2. Size distribution of the xylem elements of the three
Prunus serotina individuals. Thesize distributions were determined by directly measuring the diameters of the xylem elements. Thesolid-blue curves correspond to the power law function fitting. b represents the exponent of thepower law function f ( r ) = ar − b . The correlation coefficient (cc ) for each individual is also shown,being best for the first individual. In the three cases about 300 xylem elements were measured.TABLE II. Size distribution exponent for the Prunus serotina individuals by measuring directlythe diameters of the xylem elements. The exponent relative error is also presented. The exponentis calculated by fitting the experimental data to the power law function f ( r ) = ar − b . The relativeerror is computed by comparing the exponent obtained in the experimental work ( b ) with theexpected one ( x ) according to the field work and the statistical hydraulic model: (cid:15) = b − xx × b ) Relative error ( (cid:15) )First 3.26 4.80 %Second 2.88 7.39 %Third 3.13 0.64 % majority small size xylem elements. The experimental data fit quite well to a power lawfunction (solid-blue curve), with a size distribution exponent b ≈ .
08 and a 99 % correlationcoefficient. In addition, the relative error of the size distribution exponent is 0.96 %.10
IG. 3. Size distribution of the xylem elements of the first
Prunus serotina individual by measuringthe areas of the xylem elements according to the black-white contrast method. The solid-bluecurve correspond to the power law function fitting. In this case about 10,000 xylem elements weremeasured.
DISCUSSION
The first aspect that we want to discuss is related to the characteristics of the xylemelements. It is well known that as trees grow some xylem elements lose its conductiveability, as they fill with resin, giving place to the so-called heartwood [17]. At the sametime new xylem elements arise, keeping the tree conducting capability. The elements thatcompose the xylem have in its walls lignin [17, 21], a substance that makes them rigid11nd woody, providing mechanical stability and support to the tree. Both conducting andnon-conducting elements contribute to the thickness of the stem and branches. With themaceration process what we obtain is a disperse sample with conducting and non-conductingelements. Our first image processing method allows us to choose xylem structures withconducting capabilities such as tracheids and vessel elements. With the black and whitecontrast 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 thedetermination of the size distribution because they contribute to the thickness of the stemand branches. Under this context, we expect results in better agreement with respect to thepredicted ones with the black and white contrast method. In fact, this is what happenedwith the first individual. After considering conducting and non-conducting elements theexponent of the size distribution ( b ≈ .
08) results in close agreement with respect to thepredicted one ( x ≈ .
11) with a relative error of 0.96 %.The second aspect that we want to address it is about the type of distribution functionthat was obtained. In specific, a power law distribution function alike to the so-called L´evy-type distribution functions [22–24]. In fact, these distribution functions are characterized bymemory and correlation effects. Taking into account that xylem structures form a complexconducting network [17], in which memory and correlation effect are plausible, we considerthat it is totally reasonable to obtain a power law distribution for the xylem elements. Infact, a power law distribution, in which small size xylem elements dominate, optimized theconduction process by changing the non-ideal fluid behavior of individual xylem elements toan ideal one after summing over all xylem elements. Actually, at first sight, it is counterintu-itive 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 ratherthan its individual behavior. This intricate organization of the xylem elements it is likelyto be related to the space-filling maximization as in the case of other biological complexnetworks [25–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 smallradius 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-processingmethod. In principle, it is possible that with the maceration process a lot of conductingstructures be broken. If so, what we are obtaining with the black and white image-processing12ethod would be the distribution function of a rupture process [29, 30], as we are measur-ing automatically all structures in our sample. However, the exponent of the distributionfunction associated to a rupture process it is well known [29, 30] and it is not the one thatwe obtained experimentally. Furthermore, taking into account the similarities of the dis-tribution functions obtained with both image-processing methods we truly believe that theresults obtained with the black and white contrast method correspond to the conductingelements rather than a rupture process.
CONCLUSIONS
In summary, we have assessed the Leonardo’s rule of tree branching from the hydraulicstandpoint. We have proposed a model based on the flux conservation, the Hagen-Poseuillelaw and a xylem elements size-distribution of the L´evy-type. The model tells us that whenthe 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 elementto an idea fluid after summing over all xylem elements. We consider that this optimizationprocess is reasonable because the minerals transported by xylem have to reach all parts ofthe three to make the photosynthesis process possible. We have corroborated the Leonardo’srule by directly measuring the stem-branches diameter of four different tree species:
Quercusrugosa , Eucalyptus camaldulensis , Schinus molle and
Prunus serotina . For the latter wecarried 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 seminalhypothesis of Leonardo, the flux conservation is directly related in the trees branching, isvalid.P.V.-M. acknowledges to CONACYT-Mexico for the scholarship for graduate studies. ∗ [email protected] † [email protected][1] F. Capra, Learning from Leonardo (Barrett-Koehler Publishers, 2013).[2] J. Pevsner, Trends Neurosci. , 217 (2002).
3] H. Suh,
Leonardos Notebook (Black Dog, 2005).[4] B. Mandelbrot,
The Fractal Geometry of Nature (Freeman, 1982).[5] R. Aratsu, J. Young Investigators (1998).[6] K. Sone, A. A. Suzuki, S. I. Miyazawa, K. Noguchi, and I. Terashima, J. Plant Res , 41(2008).[7] L. Armin, K. Winfriend, and L. Douglas, Tree Physiol. , 117 (2000).[8] C. Eloy, Phys. Rev. Lett. , 258101 (2011).[9] C. Eloy, M. Fournier, A. Lacointe, and B. Moulia, Nat. Commun. , 1014 (2017).[10] K. Shinozaki, K. Yoda, K. Hozumi, and T. Kira, Jpn. J. Ecol. , 97 (1964).[11] L. Lehnebach, R. Beyer, V. Letort, and P. Heuret, Annals of Botany , 773 (2018).[12] G. B. West, J. H. Brown, and B. J. Enquist, Nature , 664 (1999).[13] V. M. Savage, L. P. Bentley, B. J. Enquist, J. S. Sperry, D. D. Smith, P. B. Reich, and E. I.von Allmen, Proceedings of the National Academy of Sciences , 22722 (2010).[14] K. Schmidt-Nielsen, Scaling: Why is Animal Size so Important? (Cambridge Univ. Press,1984).[15] K. J. Niklas,
Plant Allometry: The Scaling of Form and Process (Univ. of Chicago Press,1994).[16] G. B. West, J. H. Brown, and B. J. Enquist, Science , 122 (1997).[17] M. T. Tyree and M. Zimmermann,
Xilem Structure and the Ascent of Sap (Springer, 2002).[18] L. D. Landau and E. M. Lifshitz,
Fluid Mechanics (Elsevier Ltd., 1987).[19] R. Fonnegra and J. Santa, Actual Biol , 25 (1978).[20] P. H. Raven, R. F. Evert, and S. E. Eichhorn, Biology of Plants (W. H. Freeman, 2004).[21] L. A. Donaldson, Phytochem. , 859 (2001).[22] K.-I. Sato, Levy Process and Infinitely Divisible Distributions (Cambridge University Press,1999).[23] O. Barndorff, T. Mikosch, and S. Resnick,
Levy Processes Theory and Applications (SpringerScience, 2001).[24] W. Feller,
An Introduction to Probability Theory and Its Applications (John Wiley & Sons,1968).[25] G. M. Viswanathan, E. P. Raposo, and M. G. E. da Luz, Physics of Life Reviews , 133 (2008).
26] J. F. Gillooly, J. H. Brown, G. B. West, V. M. Savage, and E. L. Charnov, Science , 2248(2001).[27] N. E. Humphries, H. Weimerskirch, N. Queiroz, E. J. Southall, and D. W. Sims, PNAS ,7169 (2012).[28] P. R. Larson,
The Vascular Cambium (Springer, 1994).[29] O. Sotolongo-Costa, E. Lopez-Pages, F. Barreras-Toledo, and J. Marin-Antua, Phys. Rev. E , 4027 (1994).[30] O. Sotolongo-Costa, Y. Moreno-Vega, J. J. LLoveras-Gonz´alez, and J. C. Antoranz, Phys.Rev. Lett. , 42 (1996)., 42 (1996).