## Abstract

Lognormal distributions and self-similarity are characteristics associated with a wide range of biological systems. The sequential breakage model has established a link between lognormal distributions and self-similarity and has been used to explain species abundance distributions. To date, however, there has been no similar evidence in studies of multicellular organismal forms. We tested the hypotheses that the distribution of the lengths of terminal stems of Japanese elm trees (*Ulmus davidiana*), the end products of a self-similar branching process, approaches a lognormal distribution. We measured the length of the stem segments of three elm branches and obtained the following results: (i) each occurrence of branching caused variations or errors in the lengths of the child stems relative to their parent stems; (ii) the branches showed statistical self-similarity; the observed error distributions were similar at all scales within each branch and (iii) the multiplicative effect of these errors generated variations of the lengths of terminal twigs that were well approximated by a lognormal distribution, although some statistically significant deviations from strict lognormality were observed for one branch. Our results provide the first empirical evidence that statistical self-similarity of an organismal form generates a lognormal distribution of organ sizes.

## 1. Introduction

Modelling individual and ecosystem metabolism is one of the central goals of ecophysiology [1–7]. There have been two distinct approaches to modelling the metabolism (respiration or photosynthesis) of individual plants. One approach is to use biological scaling theories, which focus on allometric (i.e. power-law) relationships between organismal size and metabolism [8–15]. Among biological scaling theories, metabolic scaling theories [8–11], which are based on the self-similarity of plant vascular networks, are one of the most successful approaches because they enable modelling of metabolism from individuals to ecosystems [10,16,17]. The other approach is to use canopy optimization models [18–20], which consider whole plants or stands of plants as a heterogeneous set of terminal organs (e.g. leaves). The optimization models are based on theoretical arguments [18–21] and observations [18,22,23] that terminal organs show plasticity in terms of size, function and direction (e.g. leaf angle), and that this plasticity increases whole-plant photosynthetic rates. The use of optimization models has also been a standard approach for modelling ecosystem carbon gain [19,20].

These two approaches have been equally successful, but to date, the synergism of the two approaches has not been achieved. Previous biological scaling theories, such as metabolic scaling theories, ignore the phenotypic plasticity of terminal organs and instead make the simplifying assumption that the size and function of terminal organs (such as twigs or leaves) are invariant within individual organisms [1,3,4,8,9,13,14,16,17,24,25]. In contrast, optimization models [19,20] consider each set of terminal organs within a given space as a canopy and ignore the fact that those terminal organs are connected to a self-similar resource transportation network and that they are part of an individual plant that obeys allometric relationships. Hence, there has been a discrepancy between metabolic scaling theories and canopy optimization models. Enquist & Bentley [10] have suggested that taking account of the variability of the size of terminal organs will improve present metabolic scaling theories. Smith *et al*. [26], as well as Hunt & Savage [27], have analysed how plasticity of branching morphology affects individual metabolism. However, to date, current biological scaling theories have not successfully modelled plasticity of terminal organs.

The distributions of size in biology have often been approximated by lognormal distribution functions [28–34]. The use of lognormal distributions can be rationalized in part by a mathematical model called the sequential breakage process, which predicts a lognormal size distribution of the end products of a self-similar cascade process [35–38]. In ecology, the sequential breakage of ecological niche spaces has been proposed to explain patterns in the distributions of species abundance [39,40]. To date, however, only indirect evidence has been consistent with this model [40]. Many biological studies that have reported lognormal distributions [28,31–34] have not shown that the processes underlying the distributions were characterized by self-similar geometries. Furthermore, there has been no evidence to support a linkage between lognormal distributions and self-similarity of the forms of individual organisms. Because the relationship between self-similarity and allometry has already been established by biological scaling theories [1,3,4,9,10,14,24], demonstration of a mechanism that connects self-similarity and lognormal distributions would lead to a unified understanding of self-similarity, allometry and lognormal distributions in biology.

Plant vascular networks are a convenient model system for studying self-similarity of organismal forms [8,10,11,13,14,24–26,41,42]. Here, we applied the sequential breakage model of Kolmogorov [35] to self-similar plant forms and hypothesized that the size distribution of the end products of a self-similar process (i.e. the lengths of the terminal twigs of a tree branch) would be approximated by a lognormal distribution. The objective of this study was to provide the first empirical evidence that statistical self-similarity of the form of an individual organism generates a lognormal distribution of the size of its terminal organs.

