# A Reinforcement Learning approach to study climbing plant behaviour

Lucia Nasti, Giacomo Vecchiato, Patrick Heuret, Nick P Rowe, Michele Palladino, Pierangelo Marcati

► To cite this version:

Lucia Nasti, Giacomo Vecchiato, Patrick Heuret, Nick P Rowe, Michele Palladino, et al., A Reinforcement Learning approach to study climbing plant behaviour. Scientific Reports, 2024, 14 (1), pp.18222. 10.1038/s41598-024-62147-3 . hal-04685621

HAL Id: hal-04685621

https://hal.inrae.fr/hal-04685621v1

Submitted on 3 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# A Reinforcement Learning approach to study climbing plant behaviour

Lucia Nasti<fnref n="1" />✉, Giacomo Vecchiato<fnref n="1" />, Patrick Heuret², Nicholas P. Rowe², Michele Palladino<fnref n="1" />,³ & Pierangelo Marcati<fnref n="1" />

A plant's structure is the result of constant adaptation and evolution to the surrounding environment. From this perspective, our goal is to investigate the mass and radius distribution of a particular plant organ, namely the *searcher shoot*, by providing a Reinforcement Learning (RL) environment, that we call Searcher-Shoot, which considers the mechanics due to the mass of the shoot and leaves. We uphold the hypothesis that plants maximize their length, avoiding a maximal stress threshold. To do this, we explore whether the mass distribution along the stem is efficient, formulating a Markov Decision Process. By exploiting this strategy, we are able to mimic and thus study the plant's behavior, finding that shoots decrease their diameters smoothly, resulting in an efficient distribution of the mass. The strong accordance between our results and the experimental data allows us to remark on the strength of our approach in the analysis of biological systems traits.

Plants are living organisms coordinating a complex network of internal, e.g., nutrient concentration, and external signals, e.g., light and soil resources. As a result, plant growth is a delicate balance among different factors involving environmental and physiological conditions<fnref n="1" />.

Despite their sessile life, plants can move and react to external stimuli to search for nutrients, and avoid obstacles and dangerous conditions². In contrast to the animal kingdom, plants do not perform these movements only through “active” reversible actions, but also by expanding their organs³. Plants actually live in a complex environment with a limited amount of resources. These resources are shared between plants of the same species as well as plants of different species. In this context, the efficient use of such resources can be crucial for plant subsistence. For instance, in its growth process, a plant produces a limited amount of biomass. It can use such mass to elongate one of its stems (primary growth process) or to produce rigidity and stability of its organs (secondary growth). Efficient employment of biomass means finding the most convenient threshold between primary and secondary growth, considering that the longer the stem, the better the exploration, but the thicker, the more resistant it is to mechanical stress.

The concept of efficiency may be the key to understanding how plants interact with the environment and develop their organs. This gives a different perspective to plant modelling. Indeed, this point of view enriches a model by adding to the mathematical description of the biological system a possible interpretation. The possibility of mathematically studying the mechanics of a structure and understanding the extent of its physical limits has fascinated scientists since the time of Galileo (see in particular the work "Two new sciences"). Specifically on plants, there are studies on critical lengths (see for instance⁴⁻⁶), on the distribution of roots and branches⁷, or on how a root penetrates the soil⁸. The optimisation paradigm can be applied effectively in biological contexts⁹,¹⁰ and is what has stimulated this work.

In this work, we focus on a particular climbing plant species, known as the *Condylocarpon guianense* Desf. (see Fig. 1).

This plant species is a liana widely found in the flora of French Guiana, which twines around the branches and the trunks of neighboring plants in order to reach the canopy. Several studies on its structure, see¹¹⁻¹⁴ for instance, have revealed that in different growth stages, it changes the thickness and the nature of the layers that form its stem and consequently it changes its flexural rigidity. More specifically, the plant is more rigid during the self-supporting state, while it displays a less dense material and a thicker compliant cortex when attached to a support. Such a wide capability of *C. guianense* to adapt to the surrounding environment suggests that it is following a paradigm of efficiency, making it a suited subject for our study. In particular, we want to support the

<footnote n="1">¹Gran Sasso Science Institute, L'Aquila, Italy. ²AMAP, Univ Montpellier, CIRAD, CNRS, INRAe, IRD, Montpellier, France. ³DISIM, Department of Information Engineering, Computer Science and Mathematics, University of L'Aquila, Via Vetoio, 67100 L'Aquila, Italy. ✉email: lucia.nasti@gssi.it</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Searcher of Condylocarpon guianense in humid tropical forest canopy of French Guiana. Like many searchers the stem bears numerous expanded leaves and shows evidence of “adjustments” in stem direction along its length. In the lower part of the picture, additional searchers of the same species show the variety of searcher stem developments and orientations at different stages of development.</figcaption></figure>

hypothesis that the self-supporting organs of *C. guianense*, called searcher shoots, maximize their length avoid-
ing a maximal stress threshold. This idea is motivated by the fact that the searcher shoots are organs specialized
in finding a support to attach themselves. Hence, the longer they are, the better they explore the surrounding
environment and cross gaps. However, at the same time, they have to sustain their weight¹⁵.

To investigate this specific behavior and prove our hypothesis, we combine mechanical modelization and Reinforcement Learning (RL). Specifically, to prove that climbing plants optimize the mass distribution in their self-supporting stems, we developed a RL environment, which we called Searcher-Shoot, to study the radius along the climbing plant shoot. At the base of this environment, we considered two planar mechanical models: (MeLe) and (Me), which give us information about the development of the curvature and the mass along the shoot, considering the leaves or not, respectively.

