Looking into the Natural Log: A Mathematical Analysis of Rubber Tree Hevea Brasiliensis

Keywords: fractal scaling, reaction-diffusion, growth rate, mating system, latex flow rate

Abstract

Within the structure of Hevea brasiliensis lies a quiet yet meaningful conversation between mathematics, geometry and its survival. By taking a deeper look into these aspects, the way in which H. brasiliensis has evolved to thrive in its environment can be revealed. This essay firstly explores the patterns in vascular tissue which ultimately helps in structural support and transport. Next, the growth rate of H. brasiliensis is examined to understand its optimal height for efficient use of its energy. This is followed by an analysis of the growth of the lateral roots, first in relation to leaf growth which results in the tree allowing for better nutrient and water uptake, and second, looking at the fractal character of the roots which optimizes the root-soil interface. Furthermore, the mating system of H. brasiliensis is examined which leads to the discovery of the tree primarily undergoing cross-pollination which leads to an increase of genetic diversity within the species. Finally, the latex flow rate is examined in terms of the turgor pressure, responsible for latex exudation, and the rate of wound plugging, which is found to be greater towards the end of flow, highlighting how latex flow is a two-phase mechanism ensuring the efficient protection of the tree.

Introduction

Nature’s adaptations often follow mathematical principles, from game theory explaining interactions between different organisms, to the few variable changes in reaction-diffusion mechanisms that results in numerous pattern formation. Hevea brasiliensis is perhaps best known as the world’s primary supplier of natural rubber. In fact, the mathematical principles apply to all aspect of H. brasiliensis’ survival. The reaction-diffusion mechanisms form optimally oriented vascular tissue pattern, the time-dependent function of latex exudation protects the tree, the cross and self-pollination mechanisms ensure genetic diversity, and the fractal scaling of its root provides maximum surface area for nutrient uptake.

Among different mathematical concepts, perhaps the most universal in living organisms is the reaction-diffusion mechanisms, which accounts for many natural patterns. This phenomenon is first described by Alan Turing in 1952, as he put forward the foundation of the mechanism (Siero & Deinum, 2025). He proposes that the pattern arises from the combined effects of self-activating activators and the inhibitors that suppress them. Initially, Turing’s theory primarily described pattern formation like the zebra’s stripes (Siero & Deinum, 2025). Later, as understanding in genetic and hormone roles regulation deepened in plant development, this mechanism extended to many more functional patterns in plants, including the arrangement in roots, phyllotaxis and leaf venation. In H. brasiliensis, it is responsible for the distribution of vascular tissues.

Pattern Formation in Vascular Tissue of Hevea brasiliensis

Context on reaction-diffusion mechanism

Hevea brasiliensis’ vascular distribution pattern is attributed to the reaction-diffusion mechanisms, which facilitates differentiation of xylem, phloem and cambium instead of random distribution. Turing’s earliest model of reaction-diffusion highlights the absence of a stable, spatially homogeneous steady state (where variables are constant or tend to return to equilibrium when perturbed). In his models, diffusion destabilizes uniformity and thus drives the emergence of spatial patterns. At its simplest, the reaction-diffusion mechanism is a two-component system. Specifically, a combination of short-range activation and long-range inhibition monitored by different diffusion coefficients. Of the two, the slowly diffusing component is self-activating, amplifying the effect of local disruptions and promoting autocatalysis/activation. Meanwhile, the larger diffusion coefficient is inhibiting, reducing the effect of disruptions and attempting to return to an equilibrium (Siero & Deinum, 2025). 

Pattern formation of vascular tissue

This framework applies to plant vascular development, modeling the differentiation of procambium, phloem, and xylem from initially undifferentiated cell (Carteni et al., 2014). The model shows that the mechanism can reproduce the range of radial vascular arrangements observed in plants, from simple protosteles (solid core of vascular tissue) to complex eusteles (vascular bundles arranged in a core) like H. brasiliensis. Typically, vascular tissue distribution falls into two categories: a longitudinal pattern that describes an array of vascular bundles in different organs, or a radial pattern that describes the arrangement of xylems in the centre, cambium in between, and phloem near the periphery. This is a result of reaction-diffusion mechanism on cell differentiation. In plants, this process begins with auxin signaling during embryogenesis. Afterwards, two morphogens, centrifugally diffusing auxin and centripetally diffusing sucrose, formulate the pattern of new cambium tissues. Essentially, a concentration gradient of these two morphogens within the meristem influences procambium formation. Additionally, phloem and xylem always occur together in a fixed, coordinated pattern. The model thus hypothesizes that phloem and xylem are locally inhibiting yet laterally enhancing. This means a single cell differentiates into either xylem or phloem locally, while the facilitate the other’s development laterally (Carteni et al, 2014). 