## 2. The model

### (a) Multiplicative process

The fact that a lognormal distribution can be generated by a stochastic, multiplicative process was first formulized by Gibrat [43]. Let the size of a system (e.g. an organ, organism or population) at time *t* be *X*(*t*). Each system is assumed to grow or shrink by a randomly changing ratio *R*(*t*) [29]:
2.1

The log-transformed size of an object at final time *n* is then expressed as [29]:
2.2

The term ln*X*(0), the logarithm of the initial size, is assumed to be a fixed value. The term ln*R*(*t*) is then decomposed into two factors [44]:
2.3where *M*(*t*) is the mean of ln*R*(*t*) at time *t*, and *ɛ*(*t*) is the deviation of ln*R*(*t*) from *M*(*t*). The log-transformed final size is therefore [44]:
2.4

The second term on the right-hand side of equation (2.4) has the same value for all the final objects and hence does not affect the shape of the final distribution. If values of *ɛ*(*t*) are independent and come from the same distribution function, then the distribution of ln*X*(*n*) will asymptotically approach a normal distribution when *n* is large enough [44]. As discussed by Koch [29], the model described above is the mathematical equivalent of the sequential breakage process proposed by Kolmogorov [35] when *R*(*t*) is the size ratio of an object after a single occurrence of breakage expressed as a fraction of the size before breakage.

### (b) Application to plant form

We applied the sequential breakage process described above to the branching structure of Japanese elm trees (*Ulmus davidiana*), a temperate, deciduous tree species that grows new branches once a year. Following the plant model of Lindenmayer [45], we consider a tree branch to be the result of an iterative branching process. Let us consider a branch that has many twigs. Let the age (denoted as *T*) of the oldest stem be *n* years (*T* = *n*). The stem branches off several younger stems (child stems, age *T* = *n* − 1). Each of these child stems further branches off younger stems (*T* = *n* − 2), and so on to the terminal 1 year old twigs at the branch periphery (*T* = 1).

The stem age described above is a centripetally ordered variable, which starts at the terminal twigs (*T* = 1) and increases towards the oldest stem (*T* = *n*). To make our model consistent with the terminology used in the sequential breakage model, we convert *T* to centrifugal ordering (denoted as *t*), which increases towards the terminal twigs:
2.5

The oldest stem is characterized by order *t* = 0 (the initial single stem segment), and 1 year old stems by order *t* = *n* − 1 (the terminal stem segments). The new ordering is therefore consistent with the sequential breakage process (as is *t* in equations (2.1)–(2.4)), in which one object progressively generates several objects (figure 1). In accordance with the terminology used in river network analysis [46], we regard the order *t* as a scale parameter that represents a sequence of scales within the system. Note that a single child stem (*t* = *a* + 1) connected to a parent stem (*t* = *a*) differs in age and order from the parent stem by 1, unlike the Horton–Strahler order [47], in which the two stems have the same order.

As in the sequential breakage model, the length of a stem of order *t* is denoted by *X*(*t*), and the ratio of its length to the length of its parent is denoted by *R*(*t*) (equation (2.1)). Note that in the present model *R*(*t*) could be greater than one when a child stem is longer than its parent, unlike Kolmogorov's sequential breakage model, in which *R*(*t*) < 1. Fractal-like objects in nature are characterized by statistical self-similarity, which includes undescribed variations or errors at all scales [25,41,46]. We therefore allowed *R*(*t*) to vary among child stems of the same parent by decomposing ln*R*(*t*) into *M*(*t*), the mean of ln*R*(*t*) averaged within each order *t*, and *ɛ*(*t*), the error term that allows variation of lengths among the child stems within each order *t* (equation (2.3)). Statistical self-similarity implies that probability distributions of a stochastic variable as measured over a range of scales are similar to each other [46]. Our first hypothesis was that the branching structure of elm trees is statistically self-similar. This hypothesis implies that the distributions of the errors, *ɛ*(*t*), should be similar at all values of *t*. If this condition is satisfied, then equation (2.4) predicts that the distribution of the sum of *ɛ*(*t*) will approach a normal distribution when *n* is large enough. Our second hypothesis was therefore that the distribution of the lengths of the terminal stems would approach a lognormal distribution as *n* increases (figure 1).