The application of Artificial Intelligence (AI) to many biological problems is increasing rapidly in the analysis of plant morphology, growth, and development, or the understanding of their changing environment in con-junction with agriculture¹⁶. In particular, Machine Learning (ML) is playing a conspicuous role in developing predictive models in complex plant biological systems¹⁷ whenever possible the integration and the analysis of multidimensional omics data¹⁸. Moreover, as described in Ref. ¹⁹, with inadequate data, AI/ML applications perform poorly. When Supervised and Unsupervised methods are not able to generate a direct linear or non-linear mapping among the raw data, RL stands up for being a valid alternative method²⁰,²¹ and, as underlined in Ref. ²², it is emerging as a robust and reliable tool to face out real-world problems concerning biological systems. Example of this application are available in synthetic biology²³, metabolic engineering²⁴, chemical reaction network²⁵, and plant biology²⁶.

Our study delved into the complex and dynamic interactions that occur within plant growth behavior. Tra-
ditional numerical methods rely on deterministic models that may not accurately represent the nonlinear and
dynamic nature of plant behavior. In contrast, we utilized RL which is better suited to capture the nuances of
plant responses to various stimuli as it adapts to complex, changing environments.

RL's trial-and-error learning approach allows our model to explore and adapt to the dynamic behaviors of plants, enabling us to uncover optimal strategies that may not be evident in static, pre-defined models. Additionally, RL operates in a data-driven manner, which is particularly advantageous when studying plant dynamics. This allows us to extract patterns and behaviors that may be difficult to capture using theoretical, equation-based approaches.

To test our RL environment, we applied it to five samples of *C. guianense*, using experimental data provided by Refs.¹⁵,²⁷. We found that the optimal policy was able to reproduce the decreasing behavior that characterizes the radius of the sample in consideration. This outcome suggests that, at least for *C. guianense*, the mass distribution along the searcher shoot is distributed optimally to maximize the length while avoiding the curvature reaching a certain breakdown threshold.

Results

We develop the Searcher-Shoot environment in Python. Precisely, we employ the OpenAI Gym²⁸ and Stable-Baselines3 (SB3)²⁹ libraries, two open-source frameworks implementing several commonly used model-free deep RL algorithms. In particular, from SB3, we import the PPO algorithm³⁰. For the mathematical modelling, we develop two models for the searcher shoot: (1) a model with the mechanics, but without the leaves (Me) and (2) a model with both the mechanics and the leaves (MeLe). We perform all the simulations with the discount factor γ = 0.99 (see Supplementary Material) and we train the models setting the number of episodes to 1 million.

In Fig. 2, we plot the results of our simulations and we show the radius (Fig. 2a) and the mass distributions (Fig. 2b). We find that the radius decreases at each step, i.e. the agent chooses the actions leading to a smaller radius.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Radius and fresh mass distribution in the (Me) model. We show the radius and, consequently, the mass distribution obtained in the (Me) model after training. Specifically, in Fig. 2a, we compare the experimental radius (red) to the one obtained by our simulation (blue). Computing the relative error, we find that the discrepancy between the two radii is at most 8.55% in segment [0, 0.868]. The value 0.868 represents the total length in meters of the shoot sample. In Fig. 2b, we show the simulated mass. The mass decreases smoothly and is approximately zero at the tip.</figcaption></figure>

We compare our simulated scenario with samples S1–S5<sup>15,27</sup>. As we can notice in Fig. 2a, comparing the obtained radius distribution (in blue) with the experimental radius of sample S2 (in red), the relative error is 8.55% in the segment [0, 0.868], which represents the length of the sample S2. In the Supplementary Material, we include the simulation results and the comparison with all the other four samples.

Successively, we consider the (MeLe) model, where, in addition to the mechanics features, we model also the mass of the leaves along the shoot. In Fig. 3, we present the results of our simulations. Both the radius (Fig. 3a) and the mass distributions (Fig. 3b) are consistent with the results of the preceding model. Again, we compare the experimental radius (in red) of sample S2 with the one we derive using the (MeLe) model (in blue), and, by computing the relative error, we find that the discrepancy between the two radii is 10.28% in segment [0, 0.868].

## Study on model sensitivity w.r.t leaves configuration

To evaluate the model's adaptability to variations in leaf configuration, we conducted supplementary simulations introducing systematic changes in the length of the internode. In contrast to the original model, where the internode length is fixed at 0.13m, we explored a spectrum of lengths randomly generated within the range of [0.01, 0.5] (the length is generated at the beginning of a simulation and it doesn't change within that simulation). The central aim of this investigation was to unveil the model's sensitivity to leaf positioning, yielding nuanced insights. Noteworthy findings surfaced when the leaf cluster was positioned near the base of the shoot, revealing an average relative error of 7.43% across 100 simulation rounds. In contrast, scenarios with internode lengths

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Radius and mass distribution in the (MeLe) model. We show the radius and, consequently, the mass distribution obtained in the (MeLe) model after the training. Specifically, in Fig. 3a, we compare the experimental radius (red) to the one obtained by our simulation (blue). Computing the relative error, we find that the discrepancy between the two radii is at most 10.28% in the segment [0, 0.868]. The value 0.868 represents the total length in meters of the shoot sample. In Fig. 3b, we show the simulated mass. The mass decreases smoothly and is approximately zero at the tip.</figcaption></figure>
surpassing 0.2m displayed a significantly higher average relative error, reaching 28% over the same number of simulations.