The formation of different domains is described by the partial differential equation (Eq. 1) involving a morphogenetic signaling factor (a signal molecule acting over a long distance):

Where S₀ = the first morphogenetic signaling factor, d = distance from centre of section, μ_S0 = constant consumption rate, σ_S0 = basic production rate, r = radius of the domain, and D_S0 = diffusion coefficient (Carteni et al., 2014).

Meanwhile, additional diffusible signals of procambium interact through two diffusible substrates, S₁ and S₂ and a diffusible, autocatalytic procambium activator A_p to initiate differentiation. S₁ and S₂ are activated under different conditions: S₁ is produced when a threshold value (S^*₀) is not met, while S₂ is produced when S₀ ≥ S^*_0. As a result, S₁ is only produced in central domains, and S₂ in peripheral regions (Carteni et al, 2014).

The activation of S₁, S₂, and A_p is modelled as Eq.2, Eq. 3 and Eq. 4:

The parameters used in equations 2-4 are listed in Table 1:

Table 1. Different parameters used in equation 2-4 (Carteni et al., 2014).

The two unlisted parameters, σ_S1 and σ_S2 are spatial variables that monitor substrate production in specific locations. The basic production rate of the substrate S₁ and S₂ is spatially defined. σ_S1 is 1 (trigger production of S₁ in central region when the threshold value S^*₀ is met), 0 if threshold is exceeded. Opposingly, σ_S2 is only triggered in the peripheral region when the threshold value is exceeded, and untriggered in central region. When the activator A^*_p reaches threshold value (A^*_p), procambium differentiation occurs, otherwise the cell does not differentiate into procambium tissue (Fig. 1) (Carteni et al., 2014).

Fig.1. Schematic representation of how whether the S₀ (first morphogenetic factor) reaches the threshold value S₀^* determines a clear boundary formation between central and peripheral regions. This is also shown through the peak of A_p (procambium activator) is most activated and establishes the boundary (Carteni et al, 2014). 

Like procambium differentiation, two autocatalytic activators A_F and A_x, promote the differentiation of xylem and phloem. Unlike the cambium tissue formation, these factors are locally inhibitory through a common repressor R, yet they mutually activate each other over long distances because each produce the other’s signal/substrate (S_F and S_x) (Fig. 2) (Carteni et al, 2014).  

Fig. 2. (A) Interaction network of phloem activator A_F and xylem activator A_{x as} they autocatalyze their own production, but inhibitor the other activator, and the enhance each other over long distance through their diffusible substrates S_F and S_x . (B) Diagram of resulting spatial concentrations from interactions of activators locally (inhibitory) and long-range (enhancing) relative to the differentiation threshold A^*, resulting in formation of distinct xylem and phloem regions (Adapted from Carteni et al, 2014). 

Here, the dynamics of different activators and substrate interactions are modelled by Eq. 5, Eq. 6, Eq. 7, Eq. 8:

Numerous parameters are used to model both where vascular tissues differentiate and how their local inhibition occurs through A_X - S_F and A_F – S_X. The complete list of parameters is listed below (Table. 2)

Table. 2. List of parameters in equations 5-8 (Carteni et al., 2014).

Among the most notable factors, the spatial variables that shape phloem and xylem tissues pattern follow the constraint (Eq. 9, Eq. 10):

indicating that xylem and phloem differentiation occurs only if procambium tissues have been differentiated. In addition to this, S_F is the spatial variable representing the basic production of phloem (S_F) in the peripheral region, while S_x shows the basic production of xylem (S_x) in the central region. In the context of H. brasiliensis, these different mechanisms come together to form the collateral bundle pattern found in the tree’s vascular bundle (Fig. 3) (Sae-Lim et al., 2019; Carteni et al, 2014). 

Fig.3. Schematic representation of common arrangements of vascular tissue in H. brasiliensis, with a clear boundary between different vascular tissue, and an organized layout of xylem, cambium and phloem in pattern of collateral bundle pattern (Sae-Lim et al., 2019). 

The mechanism that causes this specific pattern formation of vascular bundles allows H. brasiliensis to optimize the function of each tissue. Xylem cells located in the centre of the stem provide both better structural support, and transport water through the tree. Phloem, positioned near the periphery of the vascular bundles, ensures the most efficient exportation of synthesized carbohydrates to sink tissues such as roots and fruits. Finally, the cambium tissue located in between the two supports secondary growth, thickening the stem and generating new secondary xylem and phloem cells (Agustí & Blázquez, 2020).