## 3. Material and methods

### (a) Measurements

The sampling site was a cool temperate natural riparian forest on the banks of the Urikari river (42°52′ N, 143°10′ E; elevation: 75 m) in the city of Obihiro in northern Japan. Three young Japanese elm trees (*U. davidiana*) were selected on 9 July 2013. One healthy branch with no significant damage was harvested from each tree (branches ‘1’, ‘2’ and ‘3’); the ages of the branches were 7, 8 and 13 years, respectively. The diameters at breast height of the trees from which the branches were harvested were 6.7, 2.4 and 2.1 cm for branches 1, 2 and 3, respectively. Their lengths (the maximum distance between the branch base and the tip) were 150, 126 and 130 cm, respectively. After harvesting, the branches were air-dried in a ventilated laboratory for three reasons: (i) to prevent rotting and deterioration during storage; (ii) to allow sufficient time for the lengths of the samples to stabilize after harvesting; and (iii) to be consistent with previous analyses of the allometry of tree forms that have used dried materials [48]. The laboratory measurements were conducted from October 2013 to June 2014.

*Ulmus davidiana* is a winter-deciduous species that elongates new shoots only once per year. We determined the order (*t*) of each stem part by counting the number of terminal bud scars from the tip to the base of each branch; a stereomicroscope (SZ61, Olympus, Tokyo, Japan) was used to help discern the oldest, most basal bud scars. With pruning shears or a saw, we decomposed each branch into ‘stem parts’, each of which was of a different order (rectangles in figure 1), while simultaneously recording the connection topology, i.e. all of the parent–child relationships among the stem parts (figure 2). We defined terminal stem parts as new branches that had elongated and matured in the preceding year. These terminal stem parts were 1 year old shoots, excluding the current-year shoots, which were still immature at the time of harvest in early summer. We recognized a stem part as being alive when at least one green leaf was attached to its descendant current-year shoots; only living stem parts were measured. We did not measure small epicormic shoots, which can elongate to fill canopy gaps when regular shoots die [49]. The healthy young branches that we sampled had no large epicormic shoots.

We measured the lengths of long stem parts (more than 10 cm) with a measuring tape and the lengths of short ones with a digital calliper (CD-15CPX, Mitutoyo, Kawasaki, Japan). We used length, rather than diameter, as an indicator of the size of a stem part for two reasons: (i) when the diameter was measured, the value depended on the small pressure of the calliper, especially for the small stem parts (less than 1 mm diameter) and (ii) the cross section of the stem was neither an ideal circle nor an ellipse. Mainly for these two reasons, when we tried to measure the diameter of the same small stem repeatedly, the values we obtained were unstable. This instability would have caused considerable variability in the estimates of the sizes of small stem parts, especially if sizes are recorded on a logarithmic scale. We therefore did not think that diameter was an appropriate metric of size in this study, which focused on the variation of size on a logarithmic scale. In contrast, length was measured without contact between the calliper and a stem part when the calliper was set parallel to the main axis of a stem part. In addition, the length of a stem part was uniquely determined by the distance between the tip (i.e. the terminal bud scar) and the base of the stem part. We measured a total of 1968 stem parts (*n* = 813 from branch 1, *n* = 471 from branch 2 and *n* = 684 from branch 3), including 827 terminal stem parts (*n* = 422 from branch 1, *n* = 200 from branch 2 and *n* = 205 from branch 3; electronic supplemental material, table S1, shows the number of data points at each order). Five small stem parts were lost before they could be measured, and their lengths were therefore not included in the analysis.

### (b) Data analysis

We calculated the ratio of the length of each stem part to the length of its parent, *R*(*t*) (equation (2.1)) and determined the value of *ɛ*(*t*) (equation (2.3)) after calculating *M*(*t*), the mean of ln*R*(*t*) averaged over all order *t* stem parts on each branch. The distribution of *ɛ*(*t*) for each order of each branch was defined as the observed relative frequency distribution of *ɛ*(*t*) pooled over each order on each branch [46]. All statistical analyses were performed with R v. 3.3.1 [50].