In Fig. 4, we show an example of the worst (Fig. 4a) and the best leaves configuration (Fig. 4b). We show on the y-axis the radius distribution, on the x-axis the length of the plant's shoot. The vertical lines in light blue represent the leaves' position. In plot 4b, the internode is 0.0548m long and the relative error in 7.29%. Instead, in plot 4a, the length of the internode is 0.2103m.

Furthermore, in the example of the worst-case scenario (4a), it becomes apparent that the mass of the leaves does not diminish to zero (see Fig. 5), indicating an interruption in the plant's growth due to the violation of the curvature threshold. This behavior can be explained by considering the influence of the weight of the leaves on the slender shoots of the plant, which causes the violation.

## Model insights

To demonstrate the effectiveness of the model, Fig. 6 depicts the median and variance calculated from 100 simulations, alongside the experimental radii corresponding to sample S1. Additional sample plots can be found in the Supplementary Material for comparison.

## Discussion

The relative error between the measured radius (data provided by¹⁵) and the simulated samples is less than 15% for all the samples and in both cases with or without leaves (see Tables 1 and 2). Such a small relative error suggests that the optimal policy successfully reproduces the radial profile of the samples. Going into more detail, in the (Me) case with generic $c_2$ and $\psi_0$ (respectively, the constant for the curvature threshold and the ratio between curvature and the fourth power of the radius at the base of the shoot), the error is between 12.9 and 16.8% (see Table 1), with the exception of sample S1. This error decreases between 7.8 and 9.8% when the constants are sample-specific. In particular, we observe that the coefficient of variation of $c_2$ is about 35% (see Table 3), while the coefficient of variation of $\psi_0$ is much greater since it is approximately 74%. This might imply that the environmental conditions have a relevant impact on the initial curvature of the shoot, while the stress threshold might be characteristic of the plant species. The difference between the length of the simulation and

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Example of different leaves configuration. We show the length and radius of the plant with different leaves' configurations (vertical lines in light blue). As we can notice, in the best scenario the leaves are close to the base of the shoot.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Mass distribution. We show the mass distribution associated to the example of the worst-case scenario of plot Fig. 4a. Here, on x-axis, we show the length of the shoot, on y-axis we show the mass of the plant. As we can notice, the mass does not go to zero, meaning that the plant's growth is interrupted by the violation of the curvature threshold.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Median and variance computed over 100 simulations. We show the comparison with the samples S1.</figcaption></figure>

<table><caption>**Table 1.** The table reports the Relative Error with and without estimated parameters (EP), without considering the mass of the leaves. To show the strength of our results, we compute the Relative Errors by comparing our results with the experimental radii of the samples$^{15,27}$. To understand the influence of two crucial parameters, $c_2$ and $\psi_0$, which are the curvature threshold and the initial curvature, respectively, we simulate our (Me) model by using their estimated values and an average value.</caption><thead><tr><td>Sample</td><td>Relative error (%) without EP</td><td>Relative error (%) with EP</td></tr></thead><tbody><tr><td>S1</td><td>8.6</td><td>13.6</td></tr><tr><td>S2</td><td>-</td><td>8.5</td></tr><tr><td>S3</td><td>16.8</td><td>8.2</td></tr><tr><td>S4</td><td>16.2</td><td>7.8</td></tr><tr><td>S5</td><td>12.9</td><td>9.8</td></tr></tbody></table>

<table><caption>**Table 2.** The table reports the Relative Error with estimated parameters, by considering the mass of the leaves. To show the strength of our results, we compute the Relative Errors by comparing our results ((MeLe) model) with the experimental radii of the samples$^{15,27}$.</caption><thead><tr><td>Sample</td><td>Relative error (%)</td></tr></thead><tbody><tr><td>S1</td><td>7.4</td></tr><tr><td>S2</td><td>10.2</td></tr><tr><td>S3</td><td>11.4</td></tr><tr><td>S4</td><td>14.1</td></tr><tr><td>S5</td><td>11.8</td></tr></tbody></table>