This comes together to explain the radially increasing trunk diameter of H. brasiliensis is the chronological series of cambial oscillations, such that each concentric ring band (xylem, phloem, and cambium) represents the temporal snapshot of hormonal (auxin) and metabolic (sucrose) wave dynamics. The pattern is functionally analogous to Liesegang rings, which arise through reaction-diffusion process involving the periodically generated, spaced precipitate bands. In H. brasiliensis, periodic surges in auxin gradients or sucrose act as biological reaction-diffusion fronts. These fronts propagate outwards through the cambial zone, and the oscillatory behaviour becomes spatially “frozen” as new xylem and phloem layers differentiate, producing the ring band repeatedly, explaining the ring formation over time, and parallels the temporal encoding mechanism seen in Liesegang systems (Nabika et al., 2019).

Growth Rate and Allometry Governing Hevea brasiliensis

Balance Between Vertical Growth and Stability

In the previous two essays on Hevea brasiliensis, there has been a focus on how the tree allocates its resources to optimize survival. In fact, allometric equations describing relationships between different growth characteristics offer insight into the distribution of energy across leaves, trunk, and roots by their growth characteristics. In addition to this, these equations show how structural stability is maintained through multiple stages of its life. 

To start, the height is a crucial component of H. brasiliensis, it must be tall enough to intercept enough sunlight, but it must also stop at a finite, optimal height (30-40 m). This balances the demand of a tall tree on a thick trunk that supplies nutrient and stability. Oyebade and Onyambo developed a general model describing the relationship between diameter at breast height (DBH) and height (H) to characterize the growth of Hevea brasiliensis across 20 different plots with varying tree densities (Table 3) (Oyebade & Oyambo, 2011). 

Table 3: Different diameter at breast height (typically 1.3 m for H. brasiliensis, being the estimated point of germination) with corresponding height for different sample sizes of H. brasiliensis in Nigeria (Oyebade & Oyambo, 2011). 

Altogether, different plots conveyed a universal pattern when DBH and H is plotted. As seen in Plot 13, and 10, with the lowest mean diameter at breast height of 29.82 cm and highest mean diameter at breast height of 58.49 cm respectively, both convey a similar, sigmoidal shape graph (Fig. 4) (Oyebade & Oyambo, 2011). 

Fig. 4. Relationship between height and diameter at breast height (dbh) among largest sample size (top), and smallest sample size (bottom) (Oyebade & Oyambo, 2011). 

This sigmoidal shape of graph is commonly found in nature, particularly with organism growth. The best-fit model for H. brasiliensis’ height and diameter at breast height is Monserud’s nonlinear equation (Eq. 11):

Where ht = height, bh = breast height, D = diameter, and β₀, β₁ are arbitrary parameters unique to different species of trees (Oyebade & Oyambo, 2011). In a similar study with hardwood species β₁ is concluded to always be a negative coefficient, and as D goes to infinity, a negative β₁ value results to zero (Colbert et al., 2002). 

From the mathematical model that depicts the relationship between diameter and height of H. brasiliensis, it reveals a pattern of growth constraint on height as diameter increases. This is because as height increases, it demands both a higher energy output to sustain its survival, as well as a more robust trunk structure to maintain its structure; but a sub-linear scaling pattern between the two that H. brasiliensis follows sets a restraint that prevents it from developing to height with an unsustainable structure and energy demand. 

Growth Characteristics with Relationship to Age and Environment

Currently, there exist numerous models that account for the different biomass accumulation of H. brasiliensis, with most of them estimating biomass using height and diameter at breast height. However, a recent study by Chen and colleges in 2020 points out a neglected variable, the effect of different age on the existing models of allometric equations.

From this, the observed partitioning of biomass is around 87.6% above ground and 12.4% below ground. Using the stand age of the tree (i.e., the average age of H. brasiliensis in a forest, in this case, forests in Hainan, China), a new equation was derived. By changing the independent variable from girth to stand age, Eq. 12 is used to approximate the total biomass as follows:

In this case, B_RG = biomass of rubber plantation, with x = stand age. This offers an alternative perspective by capturing the underlying patterns behind H. brasiliensis’ biomass accumulation through its changes with both age and environmental conditions. Using biomass as the dependent variable and stand age and environmental factors as the independent variables still produces a strong mathematical relationship. When biomass estimates come from stand age are combined with predictions of a random forest model that incorporates environmental variables such as canopy reflectance, soil properties, and other site conditions synthesize a resulting coefficient of determination ranges from 0.83 (Fig. 5) to 0.94 (as reported in study). These values are close to those obtained through traditional allometric models based on diameter at breast height and height (Chen et al., 2020). 