### (c) Testing statistical self-similarity

As shown in figure 1 (the model) and electronic supplementary material, table S1 (the data structure), the total number of stems at each order *t* increased exponentially towards the periphery on each branch. There were hence only few stem parts on the proximal positions, whereas there were large numbers of stem parts on the distal positions. According to the self-similar hypothesis, the error distribution curve should converge to a common distribution as more data are incorporated with increasing branching order. If the self-similar hypothesis is true, then each distribution would therefore be similar to the error distribution at the terminal (the most distal) position. We used the bootstrap method proposed by Clauset *et al*. [51] to evaluate the goodness of fit of the error distribution at each order (*ɛ*(*t*)) to the error distribution at the terminal position. First, for each branch, we compared the distribution of *ɛ*(*t*) at each order with the terminal error distribution by using the two-sample Kolmogorov–Smirnov (K–S) test, and we then calculated the *D* statistic (hereafter called *D*_{obs}). The results indicate that there was no significant difference between the *ɛ*(*t*) at each order and the terminal error distribution (*p* > 0.25 for all the branches). Next, random samples were taken from the terminal error distribution to obtain the same sample size as the size of the error distribution at each order by using the R function *sample* with replacement. For each set of samples, a two-sample K–S test was performed to compare the sample distribution with the terminal error distribution, and the resultant *D* statistic (hereafter called *D*_{sim}) was calculated. This procedure was repeated 100 000 times, and the *p*-value was defined as the fraction of *D*_{sim} values that were larger than *D*_{obs} [51]. We found two pairs of data with tied ranks and tested whether the effect of these tied ranks significantly affected the results. Because the true ratio of the lengths, which is a continuous variable, should be different for each pair of stem parts, we added 10^{−12} to one value of *ɛ*(*t*) from each pair of *ɛ*(*t*) values with tied ranks and repeated the above test: we obtained essentially the same results as shown in the Results section.

### (d) Testing lognormality

The two-parameter lognormal distribution is defined by equation (3.1) [44]:
3.1where *x* is the length of a terminal stem part (mm), *f*(*x*) is the density of terminal stem parts with length *x*, *N*_{total} is the total number of terminal stem parts on each branch and *m* (mean of ln *x*) and *σ* (standard deviation of ln *x*) are the curve-fitting parameters. The parameters (*m* and *σ*) in equation (3.1) were estimated by maximum-likelihood estimation (MLE) by using the R package *fitdistrplus* [52]. The cumulative distribution function (CDF) of the lognormal distribution is given by equation (3.2) [34]:
3.2where *N*(*x*) (1 ≤ *N*(*x*) ≤ *N*_{total}) is the cumulative number of all terminal stem parts as long as or longer than *x*. The function erf(*z*) is the Gauss error function and is defined by equation (3.3) [34]:
3.3

We then determined whether each lognormal distribution function and normal distribution function was significantly better than the other in terms of goodness of fit to our dataset by using the Vuong likelihood ratio test [51,53]. The CDF of the normal (Gaussian) distribution is given by equation (3.4): 3.4

The parameters of the normal distribution function (*m*′ and *σ*′) were also estimated by MLE. We then calculated the log-likelihood ratio between the lognormal and normal distributions and tested whether the log of the ratio was significantly different from zero by using the equations described in Clauset *et al*. [51].