<table><caption>**Table 3.** The table reports the values of $\psi_0$, $c_2$, with average and coefficient of variation. The values of $\psi_0$ and $c_2$ are estimated utilizing the methods described in$^{33}$.</caption><thead><tr><td>Sample</td><td>ψ<sub>0</sub></td><td>c<sub>2</sub></td></tr></thead><tbody><tr><td>S1</td><td>-5.9e-12</td><td>8.7e-4</td></tr><tr><td>S2</td><td>-3.7e-12</td><td>3.6e-4</td></tr><tr><td>S3</td><td>-7.8e-12</td><td>4.8e-4</td></tr><tr><td>S4</td><td>-2.3e-11</td><td>6.1e-4</td></tr><tr><td>S5</td><td>-1.1e-11</td><td>4.9e-4</td></tr><tr><td>Average</td><td>-1e-11</td><td>5.6e-4</td></tr><tr><td>Coeff. of variation</td><td>73.61%</td><td>34.4%</td></tr></tbody></table>
the experimental one, displayed in Fig. 2a, can be explained by the experimental error in the mass and density measurements, or by some assumptions on the model, such as a constant Young modulus all along the shoot. Regarding the simulations with the leaves, i.e. the (MeLe) model, the relative error gets worse if compared to the (Me) model (see Table 2). This drop in the accuracy of the model can be due to the fact that the agent has no control over the distribution of the leaves. Moreover, the constants $\psi_0$ and $c_2$ are estimated on a model based only on the optimisation of the mass of the main stem. Hence, effective implementation of the leaves would require a more accurate model; nevertheless, even in our approximation, we get an error between 7.4 and 14.1%. This result suggests that a functional advantage of the plant during the self-supporting stage is the optimization of the main stem's mass. This is in line with the application of the behavioural ecology theoretical framework for plants (see³¹ for a survey on the behavioural ecology of climbing plants). In other words, according to this theory, plants have the capability to place their stems and other organs in accordance with optimal economic models. For instance, depending on the external supports in the surrounding environment, the risk of herbivory and the energetic stress, a shoot may delay leaf expansion or have short internodes. In our case, the plant has a limited amount of mass and can develop limited internal stress. A longer stem means, on one hand, a greater exploration capability, but on the other hand, a greater risk, since the structure is more fragile. Our RL environment gives a quantitative answer to the trade-off that the plant has to face, in the specific case of the C. guianense in its natural habitat. However, the generality of the equations at the base makes such an environment suitable for application to other plant species. Safety factors in terms of mechanical stability of stems and other plant organs is an informative way of relating size, length mass and stiffness to mechanical stability and height. The approach was not explored in this study primarily because of the longitudinal complexity and changing tissue patterning and stiffness (E) along the searcher stems. Our observations indicate that many searchers probably function under relatively low safety factors and rely on high stiffness values over developing wide diameters. Field observations even suggest that exceeding critical buckling length after crossing maximal reach gaps might be an advantage for contacting potential host supports just after a truly self-supporting phase¹⁴. In general, the use of a RL method to study a biological phenomenon can have some disadvantages, since it depends heavily on the concept of state and value function³². However, when properly formulated, it is a valuable instrument to predict and understand the behaviour of the phenomenon in consideration.

## Conclusion

In this work, we developed a RL environment called Searcher-Shoot to optimize the mass distribution along the main stem of a climbing plant during its self-supporting phase. At the base of the environment, there is a mechanical model that treats the plant's stem as an elastic rod. We applied the Searcher-Shoot environment to the experimental data using two variants of the mechanical model: with or without leaves. The radii that we obtained from the numerical simulation are in accordance with the data, supporting the hypothesis that climbing plants apply the strategy of growing the stem as long as possible. Building an RL environment is a compelling strategy to understand complex systems biology behaviors, especially in the absence of crucial data required for the application of ML. The RL approach we employed in this work has yielded promising results, making us even more hopeful for future developments. In line with this work, we plan to build a more sophisticated model that can consider the curvature's time development in addition to the optimal mass distribution. Indeed, plant movement is partially determined by the response to external signals. This response can be optimized to maximize (or minimize) a reward based on the signals themselves.

## Methods

### Notations for the mechanics of a planar elastic rod

Let $e_1, e_2, e_3$ a basis of orthonormal vectors for $\mathbb{R}^3$. We assume that the searcher shoot behaves like an inextensible and unshearable elastic rod³⁴ confined in the plane spanned by the vectors $e_1$ and $e_2$. The centerline of this rod lies on the curve $\Gamma \subset \text{span}\{e_1, e_2\}$. We parametrize the curve $\Gamma$ with its arc length $s$, so that if the length of the curve is $L$, we have $s \in [0, L]$. We denote with $\Gamma(s)$ the position in the plane of the point on $\Gamma$ whose arc length is $s$. With this parametrisation, $\partial_s \Gamma(s)$ represents the normal tangent vector to $\Gamma$ at the point $\Gamma(s)$. We denote with $\theta(s)$ the angle between $\partial_s \Gamma(s)$ and the vector $e_2$. Consequently, $\theta'(s)$ is the curvature of $\Gamma$ at the point $\Gamma(s)$.

We refer to $\Gamma$ as the current configuration, which corresponds to the actual shape of the rod when subject to external physical forces. To study the effects of gravity acting on that rod, we need to consider the intrinsic configuration, denoted as $\hat{\Gamma}$, which corresponds to the geometric curve assumed by the rod when there are no external forces acting on it. Since the rod is inextensible, the arc length parameter of $\hat{\Gamma}$ and $\Gamma$ is the same $s \in [0, L]$.

The difference between the curvature $\partial_s \hat{\theta}(s)$ of the intrinsic configuration and the curvature $\partial_s \theta(s)$ of the current configuration is proportional to the resultant moment (of force) $m(s)$ acting at the point $\Gamma(s)$ (and directed along $e_3$) through the Euler-Bernoulli formula:

$$m(s) = -E(s)I(s)(\partial_s \theta(s) - \partial_s \hat{\theta}(s)). \quad (1)$$

This relation holds because we are considering unshearable rods. With $E$ we are denoting the Young's modulus, which expresses the stiffness of the material, and with $I$ the second moment of area of the cross-section (with respect to $e_3$). We assume that the rod is in elastic equilibrium, that is, all the internal forces and moments are in balance with the external forces and moments. In this framework, considering equation (1) and the gravity force as the only external force acting on the rod, we can write the following differential equation (see for instance³⁵)