Fig. 5. (A) Scatter plot comparing the observed biomass (B_G) against an estimated biomass using a regression model, with stand age as independent variable. (B) A Random Forest (RF) model based on both stand age and Landsat and Sentinel-2 (LS2) - based variables - which focuses on both natural and human-inflicted changes measured using multi-spectral model such as forest health and natural disaster over time. (C) A combination of the two models conveys the most accurate model between biomass and environmental factors over the period of its life (Chen et al., 2020; Burns, 2025).

From the two models, H. brasiliensis shows a common growth rate with other trees that can be modelled by allometric equations that imposes an inherent, internal restraint that optimizes its size (increasing its light interception - as mentioned in the physics paper), but also refrains from an overbearing height that demands an unrealistic strength from structural integrity and energy supply. In addition to this, the temporal pattern almost demonstrates an internal rhythm. The combined Random Forest and regression model shows a steeper early-stage increase. This aligns with the tree’s biology: rapid height growth, fast canopy expansion, and quick lateral and taproot development in its early days. Meanwhile, an underestimation (plateau) in its latter phases suggest that much of its biomass remains the same with paused growth of its structures.

Periodic and Fractal Structure of the Root System

For many species, evolution and adaptation to challenges have refined their respective characteristics, ranging from resource acquisition to resistance to mechanical stresses, often creating multifunctional organs to maximize resource utilization. H. brasiliensis is no exception. Different mathematical models show us that the root system of the rubber tree, which provides nutrients, water, and mechanical anchorage to the tree, grows in a specific architecture that helps serve these multiple functions.

Periodic Branching of Lateral Roots

The essay covering the physics of H. brasiliensis showed that its lateral root system had evolved to form superposed ring patterns around the taproot (Fig. 6). By taking this shape, the root system provided additional anchorage for the tree, with upper layers of roots participating in intra-specific root grafting with nearby H. brasiliensis specimens, and deeper layers forming a girthier Root-Soil plate.

Fig. 6. Simplified representation of the H. brasiliensis root system. Taproots (T) grow vertically from the base of the tree along axis (1), lateral roots (L_n) grow from the taproot along axis (2). λ represents the distance between two rings of lateral roots (Adapted from Priyadarshan, 2017).

This pattern can be explained or modelled mathematically through the study of the relation between shoot and root growth. In Thaler and Pagès (1996), an experiment compared the growth of leaf area to the growth of taproots, lateral roots (secondary roots, tertiary roots) every 2 days among 10 individual trees. Results show that leaves grew in flushes, short periods of accelerated growth between periods of stagnation (Fig. 7).

Fig. 7. Shoot length (solid lines) and cumulative leaf area (dotted lines) versus time for 10 numbered individuals (Thaler & Pagès, 1996).

Though the taproot growth rate did not show any significant correlation with leaf growth, continuously staying above 0.4 cm/d, lateral root growth seemed to follow the flushes. As seen in Figure 8, secondary root growth decreased within the first 15 days after a small peak of initial growth, when shoot growth started to increase. Then, it increased to its maximum, around 30 cm d^-1, right when the first flush ended, and the process was repeated for the second flush. Tertiary roots followed a similar pattern, growth decreasing to 0 during flushes.

Fig. 8. Total growth rate of secondary roots (solid lines) and leaf area growth rate (dashed lines) versus time for 10 numbered individuals. Total secondary root growth rate is the sum of all increments of secondary roots caused by branching or elongation (Thaler & Pagès, 1996).

Decreased growth during flushes also meant that fewer secondary roots grow during these periods. This forms the alternating pattern of fully grown and branched laterals and poorly grown and branched laterals along the taproot. This can be modelled mathematically by the following equation or sequence:

Where x_n is the position of lateral root rings (fully grown and branched) along the vertical axis originating from the taproot base, n the index of the lateral root ring, Δt the time between flushes and R the taproot elongation rate.

This mathematical model has been hypothesized as either the result of a chemical stimulation kickstarted by factors synthesized in young leaves such as auxin (Atzmon et al., 1994), and / or as the result of the competition for assimilates between the shoot and the root of the H. brasiliensis (Glenn & Welker, 1993). It serves as an optimized baseline for the mechanical anchorage potential of roots, a design solution that uses the chemical cycles of the tree to respond to physical challenges through a mathematically regular geometric structure.

Fractal Structure of the Root System

However, another model used to represent root system reveals another potential mathematical design solution. A fractal is a mathematical object with self-similarities that can be observed at any scale. Interpreting the root system as one corresponds to looking at root branching patterns as a multitude of smaller patterns that have approximately the same shape as the overall system (Pages et al., 2000).