Next, after log-transformation of the data, several normality indices of the empirical distributions were compared with the theoretical normal distribution by using Monte Carlo methods. For each branch, the R function *rnorm* was used to generate normal random numbers with the same mean and variance as the empirical distribution and with the same sample size as the experimental data. The distributions of normality indices were then calculated (i.e. the third- and fourth-order moments, skewness, kurtosis and *D* statistic for the Lilliefors normality test [54], which in this case was equivalent to the goodness of fit test of Clauset *et al*. [51]. For each branch, this simulation was repeated 100 000 times, and the distributions of these indices were calculated (hereafter called the simulated distributions of the indices). For each branch, we tested whether each index of the empirical data was within the 95th percentile of the simulated distribution of the index. Several terminal stem parts on the same branch had exactly the same length (e.g. 2.23 mm) because of the limited resolution of the digital calliper (0.01 mm). We made the same adjustment as described above for these tied ranks, and we obtained essentially the same results as shown in the Results section.

## 4. Results

Our first hypothesis, that the branching structure of elm trees is statistically self-similar, was supported by the similarity of the observed empirical cumulative distributions of *ɛ*(*t*) for different values of the scale parameter *t* (figure 3). The results of Clauset *et al*. [51]'s goodness-of-fit tests showed that there were no significant difference between the error distributions at each order and the terminal error distribution (*p* > 0.17 in all cases, see electronic supplemental material, table S2 for the *p*-value of each test). We tested whether the distribution of the *p*-values that we obtained differed significantly from a uniform distribution by using the two-sample K–S test and the R function *punif*. For this test, all the *p*-values from the three branches were pooled (*n* = 23). The result was not significant (*p* = 0.41), the indication being that the *p*-values were not significantly biased.

Our second hypothesis, that the distribution of the lengths of the terminal stems would approach a lognormal distribution, was also supported by the good approximation of the distributions of the lengths of terminal stem parts to a lognormal distribution (figure 4). The Vuong likelihood ratio test showed that a lognormal distribution function (the red curve in figure 4) was significantly preferred over a normal distribution function (the blue curve in figure 4; *p* < 0.01 for all the branches). However, after log-transformation of the data, we found significant deviations from a theoretical normal distribution on the basis of two statistical tests for branch 1 (table 1), which had the lowest terminal branching order (*t* = 6) with the largest number of terminal stems (*n*_{total} = 422). For branch 2 (*t* = 8) and branch 3 (*t* = 12), which had larger terminal branching orders, the distributions of the log-lengths were not significantly different from theoretical normal distributions on the basis of any normality test. These results agree with the model prediction that the distribution of the terminal twig lengths should approach a lognormal distribution as the terminal branching order increases.

## 5. Discussion

The results are consistent with the assertion that the within-scale error distribution is often approximated by a lognormal rather than normal distribution in biological allometry [30,33]. However, a strictly lognormal distribution would be approached only if the terminal branching order (*t*) approached infinity. In the case of such a distribution, the probability of observing some very small and some very large stems would increase as the sample size increases. Because this possibility is biologically unrealistic, there might be no strictly normal or lognormal distribution in real biological data if the sample size is sufficiently large and the order finite. In the present case, the terminal branching order of branch 1 (*t* = 6) may not have been large enough to satisfy the assumption of the model. In addition, in the case of elm trees, our results seem to indicate that there was a lower limit to terminal stem size; the shortest terminal stem lengths on each branch were 1.30, 1.59 and 1.23 mm for branches 1, 2 and 3, respectively. That shorter stems were not found might reflect the fact that the smaller stems did not survive or that smaller buds became dormant. The observed deviation from a theoretical lognormal distribution in branch 1 may have been caused by these biological systematic deviations from strict lognormality.

Many empirical studies have demonstrated that variations in diverse biological phenomena are well approximated by lognormal distributions [28–34], but in those studies, rigorous statistical tests were not performed. Those distributions may not have been theoretical lognormal distributions in the strict sense, even if the processes that underlay them were multiplicative. In this study, we proposed a lognormal distribution as a theoretically motivated simple model, because the underlying process was statistically self-similar and multiplicative. Including detailed biological factors, such as the limits of organ size, would further improve the predictive power of the model in future studies.

These results constitute an idealized approximation that we propose as a starting point. Plant form in nature is affected by wind, herbivory and competition with neighbours [55,56]. These factors can be expected to cause deviations from ideal self-similarity and hence from lognormality. Further study is therefore needed to improve the analytical model before it can be applied in real-world situations.

## Data accessibility

All the raw data are available from the Dryad Digital Repository: doi:10.5061/dryad.776ht [57].

## Authors' contributions

K.K. conceived of the study, designed the study and collected the data. All authors analysed the data and were involved in writing the manuscript.

## Competing interests

We have no competing interests.

## Funding

This work was supported by Japan Society for the Promotion of Science grant in aid for Scientific Research number 25891001.

## Acknowledgements

We thank the English Resource Center of the Obihiro University of Agriculture and Veterinary Medicine for help with proofreading this paper.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3649904.

- Received November 1, 2016.
- Accepted December 5, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.