$$\partial_s(EI(\partial_s \hat{\theta} - \partial_s \theta))(s) = \sin(\theta(s))g \int_s^L \rho_3(s')A(s')ds'. \quad (2)$$
In this equation, g represents the gravity acceleration constant, $\rho_3(s)$ the volume density of the shoot/rod at the point $\Gamma(s)$ and $A(s)$ the cross-section area at that point.

We are now interested in computing the internal bending stresses. We know that the internal moment $m$ is generated by the deflection from the intrinsic configuration $\Gamma$. Indeed, this deflection generates an internal pressure called stress, that we denote with $\sigma$, and a deformation $\varepsilon$ of each element of the rod, called strain. Stresses and strains vary according to the position $\Gamma(s)$ on the rod, and depend also on the position on the cross-section. Since the rod is ushearable, the cross-section is always orthogonal to the tangent vector $\partial_s \Gamma$. To describe the position of a generic point on the cross-section at $\Gamma(s)$, we name

$$ \beta(s) = \frac{\partial_s^2 \Gamma(s)}{|\partial_s^2 \Gamma(s)|} \in \operatorname{span}\{e_1, e_2\} $$

the normal vector, which is orthogonal to $\partial_s \Gamma(s)$, and we consider the binormal vector

$$ \tau(s) = \frac{\partial_s \Gamma \times \beta(s)}{|\partial_s \Gamma(s) \times \beta(s)|}, $$

Then, the cross-section at the point $\Gamma(s)$ is a subset $C(s)$ of the plane $\operatorname{span}\{\beta(s), \tau(s)\}$ with the origin on the centerline. We define

$$ C(s, z) = \{w \in \mathbb{R} : \Gamma(s) + z\beta(s) + w\tau \in C(s)\}. $$

In this framework, the maximal bending stress at $\Gamma(s)$ results to be$^{36}$

$$ \begin{aligned} \sigma_m(s) &= \max\{|w| : C(s, w) \neq \emptyset\} \cdot \frac{m(s)}{I(s)} \\ &= \max\{|w| : C(s, w) \neq \emptyset\} \cdot E(s)|\partial_s \theta(s) - \partial_s \hat{\theta}(s)|. \end{aligned} \quad (3) $$

## Formulation of the models

We assume that the searcher shoot has a circular cross-section with radius $r$ and that the Young's modulus $E$ is constant all along the shoot. So, we have

$$ \begin{aligned} A(s) &= \pi r^2(s) \\ I(s) &= \frac{\pi}{4} r^4(s) \\ r(s) &= \max\{|w| : C(s, w) \neq \emptyset\} \end{aligned} $$

and we name

$$ u(s) = r^2(s), \quad \psi(s) = u^2(s)(\partial_s \hat{\theta}(s) - \partial_s \theta(s)), \quad \mu(s) = \int_s^L \pi \rho_3(s') u(s') ds' $$

So, equation (2) can be rewritten as