Examples of fractals in nature include cauliflowers, whose little heads are a quasi-exact smaller copies of the whole head. They are formed from the intersection of Fibonacci spirals of smaller heads, which are themselves spirals of reduced heads, and so 50 times, as seen in Figure 9(Havlin et al., 1995).

Fig. 9. Photograph of Romanesco Curd (Cauliflower)’s fractal heads (Azpeitia et al., 2021).

Fractal analyses are often conducted in a two-dimensional plane, to which H. brasiliensis’ ring-shaped lateral root layers are similar. They consist of placing the overall root structure on a square and dividing that square into smaller squares of length r and counting the number N of squares the roots touch and repeating the process for smaller squares. The branching properties of roots on planes can be evaluated using the following equation, as in Tatsumi et al. (1989):

Where K is a proportionality coefficient and D is the fractal dimension. D = 1 is a line, D = 2 is a plane.

The study conducted in Tatsumi et al. (1989), which studied multiple specific root systems, came to the conclusion that roots are self-similar objects with a fractal dimension between 1.48 and 1.58, with self-similarities ranging in size from 0.3 to 20 mm. Another key characteristic of fractals, in 2 dimensions as in 3, that is crucial to understand is their total surface area and volume (for 3D systems). When calculating the volume of a fractal object when the fractal dimension diverges to infinity, it can be found that it tends to 0, whereas surface area tends to infinity. A fractal in a 3D space cannot exceed a fractal dimension of 3, but the greater it is, the more surface it covers and the less volume it takes.

This indicates that roots systems having a fractal dimension above 1 maximizes the root-soil interface by reducing the volume the roots take up in the soil while increasing the surface area of contact between the roots and the soil. This is the second mathematical design solution: less volume, and a greater root-soil interface increases soil exploration and consequently, nutrient and water absorption. For H. Brasiliensis, who lives in tropical climates susceptible to droughts and heavy competition from both other beings, it is simply the result of evolution optimizing its chances of survival.

How the Mating System Leads to Genetic Diversity

The chromosome count of H. brasiliensis is 2n=36, signifying it is biologically a diploid species. The tree is believed to have ancient amphidiploid origin which have stabilized during evolution. It is a preferentially cross-pollinating, perennial species, with a long breeding cycle (Brito da Costa et al., 2000). This has resulted in the creation of a large genetic diversity within the species (Priyadarshan, 2017). The flowers of H. brasiliensis are pollinated by insects of the Ceratopogonidae family and thrips which are tiny insects that are responsible for cross-pollination over short distances (Furlani et al., 2005).

Mating System

The mating system establishes how genetic information is transmitted from one generation to the next and can partially determine the genetic structure of the population. In order to maximize the selection process, the long breeding process has resulted in several studies of genetic parameters (Brito da Costa et al., 2000). A study was conducted to estimate the outcrossing rate, which measures how often cross-pollination occurs compared to self-pollination, in H. brasiliensis. The test included 37 mother trees with around 10 trees per plot. The seeds that were collected were restricted to a 10 to 15 m radius from where the mother tree is located. This was done to make certain that the seeds that were being collected came from the mother tree. The models used in the study were the multilocus mixed mating model which doesn’t require the assumption of inbreeding equilibrium as well as the correlated mating model which allows for the separation of progeny fractions created by random matings of the biparental matings. In the analysis these models were performed using the Multilocus “t” and “r” (parameters) (MLTR) Mating System Program. The study was done in two steps which can be seen in Figure 10 (Furlani et al., 2005).

Fig. 10. Parameters defined and found throughout the analysis (Rosa Muldoon’s own artwork).

Due to the mother trees being not directly genotyped, the mother’s genetic makeup was inferred by looking at the genes of her offspring. This was achievable by using the MLTR software in which the genotypes of the offspring from the mother trees are inputted, then the software uses the method of most likely maternal genotype to compare the observed offspring genotype to the expected genotype under different selfing or outcrossing rates. In order to check the reliability, the standard error was calculated based on 500 bootstrap (computer simulation) resamplings. The standard error was estimated at 95% or can also be expressed as 1.96 SE, which demonstrates a 95% confidence the results were accurate (Furlani et al., 2005).

To test whether mating within the population was random, an analysis was done which estimated the genetic divergence between the allele frequencies of ovules and pollen pool using F_ST statistics. The chi-square (c²) test, shown with Equations 15 and 16, can be expressed in two parts, and is used to determine whether the observed values differ significantly from the expected values.

Where n is the number of gametes in the two groups (pollen and ovules), k is the number of alleles, s is the number of groups, F_ST is the measure of genetic differentiation between the groups, and df is degrees of freedom, in other words how many independent pieces of information are used to calculate the statistic (Furlani et al., 2005).

Next, the genetic relatedness among offspring in the same family also known was estimated to further understand the mating patterns. This reveals how much the offspring shares genes and whether they come from the same or different fathers. This was achieved using the correlation of relatedness among plants within families (r_xy), shown by Equation 17.

Where F_p is the fixation index, a measure of how much inbreeding there is in a population, in the parent population, s is the selfing rate therefore how often a tree mates with itself, t_m is the multilocus outcrossing rate, meaning how often the tree mates with other trees, r_s is the correlation of selfing and r_p is the correlation of paternity, how likely offspring have the same father. In an ideal, healthy, randomly mating diploid species, it can be assumed that there is no inbreeding, then the relationship between two siblings can be expressed more simply as seen in Equation 18.

q_xy measures the probability that two individuals identical by descent. The equation presents a simple way to quantify genetic relatedness between two individuals when there is no inbreeding. Equation 19 shows how the inbreeding (F) within each family can be calculated.

Where H_e is how much genetic variation is expected under random mating, also known as the expected heterozygosity, in the Hardy-Weinberg equilibrium, Equation 20.

n is the number of sampled trees, p_i is the frequency of allele i in a specific locus, and H_o is the observed heterozygosity which is expressed in Equation 21.

In this equation P_ii² represents the frequency of homozygous genotypes ii in a particular locus (Furlani et al., 2005).

The mating system of H. brasiliensis was analyzed using six enzyme loci. Of these six only three (6-Pgdh-1, 6-Pgdh-2, and Pgi-2) of them demonstrated significant genetic differences between pollen and ovule pools which can be seen in Table 4. This indicates a heterogeneous pollen population due to outcrossing rates or the geographic distance between the trees.

Table 4. The allelic frequencies, genetic divergence in pollen and ovules (F_ST) and chi-square test (c²) of H. brasiliensis. These values demonstrate how genetic diversity is distributed among and within families. The low F_ST values show a low genetic differentiation between pollen and ovule populations which suggest random mating (Furlani et al., 2005).

The fixation index for the mother tree (F_m=0.04) was found to suggest random mating and Hardy-Weinberg equilibrium among adults due to it not being significantly far from zero. The observed inbreeding level (0.142) was greater than what was expected in the study from the selfing (0.019). This indicates that the mating among related trees is what caused homozygosity rather than self-fertilization. The multilocus and unilocus outcrossing rates were both found to be high, t_m=0.962, t_s=0.948. This shows that H. brasiliensis mainly outcrossed but did allow for some self-crossing. The small outcrossing rate between individuals (t_m- t_s= 0.014) and low correlation of selfing (r_s=0.096) signify that the mating between relatives was limited. These results are summed up in Table 5. The correlation of paternity indicates that about 35% of the offspring within a family share the same father. The results show that H. brasiliensis undergoes primarily cross-pollination with infrequent self-crossing (Furlani et al., 2005).

Table 5. Estimate of parameters of inbreeding and mating system of H. brasiliensis (Furlani et al., 2005).

The overall outcrossing rate was 0.962, meaning the trees mainly cross-pollinate rather than self-pollinate, which can be observed in Figure 11. The mating system of H. brasiliensis can be interpreted in terms of being a design solution of the tree. Cross-pollination ensures and promotes genetic diversity within the species which is important for adapting to changing environments such as diseases or climate change. H. brasiliensis has also incorporated in its mating system a safety mechanism. The occasional inbreeding that does occur allows the tree to reproduce even when pollinators or other trees are limited. This allows for population stability in isolated or disturbed areas.

Fig. 11. Distribution of frequency of the individual outcrossing rate in a H. brasiliensis population (Furlani et al., 2005).

Latex Flow Rate

As discussed in the first two essays, the physical and chemical mechanisms responsible for latex exudation and cessation of flow involve turgor pressure and latex vessel plugging. When tapped, the turgor pressure of the laticiferous system near the cut decreases from around 10 atm to atmospheric pressure (An et al., 2014; Yeang, 2005). This decrease results in laticifer wall contraction, latex exudation, and bursting of lutoids (An et al., 2014; Wang et al., 2013). The loss of turgor pressure is responsible of the rapid decrease in latex flow rate, but complete flow cessation is due to coagulation of latex by released lutoid proteins (Ng et al., 2022; Wang et al., 2013; Yeang, 2005).

Latex Vessel Plugging

In 2005, Yeang proposed a mathematical model describing latex vessel plugging over time by measuring the intensity of plugging (IP). At a given point in time during latex flow, the IP is defined by the following formula, where the yield is determined by weight and the time is counted from the time when the first drop of latex has reached the collecting vessel (Yeang, 2005):