$$ \left\{ \begin{aligned} \partial_s \psi(s) &= c_1 \sin(\theta(s)) \mu(s) \\ \partial_s \theta(s) &= \partial_s \hat{\theta}(s) - \frac{\psi(s)}{u^2(s)} \\ \partial_s \mu &= -\pi \rho_3(s) u(s) \end{aligned} \right. \quad (4) $$

with

$$ c_1 = \frac{4g}{\pi E}. $$

These equations hold for a.e. $s \in [0, L]$. At the boundary of this domain, we assume (1) to know the angle at the base of the shoot $\theta(0) = \theta_0$; (2) that at the tip of the shoot, there is not any external weight so that the intrinsic curvature equals the current curvature. Using the functions defined above, this means $\psi(L) = 0$ and $\mu(L) = 0$; (3) the mass $M$ of the whole shoot is known, so that $\mu(0) = M$.

The effectiveness of the self-sustaining behavior of the shoot can be quantitatively evaluated considering a threshold for the maximal internal stresses. In other words, we would like to find a distribution of the mass such that the maximal stress at each point of the shoot $\sigma_m(s)$ does not cross some fixed value $\bar{\sigma}$. Employing equation (Eq. 3), this means that any solution of system (Eq. 4) which satisfies the boundary conditions above discussed, has also to satisfy the condition

$$ |\psi(s)| \le c_2 u^{3/2}(s) \quad \text{for every } s \in [0, L] \quad (5) $$

with

$$ c_2 = \frac{\bar{\sigma}}{E}. $$

We group all these considerations into the following problems.
**Problem of shoot growth with Mechanics** (Me). We want to find the maximal length L of the shoot for which there exists at each point a radius r(s) such that the following system is satisfied

$$ 
\left\{
\begin{array}{l}
\partial_s \psi(s) = c_1 \sin(\theta(s)) \mu(s) \\
\partial_s \theta(s) = \partial_s \hat{\theta}(s) - \frac{\psi(s)}{u^2(s)} \\
\partial_s \mu = -\pi \rho_3(s) u(s) \\
\psi(L) = 0 \\
\theta(0) = \theta_0 \\
\mu(0) = M, \mu(L) = 0 \\
|\psi(s)| \le c_2 u^{3/2}(s)
\end{array}
\right.
\qquad (\text{Me})
$$

**Problem of shoot growth with Mechanics and Leaves** (MeLe). In problem (Me) we consider just the weight of the main stem. However, in most of climbing plant species, a relevant part of the total biomass is due to the leaves. We assume that the leaves are not uniformly distributed along the shoot. On the contrary, we assume that they are located at intervals equally spaced. Moreover, we assume that the mass $m_{lm}$ of a single leaf at the point $\Gamma(s)$ depends just on the shoot radius $r(s)$ at $\Gamma(s)$ (so $m_{lm} = m_{lm}(r(s))$). Let $d_{lm}$ the distance between two leaves locations and $n_{lm}$ the (fixed) number of leaves at each location. Then, we name

$$s_i = i \times d_{lm} \text{ for } i = 1, \dots, N_{lm},$$

where $N_{lm}$ is the total number of leaves locations. Therefore, $s_i$ denotes the arc length of the i-th leaves location. Now, we want to compute how the weight of the leaves affects the shoot at the point $\Gamma(s)$. To achieve this, we subtract from the total leaves mass $M_{lm}$ the mass of the leaves in the shoot portion between the base and $\Gamma(s)$. We define

$$q_{lm}(s) = \left\lfloor \frac{s}{d_{lm}} \right\rfloor$$

and

$$\text{RES}_{lm}(s) = M_{lm} - n_{lm} \sum_{i=0}^{q_{lm}(s)} m_{lm}(r(s_i)),$$

where the operation $\lfloor x \rfloor$ is the greatest integer lower or equal than x. Then, to take the leaves into account, we compute the gravity force acting at the point $\Gamma(s)$:

$$-g \left[ \int_{s}^{L} \rho_{3}(s') A(s') ds' + \text{RES}_{lm}(s) \right] e_{2}.$$

This leads to another problem of length maximization. Like for the (Me) case, we want to find the maximal length L of the shoot for which there exists at each point a radius r(s) such that

$$ 
\left\{
\begin{array}{l}
\partial_s \psi(s) = c_1 \sin(\theta(s))(\mu(s) + \text{RES}_{lm}(s)) \\
\partial_s \theta(s) = \partial_s \hat{\theta}(s) - \frac{\psi(s)}{u^2(s)} \\
\partial_s \mu = -\pi \rho_3(s) u(s) \\
\psi(L) = 0 \\
\theta(0) = \theta_0 \\
\mu(0) = M, \mu(L) = 0 \\
|\psi(s)| \le c_2 u^{3/2}(s)
\end{array}
\right.
\qquad (\text{MeLe})
$$

## Reinforcement Learning framework

We define the Searcher-Shoot environment as a Markov decision process problem (see Supplementary Material). Here, the agent is the liana's searchers shoot, a structure that has presumably been selected during evolution to have the greatest reach to colonize its environment by adjusting on its mechanical properties and the taper of the shoot's diameter. Specifically, the agent learns how to complete the task (i.e., the mass distribution) in the highest number of steps, choosing radius values does not generate internal stresses over a fixed threshold.

The fundamental elements of the framework are:

* **State and observations.** At each step, the agent, in the current state, can observe the mass, the radius, and the curvature before the next move. The choice of such a state is due to algorithmic and biological reasons. Indeed, in our algorithm, the agent controls the distribution of the mass by choosing the radius at each step. To this aim, the agent is able to observe the (remaining) mass and the radius. Regarding the curvature, it is a stop criterion and an index of efficiency. Moreover, plants are able to sense their own curvature through specialized fiber cells³⁷. Plants are also able to sense their own inclination³⁸, but we prefer to keep the number of observable variables as low as possible in order to reduce the computational effort of the algorithm.

* **Actions.** In this framework, the action space is discrete. In the (Me) and (MeLe) models, at each time step, the agent can select one action among eleven options: it can leave the radius value unchanged or it can increase (or decrease) the radius of a certain quantity (±1 × 10⁻⁵, ±2 × 10⁻⁵, ±3 × 10⁻⁵, ±4 × 10⁻⁵, ±5 × 10⁻⁵).
Of course, this selection will influence the mass distribution: intuitively, the larger the radius value, the larger the mass allocated in the next step.

* **Reward.** Every time the agent moves to the next step, it will receive positive feedback equal to +1. Whether the move ends with the total mass equivalent to or less than 0, or the condition on the curvature (Eq. 5) is violated, the reward is 0 and the algorithm stops.

* **Episode and Reset.** The episode does not have a fixed term. Instead, it ends whether the mass becomes zero or negative or the picked radius causes the curvature to violate condition (Eq. 5), as displayed in Algorithm 1 and 2. Then, we set the system parameters and the observation space to their initial values.

```
Require: M > 0
r > 0
(ψ, θ, μ) ← (ψ₀, π/2, M)
while μ ≥ 0 and ψ ≤ c₂ × r³ do
    r ← r or rₛ or rₛ
    ρ ← cvd + bvd · r + avd · r²
    (ψ, θ, μ) ← Solve System (Me) with u = r², ρ₃ ≡ ρ,
                    ψ(0) = ψ, θ(0) = θ, μ(0) = μ
    reward = 1
end while
reward = 0
```

**Algorithm 1.** Algorithm of Shoot growth with mechanics.

```
Require: M > 0, Mlm ≥ 0, dlm > 0, dlm > h > 0
r > 0
dbase ← 0
(ψ, θ, μ) ← (ψ₀, π/2, M)
RES ← Mlm
while μ ≥ 0 and ψ ≤ c₂ × r³ do
    r ← r or rₛ or rₛ
    ρ ← cvd + bvd · r + avd · r²
    if dbase % dlm < h and RES > 0 then
        RES ← RES - nlm · (clm + blm · r + alm · r²)
    end if
    (ψ, θ, μ) ← Solve System (MeLe) with u = r², ρ₃ ≡ ρ, RESlm = RES
    ψ(0) = ψ, θ(0) = θ, μ(0) = μ
dbase ← dbase + h
reward = 1
end while
reward = 0
```

**Algorithm 2.** Algorithm of Shoot growth with mechanics and leaves.

**Models implementation and parameters**

To begin with, we implement the System of equations (Me) in the Searcher-Shoot environment. In such a system, we consider stress and strain as factors responsible for its shaping, together with gravity, which acts as an external force on the plant's structure, affecting its curvature. Moreover, the material density is a function of the radius r, defined as follows:

$$ \rho_3 = c_{vd} + b_{vd} \cdot r + a_{vd} \cdot r^2. \qquad (6) $$

We use Algorithm 1 to clarify the approach we implemented. Starting from an initial configuration of curvature, angle and mass $(\psi_0, -\pi/2, M)$, and given an initial radius $R_0$, at each step of the algorithm the agent chooses how much to increase or decrease the current radius. Consequently, the total mass decreases and curvature and angle change accordingly to System (Me). This process is repeated until the mass vanishes or the curvature constraint is violated. The values of the constants $c_2$ and $\psi_0$ are estimated in two ways, leading to two groups of simulations. In the first group, $c_2$ and $\psi_0$ are the same for all the samples, while in the second group, they are estimated specifically for each sample by utilizing the method described in³³.

Successively, we add the leaves' mass contribution, which affects the plant's weight remarkably. We implement the System of equations (MeLe), where we can notice that the effects of the leaves are visible on the curvature and, then, in the formulation of the equation of $\psi$. As for the material density, the mass of a single leaf depends on the radius $r$ according to the following relation:

$$ m_{ml} = c_{lm} + b_{lm} \cdot r + a_{lm} \cdot r^2. \qquad (7) $$

In Algorithm 2 we clarify how we implement our model in the RL context. The Algorithm iterations are similar to Algorithm 1, and in addition, we consider the mass of the leaves and check whether at the position of the agent there is a group of leaves or not. In Table 4, we include all the parameters used in the models.
<table><caption>**Table 4.** Parameters of the models. In the case of the parameters resulting from the fitting procedure, we report the source of the data on which functions (Eqs. 6, 7) are fitted.</caption><thead><tr><td>Parameter</td><td>Description</td><td>Source</td></tr></thead><tbody><tr><td>g</td><td>Gravity acceleration constant</td><td>15,27</td></tr><tr><td>E</td><td>Young's modulus</td><td>15,27</td></tr><tr><td>M</td><td>Main stem freshmass</td><td>15,27</td></tr><tr><td>c<sub>2</sub></td><td>Curvature threshold</td><td>33</td></tr><tr><td>ψ<sub>0</sub></td><td>Initial curvature</td><td>33</td></tr><tr><td>R<sub>0</sub></td><td>Initial radius</td><td>15,27</td></tr><tr><td>(a<sub>vd</sub>, b<sub>vd</sub>, c<sub>vd</sub>)</td><td>Parameters for volume density fitting</td><td>15,27</td></tr><tr><td>(a<sub>lm</sub>, b<sub>lm</sub>, c<sub>lm</sub>)</td><td>Parameters for leaves mass fitting</td><td>15,27</td></tr></tbody></table>

## Statement on research methodology

The research methodology employed in this study adheres strictly to ethical and legal considerations associated with experimental research and field studies on plants. We affirm our compliance with relevant institutional, national, and international guidelines and legislation governing such research endeavors. Additionally, we emphasize our commitment to following the IUCN Policy Statement on Research Involving Species at Risk of Extinction and the Convention on the Trade in Endangered Species of Wild Fauna and Flora. This ensures the responsible and ethical conduct of our research with due consideration for the well-being of the studied plant species and the broader ecosystem.

## Data availability

All data generated or analyzed during this study are included in this published article and in its Supplementary Material.

Received: 13 November 2023; Accepted: 14 May 2024

Published online: 06 August 2024

## References

1. Lawlor, G. et al. *Plant Physiology* (Springer Science & Business Media, 2012).

2. Forterre, Y. Slow, fast and furious: Understanding the physics of plant movements. *J. Exp. Bot.* **64**(15), 4745–4760 (2013).

3. Rivière, M., Derr, J. & Douady, S. Motions of leaves and stems, from growth to potential use. *Phys. Biol.* **14**(5), 051001 (2017).

4. Greenhill, G. Determination of the greatest height consistent with stability that a vertical pole or mast can be made, and the greatest height to which a tree of given proportions can grow. in *Proceedings of the Cambridge Philosophical Society*. Vol. 4. p. 65. (1881)

5. McMahon, T. Size and shape in biology. *Science* **179**(4079), 1201–1204 (1973).

6. Wei, Z., Mandre, S. & Mahadevan, L. The branch with the furthest reach. *Europhys. Lett.* **97**(1), 14005 (2012).

7. Conn, A. et al. Network trade-offs and homeostasis in Arabidopsis shoot architectures. *PLoS Comput. Biol.* **15**(9), e1007325 (2019).

8. Tedone, F. et al. Optimal control of plant root tip dynamics in soil. *Bioinspiration Biomimet.* **15**(5), 056006 (2020).

9. Alexander, R. M. *Optima for Animals* (Princeton University Press, 1996).

10. Lenhart, S. & Workman, J. T. *Optimal Control Applied to Biological Models* (CRC Press, 2007).

11. Rowe, N. P. & Speck, T. Biomechanical characteristics of the ontogeny and growth habit of the tropical liana Condylocarpon guianense (Apocynaceae). *Int. J. Plant Sci.* **157**(4), 406–417 (1996).

12. Rowe, N., Isnard, S. & Speck, T. Diversity of mechanical architectures in climbing plants: An evolutionary perspective. *J. Plant Growth Regulat.* **23**(2), 108–128 (2004).

13. Rowe, N. P. & Speck, T. Stem biomechanics, strength of attachment, and developmental plasticity of vines and lianas. *Ecol. lianas* **32**–341 (2015)

14. Soffiatti, P. et al. Trellis-forming stems of a tropical liana Condy locarpon guianense (Apocynaceae): A plant-made safety net constructed by simple “start-stop” development. *Front. Plant Sci.* **13**, 1016195 (2022).

15. Hattermann, Tom et al. Mind the gap: reach and mechanical diversity of searcher shoots in climbing plants. *Front. Forests Glob. Change* **5**, 836247 (2022).

16. Soltis, P. S. et al. Plants meet machines: Prospects in machine learning for plant biology. in *Applications in Plant Sciences* Vol. 8.6 (2020).

17. Hesami, M. et al. Machine learning: Its challenges and opportunities in plant system biology. *Appl. Microbiol. Biotechnol.* **106**(9–10), 3507–3530 (2022).

18. Ma, Anjun et al. Integrative methods and practical challenges for single cell multi-omics. *Trends Biotechnol.* **38**(9), 1007–1022 (2020).

19. Yeo, H. C. & Selvarajoo, K. Machine learning alternative to systems biology should not solely depend on data. *Brief. Bioinform.* **23**(6), bbac436 (2022).

20. Elavarasan, D. & Vincent, P. M. D. Crop yield prediction using deep reinforcement learning model for sustainable agrarian applications. *IEEE Access* **8**, 86886–86901 (2020).

21. Neftci, E. O. & Averbeck, B. B. Reinforcement learning in artificial and biological systems. *Nat. Machine Intell.* **1**(3), 133–143 (2019).

22. Averbeck, B. & O'Doherty, J. P. Reinforcement-learning in fronto-striatal circuits. *Neuropsychopharmacology* **47**(1), 147–162 (2022).

23. Treloar, N. J. et al. Deep reinforcement learning for optimal experimental design in biology. *PLOS Comput. Biol.* **18**(11), e1010695 (2022).

24. Helmy, M., Smith, D. & Selvarajoo, K. Systems biology approaches integrated with artificial intelligence for optimized metabolic engineering. *Metab. Eng. Commun.* **11**, e00149 (2020).

25. Zhou, Z., Li, X. & Zare, R. N. Optimizing chemical reactions with deep reinforcement learning. *ACS Central Sci.* **3**(12), 1337–1344 (2017).

26. Hitti, Y. et al. Growspace: A Reinforcement Learning Environment for Plant Architecture. in *Available at SSRN* 4329504.
27. Vecchiato, G. et al. A 2D model to study how secondary growth affects the self-supporting behaviour of climbing plants. *PLOS Comput. Biol.* **19**(10), e1011538 (2023).

28. Brockman, G. et al. OpenAI Gym. eprint: (2016). arXiv:1606.01540.

29. Raffin, A. et al. Stable-baselines3: Reliable reinforcement learning implementations. *J. Mach. Learn. Res.* **22**(1), 12348–12355 (2021).

30. Schulman, J. et al. Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347 (2017).

31. Gianoli, E. The behavioural ecology of climbing plants. in *AoB Plants* **7** (2015).

32. Sutton, R. S. & Barto, A. G. *Reinforcement Learning: An Introduction* (MIT Press, 2018).

33. Vecchiato, G., Palladino, M., & Marcati, P. An optimal control approach to the problem of the longest self-supporting structure. In preparation.

34. Goriely, A. *The Mathematics and Mechanics of Biological Growth* Vol. 45 (Springer, 2017).

35. Agostinelli, D. et al. Nutations in growing plant shoots: The role of elastic deformations due to gravity loading. *J. Mech. Phys. Solids* **136**, 103702 (2020).

36. Goodno, B.J., & Gere, J.M. Mechanics of Materials. *Cengage Learn*. (2020).

37. Moulia, B., Douady, S. & Hamant, O. Fluctuations shape plants through proprioception. *Science* **372**(6540), eabc6868 (2021).

38. Chauvet, H. et al. Revealing the hierarchy of processes and time-scales that control the tropic response of shoots to gravi-stimulations. *J. Exp. Bot.* **70**(6), 1955–1967 (2019).

## Author contributions

Conceptualization & Formal analysis: L. Nasti, G. Vecchiato, M. Palladino, P. Marcati. Simulations & Experiments: L. Nasti. Software Development: L.Nasti. Writing original draft: L. Nasti, G. Vecchiato. Writing - review & editing: L. Nasti, G. Vecchiato. Data curation: P. Heuret, N.P. Rowe. Funding acquisition: P. Marcati, N.P. Rowe.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-024-62147-3.

**Correspondence** and requests for materials should be addressed to L.N.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024