Therefore, re-tapping before cessation of the flow from the original tapping enables the extent of cumulative latex vessel plugging up to the time of re-tapping to be estimated as the intensity of plugging. It is important to note that, since re-tapping before stoppage of flow affects the subsequent flow characteristics, it is only possible to re-tap a tree once per recording day. Therefore, to obtain a representation of the plugging behavior for the entire duration of flow, re-taping must be done on different trees at different times. On each recording day, 8 rubber trees were subjected to one of 8 treatments that involved their being re-tapped half-spiral or not re-tapped. The scheme of experimental treatments in one treatment cycle, where all 8 treatments were applied in rotation over 8 recording days, is shown in Table 6 (Yeang, 2005).

Table 6. Treatment cycle of the 8 rubber trees with medium flow durations (Yeang, 2005).

As shown in Table 7, six other trees were tapped alternatively daily with a shortened tapping cut of one-eighth spiral to increase the rate of latex vessel plugging for shorter flow duration, and another set of six trees were tapped half-spiral on alternate days, while being stimulated with 2.5% chloroethylphosphonic acid, ethephon, painted on the tapping cut once every four days, to decrease rate of plugging for longer flow duration. In all the experiments, the total duration of latex flow was determined from the tree that was not re-tapped on each recording day. The average flow durations for the short, medium, and long flows were respectively 64, 126, and 217 min. The value of the medium flow is a mean of 24 readings, or 3 cycles, and the values for the short and long flows are the means of 12 readings, or 2 cycles.

Table 7. Treatment cycle of the two sets of 6 trees with either short or long flow durations (Yeang, 2005).

Figure 12 represents the combined results of IP across the short, medium, and long flow durations. In each case, IP, or the percent cumulative latex vessel plugging, is proportional to the square root of time:

where y is IP, x is the time from original tapping, and a, b are constants.

Fig. 12. Intensity of plugging (IP) in time. (A) Short flow with curve y_A1 = 12.41√x + 0.16. (B) Medium flow with curve y_B1 = 9.00√x - 3.20. (C) Long flow with curve y_C1 = 7.03√x - 2.85 . Each point represents the means of 24 readings in (A), and 12 readings in (B) and (C) (Yeang, 2005).

At x = 0, cumulative plugging y should be 0 if all latex vessel plugs have been removed upon tapping. Thus, a can be dispensed from equation (24):

At the time of flow cessation t, all the latex vessels are plugged, and the cumulative plugging on the tapping cut y is 100%. Thus, at t, 100 = b√t, and equation (24) becomes:

With equation (25), the cumulative plugging at any time during latex flow can be estimated. For example, at the mid-point of the total flow duration, x = t/2, y = 70.7, which means that 70.7% of the latex vessels at the tapping cut are plugged (Yeang, 2005).

Plugging Rate

By computing the first differential of equation (26), it is possible to obtain an equation for the rate of latex vessel plugging over the entire tapping cut:

Equation (26) demonstrates that the rate of plugging decreases as latex flow progresses, since less and less vessels are left to plug. To obtain the tendency of a yielding vessel to plug, the average rate of plugging of the remaining unplugged vessels can be calculated. Using equation (25), at time x, the percentage of latex vessels that have been plugged y is 100√(x/t) , so the percentage of latex vessels still yielding is 100 - 100√(x/t). Hence, the relative plugging rate y_r for x > 0, is given by the rate involving all the latex vessels severed at the tapping cut divided by the portion of latex vessels still yielding:

Figure 13 displays y_r for short, medium, and long flow durations. For each flow duration, y_r is relatively low in early and mid-flow, before increasing sharply towards the end of flow, especially just before flow cessation (Yeang, 2005). In all cases, the positive feedback loop regulation mechanism of the relative plugging rate during the critical event of wound plugging can be observed, where the plugging rate follows exponential growth after a plateau leading to an effective change of the status quo, from latex exudation up to flow cessation.

Fig. 13. Plugging rate of yielding latex vessels over time. Curves done with equation (27) for the short, medium, and long flow durations. Correction: highest value on the x-axis should read “250” (Yeang, 2005).

Turgor Pressure and Latex Exudation

Turgor pressure is instrumental in expelling latex from a tapping cut. The latex flow rate from the whole tapping cut declines to a slow trickle before stoppage, because more and more latex vessels are plugged. However, this does not imply that the flow from unplugged vessels also decreases (Yeang, 2005).

Immediately after re-tapping, any increase in flow rate is deemed to be due to the removal of plugs that have formed, and all latex vessels at the tapping cut would be actively yielding. Thus, by assuming random plugging of latex vessels and no change in turgor pressure in the 30 s before and after re-tapping, the proportion of latex vessels still yielding at x min would be:

By measuring the flow rate from the entire tapping cut at x min in the field, the relative flow rate from yielding latex vessels would then be proportional to:

By combining equations (28) and (29), it can be determined that the average flow rate from individual unplugged latex vessels would be proportional to the flow rate from the entire tapping cut after re-tapping. Figure 14 shows the flow rate after re-tapping in relation to the tapping panel turgor pressure (TP) for the three flow durations. The curves for flow rates after re-tapping also represent the flow rates from unplugged latex vessels, since the former is proportional to the latter. Immediately after tapping, the flow rate per yielding latex vessel dropped markedly, but a gradual recovery to 65-75 % of the initial rate followed. This means that the slow drip from the whole tapping cut near flow cessation does not reflect the flow from individual latex vessels. Moreover, the flow rate from the tapping cut after re-tapping, which is proportional to the flow rate from unplugged vessels, changes in tandem with TP. Therefore, TP is primarily responsible for the rate of latex exudation from a yielding latex vessel. In all cases, TP decreased up to 2-3 atm shortly after tapping, before recovering to 5-7 atm at flow cessation (Fig. 14). This further indicates that stoppage of latex flow cannot be attributed entirely to TP loss (Yeang, 2005).

Fig. 14. Panel turgor pressure (TP) and flow rate after re-taping in time. Short (A), medium (B), and long (C) flow durations. Each point represents the means of 24 readings in (A), and 12 readings in (B) and (C). Correction for (B): lowest values on the y-axis should read “0”, “2”, and “4” (Yeang, 2005).

Latex Flow Rate: A Function of Time

As shown in Figures 12 and 14, TP and IP are major determinants in regulating the latex flow rate. Based on equation (25), IP is proportional to √(x/t). Therefore, TP and IP data can be fitted into the following linear multiple regression model:

Where r is the latex flow rate, and a, b₁ and b₂ are constants. The parameter b₂ incorporates √(1/t) since t is constant for a particular tapping. Figure 15 shows that model (30) explained 99% of the variation in flow rate for all flow durations (Yeang, 2005).

Fig. 15. Measured and predicted flow rates. Short (A), medium (B), and long (C) flow durations. Left panel: Predicted and measured rates of latex flow at intervals after tapping. The multiple regressions for the flow rate predictions are r_A4 = 2.448 + 0.163(TP) - 0.455√x, r_B4 = 3.351 + 0.440(TP) - 0.538√x, and r_C4 = 5.365 + 0.158(TP) - 0.459√x. Right panel: Linear regressions of the measured flow rate with the predicted flow. TP in atm; x in min (Yeang, 2005).

Based on the results of Figure 15, it can be established that both cumulative latex plugging and turgor pressure are the parameters that regulate latex flow rate over the entire duration of latex flow. It is important to note that, since IP is itself a function of time, latex flow rate can be expressed as a function of TP and time, without having to invoke separately the considerations of fluid dynamics. As shown in Figures 12, 13, and 14, loss of TP is responsible for the rapid decrease in latex flow after tapping, while latex vessel plugging, which is more active in the late flow, leads to flow cessation (Yeang, 2005). This indicates that the mechanism behind the rate of latex exudation is a two-phase mechanism, where turgor pressure allows for the rapid release of latex to the site of injury, and lutoid proteins plug the latex vessels to prevent further metabolite loss and pathogen invasion (Gidrol et al., 1994; Wititsuwannakul et al., 2008; Yeang, 2005).

Conclusion

H. brasiliensis has no shortage of design solutions to the challenges its environment imposes, be they physical, chemical or mathematical. Reaction-diffusion phenomena in the tree cause the formation of a patterned vascular tissue, aiding in structural support as much as in the transport of metabolites throughout the tree. This also allows for the necessary formation of secondary xylem and phloem cells. To complement this, H. brasiliensis’ growth rate follows allometric equations, which allow it to experience exponential growth in its early years in relation to its diameter. However, this eventually plateaus at a certain height, effectively preventing itself from extending past its capacity to maintain itself mechanically or energetically. Growth also takes specific turns underground: the periodic growth of lateral roots in relation to leaf growth allows it to create a pattern of roots that further reinforces mechanical resistance; and the fractal organization of these same roots optimizes the root-soil interface, allowing for better nutrient and water absorption. Additionally, H. brasiliensis’ preference for cross-pollination enhances the genetic diversity of the species, enabling better adaptation to its environment. Finally, the integrity of the H. brasiliensis is protected by latex, whose flow rate depends on turgor pressure for efficient release of latex, and on wound plugging for cessation of flow by coagulation.