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

HAL
open science

# Stronger wind, smaller tree: Testing tree growth plasticity through a modeling approach

Haoyu Wang, Jing Hua, Mengzhen Kang, Xiujuan Wang, Xing-Rong Fan, Thierry Fourcaud, Philippe de Reffye

► To cite this version:

Haoyu Wang, Jing Hua, Mengzhen Kang, Xiujuan Wang, Xing-Rong Fan, et al. Stronger wind, smaller tree: Testing tree growth plasticity through a modeling approach. Frontiers in Plant Science, 2022, 13, 10.3389/fpls.2022.971690. hal-03888028

HAL Id: hal-03888028

https://hal.inrae.fr/hal-03888028

Submitted on 7 Dec 2022

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_2.png" /></figure>

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

OPEN ACCESS

EDITED BY

Hartmut Stützel,
Leibniz University Hannover, Germany

REVIEWED BY

Jan Graefe,
Leibniz Institute of Vegetable and
Ornamental Crops, Germany
Félix P. Hartmann,
Institut National de recherche pour
l'agriculture, l'alimentation et
l'environnement (INRAE), France

*CORRESPONDENCE

Mengzhen Kang
mengzhen.kang@ia.ac.cn
Xiujuan Wang
xiujuan.wang@ia.ac.cn

SPECIALTY SECTION

This article was submitted to
Plant Biophysics and Modeling,
a section of the journal
Frontiers in Plant Science

RECEIVED 17 June 2022

ACCEPTED 05 October 2022

PUBLISHED 10 November 2022

CITATION

Wang H, Hua J, Kang M, Wang X, Fan X-R, Fourcaud T and de Reffye P (2022) Stronger wind, smaller tree: Testing tree growth plasticity through a modeling approach.
*Front. Plant Sci.* 13:971690.
doi: 10.3389/fpls.2022.971690

COPYRIGHT

© 2022 Wang, Hua, Kang, Wang, Fan, Fourcaud and de Reffye. This is an open-access article distributed under the terms of the [Creative Commons Attribution License (CC BY)](https://creativecommons.org/licenses/by/4.0/). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Stronger wind, smaller tree:
Testing tree growth plasticity
through a modeling approach

Haoyu Wang<fnref n="1" />,², Jing Hua<fnref n="1" />,², Mengzhen Kang<fnref n="1" />,³*,
Xiujuan Wang<fnref n="1" />,²*, Xing-Rong Fan⁴, Thierry Fourcaud<fnref n="5" />
and Philippe de Reffye<fnref n="5" />

Plants exhibit plasticity in response to various external conditions, characterized by changes in physiological and morphological features. Although being non-negligible, compared to the other environmental factors, the effect of wind on plant growth is less extensively studied, either experimentally or computationally. This study aims to propose a modeling approach that can simulate the impact of wind on plant growth, which brings a biomechanical feedback to growth and biomass distribution into a functional-structural plant model (FSPM). Tree reaction to the wind is simulated based on the hypothesis that plants tend to fit in the environment best. This is interpreted as an optimization problem of finding the best growth-regulation sink parameter giving the maximal plant fitness (usually seed weight, but expressed as plant biomass and size). To test this hypothesis in silico, a functional-structural plant model, which simulates both the primary and secondary growth of stems, is coupled with a biomechanical model which computes forces, moments of forces, and breakage location in stems caused by both wind and self-weight increment during plant growth. The Non-dominated Sorting Genetic Algorithm II (NSGA-II) is adopted to maximize the multi-objective function (stem biomass and tree height) by determining the key parameter value controlling the biomass allocation to the secondary growth. The digital trees show considerable phenotypic plasticity under different wind speeds, whose behavior, as an emergent property, is in accordance with experimental results from works of literature: the height and leaf area of individual trees decreased with wind speed, and the diameter at the breast height (DBH) increased at low-speed wind but declined at higher-speed wind. Stronger wind results in a smaller tree. Such response of trees to the wind is

<footnote n="1">¹The State Key Laboratory of Management and Control for Complex Systems, Institute of Automation, Chinese Academy of Sciences, Beijing, China, ²Beijing Engineering Research Center of Intelligent Systems and Technology, Institute of Automation, Chinese Academy of Sciences, Beijing, China, ³School of Artificial Intelligence, University of Chinese Academy of Sciences, Beijing, China, ⁴Engineering Research Centre for Waste Oil Recovery Technology and Equipment, Ministry of Education, Chongqing Technology and Business University, Chongqing, China,</footnote>

<footnote n="5">⁵CIRAD, AMAP, Univ Montpellier, CIRAD, CNRS, INRAE, IRD, Montpellier, France</footnote>
realistically simulated, giving a deeper understanding of tree behavior. The
result shows that the challenging task of modeling plant plasticity may be
solved by optimizing the plant fitness function. Adding a biomechanical model
enriches FSPMs and opens a wider application of plant models.

KEYWORDS

functional-structural plant model, mechanical model, critical wind speed, tree
breakage, optimization, thigmomorphogenesis

# 1 Introduction

Environmental conditions, such as light, temperature, and
humidity, influence the physiological processes and structural
development of trees. Wind, being almost ubiquitous in nature,
affects tree growth, reproduction, and even survival (Moulia et al.,
2006; de Langre, 2008; Lopez et al., 2011). The “compensating
mechanism” whereby advantageous modifications of a tree’s
phenotype can emerge to adapt to the wind conditions has
already been discovered (Whitehead, 1962; Whitehead & Luti,
1962). The term thigmomorphogenesis has been coined, which
refers to a phenomenon that trees demonstrate plasticity in
reaction to the external mechanical stimulus (Jaffe, 1973; Jaffe &
Forbes, 1993; Telewski & Pruyn, 1998). Considerable support to
this hypothesis has been provided by experimental observations:
when trees were exposed to wind forces or other biomechanical
perturbations, total dry weight, in particular stem dry weight,
decreased, and biomass allocation among organs altered
(Bonnesoeur et al., 2016). A reduction in both stem height and
leaf area, but an increase in diameter (Whitehead & Luti, 1962;
Jaffe & Forbes, 1993; de Langre, 2008), has been observed.
Allometry of axes, e.g., the ratio between length and radius, is
also affected. Wind sways induce strain force on trees, and it is
demonstrated that plants can sense strains they are subjected to
(Coutand & Moulia, 2000; Coutand et al., 2014). In some
circumstances, the secondary growth (increase in the stem
diameter) is more sensitive to changes in environmental
conditions than the primary growth (increase in the stem
length) (Collet et al., 2001).

How does a plant adjust its behavior accordingly? A plant
simulation system can be a helpful way to support or test
biological hypotheses (Trewavas, 2017; Calvo et al., 2019). The
functional-structure plant model (FSPM) is a key tool for
biomechanical computation as it provides detailed information
on plant structure and takes into account the interaction
between plant architecture and physiological processes during
plant growth (Fourcaud et al., 2008). Typically, FSPMs include
two components: a structural model and a functional model. The
one-directional or bidirectional interaction between these two
sub-models is necessarily a feature of FSPMs (Perttunen et al.,

1996; de Reffye et al., 2020). The bending of a tree stem has been
under study when a very early coffee tree structure was simulated
(de Reffye, 1976). The mechanics of orientation of a plant stem
due to the changes in loading and stresses have been studied with
biomechanical models (Fournier et al., 1994; Almeras et al.,
2002; Taylor-Hell et al., 2005). Considering the dynamic growth
of trees, an incremental biomechanical model has been
introduced to calculate the deformation of a single growing
stem (Fourcaud et al., 2003; Fourcaud & Lac, 2003), simulating
the influence of wind and self-weight. The modeling work has
been done at the scale of the stem (Coutand et al., 2011; Guillon
et al., 2012) or even the whole plant (Fourcaud et al., 2003).
However, the simulation of tree thigmomorphogenesis under a
wind environment is hardly seen. The feedback between
mechanics and plant growth has rarely been taken into
account in growth models at the architectural scale (Eloy
et al., 2017).

This study aims to explore tree growth plasticity under a
wind environment in a mechanistic way through a modeling
approach. To achieve this goal, firstly, a FSPM GreenLab (de
Reffye et al., 2020) is coupled with a biomechanical model
(Ancelin et al., 2004b) to simulate the effect of wind on
biomechanical variables inside a full tree structure; secondly, a
sensitivity analysis is performed to evaluate the effect of key
parameters controlling the secondary growth; thirdly, tree
plasticity is simulated through parameter optimization for
different wind speeds searching the best parameter values that
benefit tree growth. Thigmomorphogenesis is modeled as an
emergent result of optimization in biomass reallocation.

The validation of this theoretical result can be made by
experiments with different wind environments since parameters
for secondary growth can be estimated inversely from measured
data (Letort et al., 2008). However, as the sampling in such
experiments is very costly, an in silico experiment remains an
appealing tool to test hypotheses and inspire new thought. In a
former study, the optimization of wood yield with the constraint
of stability has been conducted (Qi et al., 2009), which also
attempts to find the best parameter for secondary growth. The
current work differs from it in several aspects: first, the wind
environment is considered, while the former is limited to gravity.
Secondly, this work couples the biomechanical and growth models, so that the incremental growth of trees is considered, while the former evaluates the mechanical properties of a straight adult tree. Thirdly, the biomechanical properties for all stems, including branches, are computed, while the previous work was limited to the main stem. The difference also lies in the constraint of the biomechanical stability of the tree.

# 2 Material and methods

## 2.1 System overview

The components of the whole system are shown in Figure 1. The plant model simulates step by step the tree structure and growth driven by a source-sink balance, giving the topological and geometrical tree structure with detailed information of the size and weight of individual organs (leaves, internodes, etc.). The wind environment model provides the wind profile and corresponding drag force. As this work is dedicated to wind, other environmental factors like light or temperature are not considered here, regarded as non-limiting. The biomechanical model computes the average force, moment, and new position of each internode under wind and gravity, for each growth step. The Factor of Safety Model checks whether the critical point of breakage is reached for every stem. The optimization model simulates plant plasticity by searching parameter values to maximize the objective function of tree growth while keeping the main stem unbroken.

For a plant growth model, it is necessary to clarify at which level of detail the model is built. Here, the temporal scale corresponds to a growth cycle (GC), i.e., the time it takes to create a new growth unit (GU, Figure 2A), which is a year for temperate trees. The spatial scale corresponds to a phytomer,

composed of an internode and its axillary leaves. A GU is composed of one or more phytomers of the same feature.

The integration of a dynamic plant growth process for plants and biomechanical model is as follows. At each GC, the organogenesis and organ growth take place, which creates new organs and increases the existing stem diameter (Figure 2B). A biomechanical model is called to compute force and moment as well as the new stem position caused by the increment in self-weight (Figure 2C). The deformation caused by self-weight is regarded as permanent. Given wind load, the wind-induced force leads to new destabilizations, and the biomechanical model is called again until a new equilibrium status, as well as the displacement of each internode, is reached (Figure 2D). However, after the withdrawal of wind, the plant recovers its original position as one can see in nature. This is achieved by running again the biomechanical module after the removal of wind force. This step is necessary as branch breakage induced by wind may have taken place, which influences self-weight distribution. In the next growth cycle, the plant development continues based on the deformed or broken structure of the last cycle.

## 2.2 Biomechanical model

Stem deformation is calculated at each internode k using the Incremental Transfer Matrix Method (ITMM) (Ancelin et al., 2004b), as in Eq. 1.

$$ \left\{ \begin{aligned} D_{k+1} &= T_{11}D_k + T_{12}(C_k - S_k) + M_1F + M_a \\ S_{k+1} &= T_{21}D_k + T_{22}(C_k + S_k) + M_2F \end{aligned} \right. \qquad (1) $$

where D is the displacement vector, which includes the positive translations of the cross-section center (m) and the

<figure><img src="image_4.png" /><figcaption>FIGURE 1
Schematic diagram of the modeling system.</figcaption></figure>
<figure><img src="image_5.png" /></figure>

FIGURE 2

Illustration of tree development, primary and secondary growth, as well as the simulation of wind and gravity effect during a growth cycle.
(A) plant top with two physiological ages (shown with different textures) at GC *n* - 1; (B) organogenesis, primary growth (blank rectangles) and secondary growth (new rings on existing rectangles) at GC *n*; (C) permanent deformation caused by increment in gravity; (D) temporary deformation caused by wind.

positive rotations of the cross-section normal (radians) (see components in **Supplement A**), S is the stress and moment vector; C is the vector of the sum of all concentrated forces and moments at the branching node, with the incremental weight of the branch computed by the plant model; F is the vector of linear forces, caused by both the incremental self-weight and wind load (Eq. 9); and T and M are transformation matrices (see Eqs. S3-S6 and Eqs. S7-S12 in **Supplement A**). Ma denotes the maturation strains (MS) at the periphery of stems, associated with the formation of reaction wood (see Eq. S13 in **Supplement A**). The link between the biomechanical model and plant model lies in the concentrated forces and moments at the branching node (C), and the self-weight of internode (F), as the FSPM gives the 3D structure and the concentrated force of each sub-structure.

During the dynamic growth of plants, in each growth cycle, tree secondary growth takes place by adding a new layer to the existing stems. In the biomechanical model, a tree stem is thus regarded as a sequence of multilayer three-dimensional beam elements (Ancelin et al., 2004b). The biomechanical calculation consists of two steps. Firstly, compute the new incremental loads and the moments (caused by self-weight) on each internode. This is done using the idea of sub-structure (Kang et al., 2008), starting from the twigs of the highest branching order until the trunk (main stem). In ITMM, only the increment of weight must be considered in a growth step. Leaves, fruits, and flowers are considered as concentrated forces and moments for each internode. Dead or pruned branches are removed from the insertion point on the parent branch. Secondly, compute the new incremental displacement based on the last equilibrium, until the top of the main stem. Once this process is completed, the result is validated by examining the biomechanical boundary condition (BC) at the free top of the tree (Ancelin et al., 2004a),

so that there is hardly any force or moment at each tip. If the BCs are not met, the last two steps are repeated.

When the wind is imposed, the displacement of the internode will take place, which changes the gravity center. With this increment of load, at each growth cycle, ITMM needs to be run multiple times to reach a new equilibrium, as if the plant swings before becoming stable. Similarly, when the wind calms down, iterated computation is called again. As a result, at each growth cycle, the equilibrium status of the tree with and without wind is updated. The iteration number can be set respectively for the main stem and branches.

## 2.3 Incremental weight in plant structure

The incremental weight of individual organs during plant growth is computed using the GreenLab model. GreenLab is a generic FSPM that has been applied for beech tree (Fagus sylvatica L.) (Letort et al., 2008), poplar tree (Salicaceae) (Yang et al., 2011), and Chinese pine (Pinus tabulaeformis Carr.) (Guo et al., 2012). For trees, a model with a finer temporal scale with stochastic behavior has been proposed (Kang et al., 2018), but here a deterministic model based on the growth cycle is adopted, where internodes created during the same GU have the same size. We recall here mainly components for secondary growth. See Supplement B for a more complete presentation, especially the primary growth.

For biomechanical computation, the output needed from GreenLab include the length and diameter of individual internodes; the weight, position, and gravity center of the aggregated branch; and the initial geometrical shape of the tree structure for computing the gravity center. The last is computed
from the topological structure, the size of organs, and the insertion and phyllotaxy angle, which is not detailed here. The poplar (*Populus*) tree is selected for this study. Parameter settings used in this study are shown in Table 1.

### 2.3.1 Weight of organs and substructure

The plant growth is based on source-sink regulation. Based on the pipe model, the demand for the secondary growth ($D_{sec}^{layer}$) is proportional to the number of leaves in the last cycle $N_p^{b}(n-1)$, using Eq. 2:

$$D_{sec}^{layer}(n) = S_{layer} \sum_{p=1}^{p_m} N_p^{b}(n-1) \quad (2)$$

where $S_{layer}$ denotes the sink strength of a new layer contributed by an individual leaf. It controls the proportion of produced biomass allocated to the secondary growth at each GC. Subscript *p* is the physiological age (PA), from 1 to maximum physiological age $P_m$. Superscript *b* denotes the leaf blade. Symbol *n* denotes growth cycle *n*. In GreenLab, new rings of wood formed around the stem periphery. Each ring on a stem is contributed by numerous layers. Rings are modeled in two steps. Firstly, calculate the amount of biomass allocated to the secondary growth at GC *n* at the whole-plant scale, denoted by $Q_{sec}^{layer}(n)$, as in Eq. 3:

$$Q_{sec}^{layer}(n) = D_{sec}^{layer}(n) \frac{Q(n)}{D(n)} \quad (3)$$

Secondly, $Q_{sec}^{layer}(n)$ is allocated to each internode with two possible modes. Let parameter $\lambda$ describe the way of biomass allocation associated with the secondary growth for each phytomer; the incremental biomass at GC *n* for an internode of PA *p* that showed GC *j* is described in Eq. 4:

$$q_{sec,p}^{layer}(j, n) = \left( \frac{1 - \lambda}{D_1(n)} + \frac{\lambda \cdot s_p^a(j, n)}{D_2(n)} \right) \cdot l_p(j) \cdot P_p^r \cdot Q_{sec}^{layer}(n) \quad (4)$$

See Eq. S22 for demand for primary growth $D_1(n)$ and secondary growth $D_2(n)$, as well as other parameters. If $\lambda = 0$, it means that at each cycle the biomass for the layer is uniformly allocated to internodes of the same PA; tapering can appear because older internodes have more layers. If $\lambda = 1$, the allocation is also proportional to the total number of living leaves above the internode; since upper internodes have fewer leaves above, thus tapering is more obvious compared to the case of $\lambda = 0$. The weight of an individual internode that is born in cycle *j* is thus the result of primary and secondary growth:

$$q_p^i(j, n) = q_{pri,p}^i(j) + \sum_{i=j+1}^{n} q_{sec,p}^{layer}(j, i) \quad (5)$$

Two key parameters that control the biomass allocation to the secondary growth are (1) parameter $S_{layer}$ controlling the amount of biomass allocated to the secondary growth at the whole-plant level (Eq. 2) and (2) parameter $\lambda$ controlling the allocation of biomass for secondary growth among all internodes (Eq. 4) and, accordingly, stem shape, as shown in Supplementary Figure S1.

### 2.3.2 Weight of plant and substructures

The weight of a branch (sub-structure) of PA *p* is the total weight of its bearing axis (stem) and all axillary axes (Kang et al., 2008). With the sub-structure method, the summing up starts from the branch of the highest PA (no sub-branch) until the main stem of PA 1. Similarly, the incremental weight of the branch can be obtained.

The total stem weight of the structure of the plant is the sum of the weights of the main stem and all branch stems, which is one of the object functions for optimization:

<table><caption>TABLE 1 Parameter values of poplar tree.</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Value</th></tr></thead><tbody><tr><td>S<sub>p</sub></td><td>Total ground projection area available of the crown of plant; see Eq. S15 in Supplement B.</td><td>12 m<sup>2</sup></td></tr><tr><td>P<sub>m</sub></td><td>Maximum physiological age (PA); see Eq. 2.</td><td>4</td></tr><tr><td>e</td><td>Specific leaf weight; see Eq. S23 in Supplement B.</td><td>0.023 g·cm<sup>-2</sup> (Yang et al., 2011)</td></tr><tr><td>P<sub>p</sub><sup>b</sup></td><td>Sink strength of leaf blade of PA p; see Eq. S16 in Supplement B.</td><td>1, 0.66, 0.32, 0.2, from PA 1 to PA P<sub>m</sub> (Yang et al., 2011)</td></tr><tr><td>P<sub>p</sub><sup>i</sup></td><td>Sink strength of internode pith of PA p; see Eq. S16 in Supplement B.</td><td>0.875, 0.5, 0.3, 0.1, from PA 1 to PA P<sub>m</sub></td></tr><tr><td>P<sub>p</sub><sup>r</sup></td><td>Relative sink strength of layer with PA p; see Eq. S22 in Supplement B.</td><td>0.8, 0.05, 0.05, 0.05</td></tr><tr><td>α</td><td>Internode allometry parameter; see Eq. S19 in Supplement B.</td><td>-0.28, -0.36, -0.26, -0.06 from PA 1 to PA P<sub>m</sub> (Yang et al., 2011)</td></tr><tr><td>β</td><td>Internode allometry parameter; see Eq. S19 in Supplement B.</td><td>6.51, 4.22, 3.77, 5.20 from PA 1 to PA P<sub>m</sub> (Yang et al., 2011)</td></tr><tr><td>P<sub>p</sub><sup>pih</sup></td><td>Density of pith; see Eq. S18 in Supplement B.</td><td>0.36 g·cm<sup>-3</sup> (Kord et al., 2010)</td></tr><tr><td>P<sub>p</sub><sup>layer</sup></td><td>Density of layer; see Eq. 7</td><td>0.36-0.48 g·cm<sup>-3</sup> (Kord et al., 2010)</td></tr><tr><td>MOR</td><td>Modulus of rupture</td><td>45 MPa (Grotta et al., 2005b)</td></tr></tbody></table>
$$Q_{\text{stem}} = \sum_{p=1}^{P_m} \sum_{j=1}^{n} N_p^i(j) q_p^i(j, n) \quad (6)$$

The number of internodes that are created at cycle j, , is updated in each cycle considering both organogenesis and possible stem breakage.

## 2.4 Size of organs

Considering that wood density is dependent on the radial position and height of the internode, parameter $ρ_p^{\text{layer}}(j, i)$ is set to describe the density of the layer in GC i for internode of PA p created at GC j. Then, at cycle i, the section area of layer added to such internodes is computed as

$$s_{\text{sec},p}^{\text{layer}}(j, i) = \frac{q_{\text{sec},p}^{\text{layer}}(j, i)}{\rho_p^{\text{layer}}(j, i) * l_p(j)} \quad (7)$$

The radius of internode is finally the result of layer accumulation:

$$r_p^i(j, n) = \sqrt{s_p^i(j, n)/\pi} = \sqrt{[s_p^{\text{pith}}(j) + \sum_{i=j+1}^n s_{\text{sec},p}^{\text{layer}}(j, i)]/\pi} \quad (8)$$

where $l_p(j)$ and $s_p^{\text{pith}}(j)$ are the length and section area of pith from primary growth (see Supplement B).

## 2.5 Wind-induced force

According to literature (Mayhead, 1973; Niklas & Spatz, 2000; Ancelin et al., 2004a; Diener et al., 2010), the wind-induced drag force $Fw$ by wind on the kth element at the height z aboveground can be described as in Eq. 9:

$$F_{w_k}(z) = \frac{1}{2} \rho_{\text{air}} C_d A_k u(z)^2 \\ A_k = 2l_k(z)r_k(z) |\sin \chi| \quad (9)$$

where $ρ_{\text{air}}$ is the air density (1.226 kg·m⁻³ when the air temperature is 15°C). $C_d$ (dimensionless) is the drag coefficient, assumed to be 0.25 according to Koizumi et al. (2010). Differing from Sellier et al. (2008) where the exposed area is estimated from the biomass ratio, here $A_k$ (m²) is the frontal area of the kth internode, i.e., the projected area of the internode against the wind direction. It is computed with the length, diameter of the internode, and angle to the vertical direction of the internode, where l and r are the length and radius of the internode, respectively; X is the angle between this internode with the wind direction. $u(z)$ is the wind speed at vertical height z which is given by a logarithmic wind profile.

Wind profile refers to the vertical distribution of wind speed. There are different wind profiles in the literature describing the horizontal mean wind speed according to the height

aboveground (Niklas & Spatz, 2000). Some researchers (Irvine et al., 1997; Ancelin et al., 2004a) distinguished the wind profiles in forests located inside the stand and at the stand edge, including above and below the canopy. For trees standing in an open terrain, as assumed in this paper, a logarithmic wind profile is used (Ancelin et al., 2004b), which has provided the best fit for wind speeds empirically recorded within and just above the canopy (Niklas & Spatz, 2000).

$$u(z) = u(h_0) \frac{\ln(z/z_o)}{\ln(h_0/z_o)} \quad (10)$$

where $u(z)$ (m·s⁻¹) is the wind speed at vertical height z; $h_0$ is the standard reference height of 10 m in Meteorology (Davenport, 1962); $z_o$ is the roughness length at stand edge (m). Ratio $z_o/h_0$ is set to 0.06 (Peltola & Kellomäki, 1993; Ancelin et al., 2004a).

## 2.6 Factor of safety

It is assumed that stem breakage takes place when its maximum longitudinal stress exceeds the Modulus of Rupture (MOR) of wood fibers (Jones, 1983; Gardiner et al., 2000). According to Jones (1983); Gardiner et al. (2000); Grotta et al. (2005), and Fjeld (2012), it is found that individual trees differed significantly in their MOR (range: 37–148 MN/m²). The maximum absolute value of longitudinal stress is found at the periphery of stem cross sections and depends on the bending moment (M, N·m) of the stem and the third power of the stem diameter (Jones, 1983; Gardiner et al., 2000), as in Eq. 11:

$$\text{Stress}_k = \frac{32M_k}{\pi d_k^3} \quad (11)$$

Differing from Gardiner et al. (2000) who assumed that the stress is constant and thus calculated the stress only at breast height, here for each internode this stress is computed. The bending moment $M_k$ caused by the wind results from Eq. 1 (a component of $S_k$). $d_k$ is the diameter of the kth stem internode as calculated in Eq. 8, and Stress$_k$ (Pa) is the longitudinal stress of internode k.

The critical wind speed (CWS) is defined as the minimum wind speed causing the breakage of the trunk, with the given parameter set. Higher CWS means that the tree is more wind-resistant.

The wind resistance can also be indicated by the top stem deflection of the trunk, defined as the angle between the trunk's top and vertical direction.

## 2.7 Objective function

As mentioned above, the adaptation of trees to varying circumstances is modeled as an optimization problem. The
optimization objective of the acclimatization is set to be the wood fresh weight $Q_{\text{stem}}$ (total internode weight, Eq. 6) and the height of tree $H_{\text{tree}}$. The reason for adding tree height as one of the objectives is that reaching a certain height is beneficial for light interception. Since in stronger winds the tree is more prone to breakage, thus it is expected that the tree's reaction to the wind is to increase the radial growth. However, over-allocation to secondary growth will diminish the photosynthesis production and tree growth; thus, there is a balance between primary and secondary growth. The aim of optimization is to find the best balance. We choose the sink strength parameter for the secondary growth $S_{\text{layer}}$ (Eq. 2) as the control variable.

Stem breakage under loads includes both trunk breakage and branch breakage (Lopez et al., 2011). Breakage on the trunk is likely to be fatal, while tree survival is often possible after branch breakage. Thus, here trunk breakage is considered a constraint condition.

In summary, this optimization problem with constraints is given by Eq. 12:

$$ \left\{ \begin{array}{l} \text{Max } H_{\text{tree}}(S_{\text{layer}}) \\ \text{Max } Q_{\text{stem}}(S_{\text{layer}}) \\ \text{subject to } \text{Stress}_{\text{trunk}} < \text{MOR} \end{array} \right. \qquad (12) $$

The value of MOR is 45 MPa (Table 1). The multi-objective optimization algorithm of the Non-determined Sorting Genetic Algorithm II (NSGA-II) (Deb et al., 2000) from PAGMO (a C + scientific library for massively parallel optimization) (Biscani & Izzo, 2020) is adopted. NSGA-II is one of the most popular multi-objective genetic algorithms. It has the advantages of fast-running speed and good convergence of the solution set, making it a benchmark of the performance of other multi-objective optimization algorithms. The population size is set as 100, the

evolutionary algebra as 100, the crossover probability as 0.7, and the mutation probability as 0.4.

# 3 Results

## 3.1 The equilibrium of tree under wind load

As mentioned above, the biomechanical calculation (Eq. 1) consists of two steps and requires multiple computational iterations to reach equilibrium. Under wind, the stem deformation stabilizes after several iterations or simply breaks under strong wind (Figure 3). At a wind speed of 13 m/s, the tree sways during about 10 iterations before the stabilization. When the wind speed is 15 m/s, trunk breakage takes place after the first iteration. With more iteration steps, the computational time increases linearly. With 10 iterations, at age 8 of the plant which means 8 years, it takes about 20 s for running the whole model in a PC of Core™ i7-8550U CPU 1.80 GHz, 16 GB RAM. Let $N_p$ be the iteration number for the substructure of PA p. According to the result, in the following computations, the iteration number of ITMM is set to 10 since the displacement error becomes negligible.

## 3.2 Numerical results

### 3.2.1 Force and moment

To visualize how the system works, the tree structure of GC 8 is simulated with parameter set $S_{\text{layer}} = 4$, $\lambda = 0.1$, under wind speed 15 m/s. The result is shown for the first iteration before trunk breakage. The displacement, force (Figures 4A-C, the first row),

<figure><img src="image_6.png" /><figcaption>FIGURE 3</figcaption></figure>

The stem deflection angle of the trunk under wind speeds of 13 m/s and 15 m/s. When wind speed is 13 m/s, it takes about 10 iterations before the stabilization. When wind speed is 15 m/s, trunk breakage takes place after the first iteration.
<figure><img src="image_7.png" /><figcaption>FIGURE 4
Computed force (A-C), moment (D-F) and stress<sub>k</sub> (G-I) of x-, y-, and z-component for the plant under wind speed 15 m/s, shown with different colors. S<sub>layer</sub> = 4, λ = 0.1.</figcaption></figure>

moment (Figures 4D-F, the second row), and stress<sub>k</sub> (Figures 4G-I,
the third row) of each stem under a wind environment (which
blows in the x-direction from left to right) are shown through tree
shape and color, for the x-, y-, and z-components, respectively. The
force is a vector. When the direction of the force is inconsistent with
the direction of the coordinate axis, the force takes a negative value.

3.2.2 Effect of wind speed

To understand how wind speed can influence tree breakage,
the tree shapes are shown in Figure 5. Here, S<sub>layer</sub> and λ are set to
6 and 0.1, respectively. With increasing wind speed, more
branches are lost because of breakage. The effect is evident on
the top of the trunk.
<figure><img src="image_8.png" /><figcaption>FIGURE 5</figcaption></figure>

The effect of wind speed on tree breakage: 13m/s (A), 15m/s (B), 18m/s (C), respectively. Result is shown under calm conditions. The remaining trunk biomass is 170.78, 159.85 and 120.58 kg, respectively. The circles indicate the points where the trunk breakage takes place. The plant age is 8, S$_{layer}$ = 6, λ = 0.1.

3.3 Parameter sensitivity

3.3.1 Parameter S$_{layer}$

Before parameter optimization, it is useful to understand the
model behavior using a sensitivity analysis. According to the
GreenLab model, the biomass allocation to the secondary growth,
regulated by parameter S$_{layer}$ (Eqs. 2, 3), has twofold effects. On the

one hand, the increase in biomass allocation to secondary growth
increases the stem biomass. On the other hand, it can inhibit
primary growth, thus decreasing the total leaf area and biomass
production. In Figure 6, the stem biomass and top stem deflection
for a trunk with different Slayer values are simulated. As expected, the
peak value of the biomass of tree and trunk exists, which is
independent of the wind. With a lower Slayer value, stem breakage

<figure><img src="image_9.png" /><figcaption>FIGURE 6</figcaption></figure>

Influence of the sink strength for the secondary growth (S<sub>layer</sub>) on the wood biomass Q<sub>stem</sub> and top step deflection for the trunk, under wind speed 15 m/s with λ = 0.4.
is more evident; a sudden drop for S<sub>layer</sub><8 means trunk breakage. Supplementary Figure S2 shows 3D tree shapes with different S<sub>layer</sub> values, where broken branches are removed. With a higher S<sub>layer</sub> value, the tree structure is more complete but smaller, because of the smaller leaf area for photosynthesis.

The critical value of wind speed (CWS) gradually increases with S<sub>layer</sub> (Figure 7). It is because when more biomass was allocated to the secondary growth, the tree is easier to resist the breakage induced by wind. This shows from another point of view than the step deflection angle that with a bigger S<sub>layer</sub>, the tree is more resistant to wind.

### 3.3.2 Parameter $\lambda$

Differing from the conclusion that $\lambda$ should be equal to 1 to obtain the maximal wood biomass (Qi et al., 2009), here the influence of $\lambda$ is more complex, as shown in Figure 8. At a lower S<sub>layer</sub> value, e.g., S<sub>layer</sub> = 0.5, the trunk is thin and the diameter of the tree top is too slim to support itself. Even very small wind can cause a breakage due to self-weight. The position with the maximum value of the breakage value also moves downward along the trunk. At a higher S<sub>layer</sub> value, e.g., S<sub>layer</sub> = 4, the critical wind bore by the trunk becomes stronger and becomes indifferent to the $\lambda$ value. Table 2 shows the force and moment given MOR = 45 MN/m². Supplementary Figure S3 shows the 3D tree shapes under different $\lambda$, where the largest stress value on the trunk Stress<sub>Trunk</sub> is marked with a black circle. Only the trunk and the primary branch structure are shown. $\lambda$ does not change the height of the tree as it matters only the allocation of biomass for secondary growth, but it influences the stem shape and accordingly breakage. According to the numerical results, $\lambda$ has little effect on optimization results; thus, it is not considered as a control parameter for tree reaction (set to $\lambda$ = 0.1).

## 3.4 Plant plasticity simulation

### 3.4.1 Pareto optimal frontier

Tree plasticity is simulated by searching the best parameter value of S<sub>layer</sub> to maximize the tree biomass and tree height while keeping the trunk unbroken, given the wind speed. Because of multi-objective optimization, the result is no more a single point but a point set, called Pareto set or Pareto frontier (Figure 9). In general, when the wind speed is smaller than 16 m/s, the distribution of points in the Pareto boundary coincides with each other, which means the reaction tree is close to the relatively small wind. When the wind speed is high, as seen at the wind speed of 16 m/s, the points in the Pareto optimal front are separated from the previous ones, which means a different reaction. In both situations, according to the Pareto sets, the optimal tree can be a taller tree with smaller biomass, or a shorter tree with bigger biomass. In general, the tree height shows an obvious downward trend with wind speed. The total weight does not change significantly when the wind speed is low but has wider distribution with the increase in wind speed.

### 3.4.2 Optimal parameter values

Each dot in the Pareto frontier corresponds to an optimal parameter, so the box diagram is used to show the dispersion degree of the optimal S<sub>layer</sub> for each wind speed (Supplementary Figure S4). When the wind speed is small, the box is small, which means the parameters are more concentrated. When the wind speed increases to 16 m/s or bigger, the maximum deviation of the “outlier” from the box gradually increases, indicating that the tree tends to allocate more biomass to the secondary growth to reinforce its structure under the wind. With the increase in wind speed, the median value of each box has a gradual upward trend.

<figure><img src="image_10.png" /><figcaption>FIGURE 7
The critical wind speed (CWS) increases with S<sub>layer</sub> values with $\lambda$ = 0.4.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 8
Effect of $\lambda$ on the critical wind speed (CWS) with different $S_{layer}$ values.</figcaption></figure>

Table 3 shows the optimal parameter values under different wind speeds according to the box diagram (Supplementary Figure S4) by taking representative values (the maximum, minimum, median, and the ratios of 10%, 25%, 75%, and 90%) of the box. In general, the trend is that the parameter value increased with wind speed.

### 3.4.3 Morphological result

To better analyze the morphological changes of trees reacting to different wind conditions, three points (representing three strategies) at the two ends and the middle of the Pareto frontier are selected. For each wind speed, the corresponding $S_{layer}$ value of the same strategy is shown in Table 4. The strategy I (0% H, 100% Q) in Table 4 represents

the point corresponding to the minimal tree height and the maximal biomass in the Pareto optimal front. This point can be considered as the tree focuses on pursuing weight goals. In this case, when the wind speed is small (less than 12 m/s), the change of $S_{layer}$ is not obvious, so it can be considered that it is not sensitive to the stimulation of the breeze. Strategy III (100% H, 0% Q) corresponds to the point with the largest tree height and the smallest biomass in the Pareto optimal frontier. It can be considered that in this case, trees focus on pursuing tree height but not the tree biomass. In this case, the $S_{layer}$ value is taken in a small range. Strategy II (50% H, 50% Q) corresponds to a middle point in the Pareto front representing this set of optimal solutions. Its selection method is as follows: (1) average the optimal solution set corresponding to a Pareto front; (2) select

<table><caption>TABLE 2 The force and moment given MOR 45 MN/m².</caption><thead><tr><th>λ</th><th>Stress<sub>tree</sub></th><th>Stress<sub>trunk</sub></th><th>POT</th><th>H</th><th>D</th><th>f<sub>x</sub></th><th>f<sub>y</sub></th><th>f<sub>z</sub></th><th>M<sub>x</sub></th><th>M<sub>y</sub></th><th>M<sub>z</sub></th></tr></thead><tbody><tr><td>0</td><td>16.06</td><td>3.94</td><td>1</td><td>0.00</td><td>20.37</td><td>245.53</td><td>0</td><td>-485.59</td><td>3.38</td><td>3271.61</td><td>0.20</td></tr><tr><td>0.1</td><td>17.44</td><td>3.12</td><td>77</td><td>1.62</td><td>20.47</td><td>218.73</td><td>0</td><td>-450.67</td><td>2.33</td><td>2623.95</td><td>-2.41</td></tr><tr><td>0.2</td><td>19.05</td><td>1.86</td><td>74</td><td>1.54</td><td>24.40</td><td>218.62</td><td>0</td><td>-468.49</td><td>3.43</td><td>2569.46</td><td>0.05</td></tr><tr><td>0.3</td><td>21.44</td><td>1.43</td><td>74</td><td>1.54</td><td>26.22</td><td>216.16</td><td>0</td><td>-477.19</td><td>2.96</td><td>2455.43</td><td>-0.05</td></tr><tr><td>0.4</td><td>26.26</td><td>1.14</td><td>74</td><td>1.54</td><td>27.92</td><td>213.35</td><td>0</td><td>-485.90</td><td>2.46</td><td>2340.48</td><td>-0.19</td></tr><tr><td>0.5</td><td>34.33</td><td>0.91</td><td>74</td><td>1.54</td><td>29.52</td><td>210.15</td><td>0</td><td>-494.60</td><td>1.90</td><td>2224.00</td><td>-0.42</td></tr><tr><td>0.6</td><td>49.20</td><td>0.74</td><td>74</td><td>1.54</td><td>31.04</td><td>206.52</td><td>0</td><td>-503.30</td><td>1.26</td><td>2105.41</td><td>-0.77</td></tr><tr><td>0.7</td><td>81.33</td><td>0.61</td><td>74</td><td>1.54</td><td>32.49</td><td>202.36</td><td>0</td><td>-512.01</td><td>0.51</td><td>1983.95</td><td>-1.38</td></tr><tr><td>0.8</td><td>167.62</td><td>0.50</td><td>74</td><td>1.54</td><td>33.88</td><td>197.50</td><td>0</td><td>-520.71</td><td>-0.41</td><td>1858.54</td><td>-2.53</td></tr><tr><td>0.9</td><td>456.40</td><td>0.42</td><td>74</td><td>1.54</td><td>35.21</td><td>191.60</td><td>0</td><td>-529.42</td><td>-1.46</td><td>1727.47</td><td>-4.91</td></tr><tr><td>1</td><td>140.26</td><td>0.12</td><td>39</td><td>0.61</td><td>39.01</td><td>206.62</td><td>0</td><td>-670.55</td><td>2.84</td><td>1017.44</td><td>3.84</td></tr></tbody></table>

The wind speed is 13.5 m/s. Stress<sub>tree</sub> represents the maximum stress value (MN/m² or MPa) inside the whole tree. POT represents the position that bears the maximum stress value on the trunk (Stress<sub>trunk</sub>, MN/m²), which is the number of internodes counted from the base to the top. H and D represent the height (m) and diameter (cm) of POT, respectively. $f_x$, $f_y$, and $f_z$ represent the x-, y-, and z-components of force (N) at POT, respectively. $M_x$, $M_y$, and $M_z$ represent the x-, y-, and z-components of the moment (N·m) at POT, respectively.
<figure><img src="image_12.png" /><figcaption>FIGURE 9
Pareto optimal front for different wind speeds.</figcaption></figure>

an optimal solution belonging to the Pareto set that is closest to the average value to represent the whole set of optimal solutions. This strategy is a compromise between tree height and weight. The $S_{layer}$ value gradually increases in strategy II.

Through the representative optimal result of strategy II, the plasticity in tree morphological features can be seen. With the increase in wind speed, the diameter at breast height (DBH) increased slightly with wind speed but declined significantly at a strong wind speed; the height and leaf area of the tree progressively declined, which is consistent with the experiments in the literatures (Whitehead, 1962; Telewski & Pruyn, 1998). The partitioning of biomass to the secondary growth increased,

especially those for the trunk. The morphology of the tree significantly varies with the gradual acceleration of wind speed. The tree structure corresponding to strategies I, II, and III changes with the wind speed, as shown in Figure 10. The corresponding optimal 3D shapes of simulated poplar trees for strategy II under different wind environments are shown in Figure 11.

## 4 Discussion

The interest in this work was initially inspired by the phenomena of thigmomorphogenesis (Hamant, 2013), which

<table><caption>TABLE 3 Optimized parameter values for different wind speeds.</caption><thead><tr><th>Wind speed (m/s)</th><th>Minimum</th><th>Error bars below box (10th percentile)</th><th>Lower boundary of the box (25th percentile)</th><th>Media (50th percentile)</th><th>Upper boundary of the box (75th percentile)</th><th>Error bars above box (90th percentile)</th><th>Maximum</th></tr></thead><tbody><tr><td>0</td><td>0.25</td><td>0.28</td><td>0.34</td><td>0.47</td><td>0.66</td><td>0.84</td><td>0.96</td></tr><tr><td>2</td><td>0.26</td><td>0.29</td><td>0.34</td><td>0.45</td><td>0.65</td><td>0.84</td><td>0.97</td></tr><tr><td>4</td><td>0.3</td><td>0.33</td><td>0.37</td><td>0.48</td><td>0.68</td><td>0.82</td><td>0.97</td></tr><tr><td>6</td><td>0.32</td><td>0.34</td><td>0.39</td><td>0.54</td><td>0.72</td><td>0.85</td><td>0.97</td></tr><tr><td>8</td><td>0.35</td><td>0.38</td><td>0.44</td><td>0.57</td><td>0.74</td><td>0.87</td><td>0.97</td></tr><tr><td>10</td><td>0.44</td><td>0.46</td><td>0.51</td><td>0.63</td><td>0.78</td><td>0.88</td><td>0.97</td></tr><tr><td>12</td><td>0.61</td><td>0.63</td><td>0.67</td><td>0.73</td><td>0.82</td><td>0.91</td><td>0.97</td></tr><tr><td>14</td><td>0.88</td><td>0.89</td><td>0.89</td><td>0.91</td><td>0.93</td><td>0.95</td><td>1.13</td></tr><tr><td>16</td><td>0.51</td><td>0.53</td><td>0.75</td><td>0.88</td><td>0.95</td><td>1</td><td>2.97</td></tr><tr><td>18</td><td>0.86</td><td>0.86</td><td>0.87</td><td>0.89</td><td>0.97</td><td>1.45</td><td>4.89</td></tr><tr><td>20</td><td>0.54</td><td>0.55</td><td>0.57</td><td>0.6</td><td>0.65</td><td>2.28</td><td>6.93</td></tr><tr><td>22</td><td>0.65</td><td>0.65</td><td>0.65</td><td>1.54</td><td>8.6</td><td>9.22</td><td>9.29</td></tr><tr><td>24</td><td>0.77</td><td>5.84</td><td>12.48</td><td>12.48</td><td>12.48</td><td>12.48</td><td>12.48</td></tr><tr><td>26</td><td>1.33</td><td>1.33</td><td>1.33</td><td>18.21</td><td>18.21</td><td>18.21</td><td>18.29</td></tr></tbody></table>
<table><caption>TABLE 4 Corresponding S<sub>layer</sub> value with different wind speeds under strategies I, II, and III.</caption><thead><tr><th rowspan="2">Wind speed (m/s)</th><th colspan="3">Strategy I:0% H, 100% Q</th><th colspan="3">Strategy II:50% H and 50% Q</th><th colspan="3">Strategy III:100% H, 0% Q</th></tr><tr><th>H<sub>tree</sub> (m)</th><th>Q<sub>stem</sub> (kg)</th><th>S<sub>layer</sub></th><th>H<sub>tree</sub> (m)</th><th>Q<sub>stem</sub> (kg)</th><th>S<sub>layer</sub></th><th>H<sub>tree</sub> (m)</th><th>Q<sub>stem</sub> (kg)</th><th>S<sub>layer</sub></th></tr></thead><tbody><tr><td>0</td><td>15.73</td><td>150.75</td><td>0.96</td><td>17.96</td><td>143.11</td><td>0.52</td><td>21.73</td><td>75.83</td><td>0.25</td></tr><tr><td>2</td><td>15.69</td><td>150.75</td><td>0.97</td><td>18.02</td><td>142.71</td><td>0.51</td><td>20.94</td><td>89.75</td><td>0.26</td></tr><tr><td>4</td><td>15.69</td><td>150.75</td><td>0.97</td><td>17.85</td><td>143.86</td><td>0.54</td><td>20.05</td><td>95.62</td><td>0.3</td></tr><tr><td>6</td><td>15.69</td><td>150.75</td><td>0.97</td><td>17.67</td><td>144.91</td><td>0.57</td><td>19.74</td><td>96.14</td><td>0.32</td></tr><tr><td>8</td><td>15.69</td><td>150.75</td><td>0.97</td><td>17.56</td><td>145.54</td><td>0.59</td><td>19.07</td><td>102.53</td><td>0.35</td></tr><tr><td>10</td><td>15.69</td><td>150.75</td><td>0.97</td><td>17.23</td><td>147.19</td><td>0.65</td><td>18.46</td><td>116.05</td><td>0.44</td></tr><tr><td>12</td><td>15.69</td><td>150.75</td><td>0.97</td><td>16.71</td><td>149.18</td><td>0.75</td><td>17.45</td><td>146.13</td><td>0.61</td></tr><tr><td>14</td><td>15.03</td><td>149.97</td><td>1.13</td><td>15.93</td><td>150.68</td><td>0.92</td><td>16.09</td><td>150.52</td><td>0.88</td></tr><tr><td>16</td><td>9.63</td><td>85.68</td><td>2.97</td><td>15.94</td><td>150.67</td><td>0.91</td><td>18.02</td><td>142.71</td><td>0.51</td></tr><tr><td>18</td><td>6.46</td><td>28.04</td><td>4.89</td><td>15.65</td><td>150.75</td><td>0.98</td><td>16.18</td><td>150.39</td><td>0.86</td></tr><tr><td>20</td><td>4.76</td><td>10.34</td><td>6.93</td><td>14.15</td><td>146.62</td><td>1.36</td><td>17.85</td><td>143.86</td><td>0.54</td></tr><tr><td>22</td><td>3.73</td><td>4.75</td><td>9.29</td><td>5.96</td><td>21.73</td><td>5.35</td><td>17.23</td><td>147.19</td><td>0.65</td></tr><tr><td>24</td><td>2.99</td><td>2.45</td><td>12.48</td><td>3.1</td><td>2.72</td><td>11.87</td><td>16.61</td><td>149.47</td><td>0.77</td></tr><tr><td>26</td><td>2.34</td><td>1.25</td><td>18.29</td><td>2.68</td><td>1.8</td><td>14.67</td><td>14.26</td><td>147.18</td><td>1.33</td></tr></tbody></table>

H<sub>tree</sub> represents the minimal tree height, and Q<sub>stem</sub> represents the maximal biomass in the Pareto optimal front for different strategies.

is the plant response to biomechanical signals. Such knowledge is applied in a greenhouse by posing artificial wind to have stronger plants. Under lower light, the plant tends to be thinner (with lower biomass) yet keeps tall to catch the light (with modified allometry, i.e., the ratio between length and diameter for the internode). This is a usually unexpected phenotype as plants will be prone to breakage. Blowing the plant with wind can compensate for plants grown in a greenhouse, where light is less intensive than in an open field. An integrated modeling system with light and wind can support quantitative evaluation. However, most modeling studies on plant plasticity concentrate on the response of plants to other environmental conditions like temperature (Liu et al., 2016), light (Kang et al., 2012), or irrigation (Wu et al., 2012), without mechanical constraint. Combining the effect of both light and wind is feasible as light environment simulation has become very common in the FSPM community (Wang et al., 2012). The balance between light and wind can be simulated: to have a better light interception, the leaves need to be more dispersed, while the biomechanical constraint may require a more compact distribution of leaves. An example can be found in Eloy et al. (2017) who explored tree fitness under both light and biomechanical signals. Light competition influences plant height through both the whole tree biomass and architecture, as simulated in Cournède et al. (2008) as well as the allometries; in some cases, stronger light competition leads to higher individual plants (Chen et al., 2010). The current paper deals with trees grown in the open field where the light condition is uniform, and it is not expressed explicitly in the environmental variable E(n) (Eq. S15). Its competition level is expressed with the ground projection area S<sub>p</sub> (Eq. S15) as described in Cournède et al. (2008). For the indoor environment as in the greenhouse, the light effect can be considered in the

integrated variable E(n) as done in Fan et al. (2015). The sophisticated interplay between light and wind can be better understood by introducing a realistic light component.

Tree reaction to a biomechanical stimulus can be split into two complementary processes: (1) the reallocation of biomass within the whole plant and stem allometry changing and (2) the differentiation of wood cells (e.g., formation of reaction wood) and active biomechanical control of stem shape due to wood fiber maturation strains (Moulia et al., 2006). This paper concentrates on the first complementary mechanism by finding the parameter controlling the reallocation of biomass to secondary growth. It may be called a goal-seeking approach (Takahara & Mesarovic, 2003) in optimizing dry matter distribution with mechanical constraints: instead of simulating directly the impact of wind on plant growth, it searches for suitable growth parameters that lead to optimal biomass and height at given wind load. A similar idea of simulating plant acclimation through optimization can be found in Eloy et al. (2017), with a conceptual model. In the current work, the reaction of the tree is reflected by a single parameter S<sub>layer</sub>. Even with a constant S<sub>layer</sub>, the allocation changes because of tree organogenesis (Figure 10E). The results fit the situation where trees receive stable wind. When trees grow in varying wind environments, a variable S<sub>layer</sub> is expected. Such variable value has been found on maple trees (De Reffye et al., 2017).

The numerical results show that the response of the plant to the wind is mild for a breeze. The response takes place gradually according to Figure 10. It is applicable for the situation where plants receive frequent wind loads for a long period. However, when the tree starts to respond to the wind load is not studied here. How plants react instantly to critical wind load that causes top stem deflection is uninvestigated. Due to the limit of
<figure><img src="image_13.png" /><figcaption>FIGURE 10
Simulated tree morphological structures reacting to wind with strategies I, II, and III. (A) Plant biomass (kg); (B) plant height (m); (C) trunk diameter at breast height (cm); (D) tree leaf area (m²); (E) biomass partitioning to secondary growth (%); and (F) ratio for trunk biomass to the sub-branch biomass (%).</figcaption></figure>

experiment conditions, most experiments were done for low wind speeds. The *in silico* experiment gives simulated results under strong wind that is hard to achieve and needs further verification. For trees, through measurement and model calibration, all sink-source parameters can be estimated, which is the basis of the *in silico* experiment. In this work, the model parameters are based on both literature and previous study. The role of wind is not limited to a biomechanical stimulus; it also influences the photosynthesis process by modifying leaf boundary layer conductance, which in turn influences leaf surface CO₂ concentration and humidity. As wind speed increases, the increasing transpiration rate can lead to a

decrease in the water content of leaves (Yabuki, 2004). The simulations of such physiological processes are expected to bring more realistic results.

## 5 Conclusion

The compensating mechanism of plants in reacting to the wind environment is simulated using a system coupling a FSPM with a biomechanical model. A smaller tree (with smaller plant height and leaf area) with higher wind resistance is obtained at higher wind speed, as an emergent property of the system.
<figure><img src="image_14.png" /><figcaption>FIGURE 11</figcaption></figure>

Simulated optimal tree plasticity for strategy II under different wind speeds. (A) 0.0 m/s; (B) 8 m/s; (C) 16 m/s; (D) 20 m/s (E) 22 m/s. Images are created under calm conditions.

Secondary growth increases to avoid trunk breakage. This work allows investigating the compensating mechanism of biomass computationally and increases the potential of FSPM for a wider application.

## Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

## Author contributions

MK, TF, and PDR planned and designed the research. HW performed the experiments and analyzed the data. MK, HW, JH, XW, XF and TF wrote and revised the manuscript. All authors contributed to the article and approved the submitted version.

## Funding

This work is supported by the Major S&T project (Innovation 2030) of China (2021ZD0113704), the National Natural Science Foundation of China (62076239), the CAS-NSTDA Joint Research Program (GJHZ2076), the Natural Science Startup Foundation of Chongqing Technology and Business University (2056019) and the Science and Technology Research Program of Chongqing Municipal Education Commission (KJQN201900833).

## Acknowledgments

The authors thank LIAMA for its support within the cPlant project framework.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.971690/full#supplementary-material
## References

Almeras, T., Gril, J., and Costes, E. (2002). Bending of apricot tree branches under the weight of axillary growth: test of a mechanical model with experimental data. *Trees* 16, 5–15. doi: 10.1007/s00468-001-0139-1

Ancelin, P., Courbaud, B. I. T., and Fourcaud, T. (2004a). Development of an individual tree-based mechanical model to predict wind damage within forest stands. *For. Ecol. Manage.* 203, 101–121. doi: 10.1016/j.foreco.2004.07.067

Ancelin, P., Fourcaud, T., and Lac, P. (2004b). Modelling the biomechanical behaviour of growing trees at the forest stand scale. part I: Development of an incremental transfer matrix method and application to simplified tree structures. *Ann. For. Sci.* 61 263–275. doi: 10.1051/forest:2004019

Biscani, F., and Izzo, D. (2020). A parallel global multiobjective framework for optimization: pagmo. *J. Open Sour. Softw.* 5, 2338. doi: 10.21105/joss.02338

Bonesneueur, V., Constant, T., Moulia, B., and Fournier, M. (2016). Forest trees filter chronic wind-signals to acclimate to high winds. *New Phytol.* 210, 850–860. doi: 10.1111/nph.13836

Calvo, P., Gagliano, M., Souza, G. M., and Trewavas, A. (2019). Plants are intelligent, here's how. *Ann. Bot.* 125, 11–28. doi: 10.1093/aob/mcz155

Chen, S., Zhang, J., Jia, P., Xu, J., Wang, G., and Xiao, S. (2010). Effects of size variation and spatial structure on plastic response of plant height to light competition. *Chin. Sci. Bull.* 55, 1135–1141. doi: 10.1007/s11434-010-0107-5

Collet, C., Lanter, O., and Pardos, M. (2001). Effects of canopy opening on height and diameter growth in naturally regenerated beech seedlings. *Ann. For. Sci.* 58, 127–134. doi: 10.1051/forest:2001112

Cournède, P.-H., Mathieu, A., Houllier, F., Barthélémy, D., and de Reffye, P. (2008). Computing competition for light in the GREENLAB model of plant growth: A contribution to the study of the effects of density on resource acquisition and architectural development. *Ann. Bot.* 101, 1207–1219. doi: 10.1093/aob/mcm272

Coutand, C., Mathias, J. D., Jeronimidis, G., and Destrebecq, J. F. O. (2011). TWIG: A model to simulate the gravitropic response of a tree axis in the frame of elasticity and viscoelasticity, at intra-annual time scale. *J. Theor. Biol.* 273, 115–129. doi: 10.1016/j.jtbi.2010.12.027

Coutand, C., and Moulia, B. (2000). Biomechanical study of the effect of a controlled bending on tomato stem elongation: local strain sensing and spatial integration of the signal. *J. Exp. Bot.* 352, 1825–1842. doi: 10.1093/jexbot/51.352.1825

Coutand, C., Pot, G., and Badel, E. (2014). Mechanosensing is involved in the regulation of autostress levels in tension wood. *Trees* 28, 687–697. doi: 10.1007/s00468-014-0981-6

Davenport, A. G. (1962). The spectrum of horizontal gustiness near the ground in high winds. *J. R. Meteorol. Soc* 88, 197–198. doi: 10.1002/qj.49708837618

Deb, K., Agrawal, S., Pratap, A., and Meyarivan, T. (2000). "A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II," in *Parallel problem solving from nature* PPSN VI. Eds. M. Schoenauer, K. Deb, G. Rudolph, X. Yao, E. Lutton, J. J. Merelo and H.-P. Schwefel (Berlin, Heidelberg: Springer Berlin Heidelberg), 849–858.

de Langre, E. (2008). Effects of wind on plants. *Annu. Rev. Fluid Mech.* 40, 141–168. doi: 10.1146/annurev.fluid.40.111406.102135

de Reffye, P. (1976). Modélisation et simulation de la verse du caféier, l'aide de la théorie de la résistance des matériaux. *Cafè Cacao Thè* 20, 251–272. Available at: https://agritrop.cirad.fr/454749/

de Reffye, P., Hu, B.-G., Kang, M., Letort, V., and Jaeger, M. (2020). Two decades of research with the GreenLab model in agronomy. *Ann. Bot.* 127, 281–295. doi: 10.1093/aob/mcaa172

De Reffye, P., Jaeger, M., Barthélémy, D., and Houllier, F. (2017). *Architecture des plantes et production végétale : Les apports de la modélisation mathématique* (Versailles: Quae Edition).

Diener, J., Rodriguez, M., Baboud, L., and Reveret, L. (2010). Wind projection basis for real-time animation of trees. *Comput. Graphics Forum* 28, 533–540. doi: 10.1111/j.1467-8659.2009.01393.x

Eloy, C., Fournier, M., Lacointe, A., and Moulia, B. (2017). Wind loads and competition for light sculpt trees into self-similar structures. *Nat. Commun.* 8, 1014. doi: 10.1038/s41467-017-00995-6

Fan, X.-R., Kang, M.-Z., Heuvelink, E., de Reffye, P., and Hu, B.-G. (2015). A knowledge-and-data-driven modeling approach for simulating plant growth: A case study on tomato growth. *Ecol. Model.* 312, 363–373. doi: 10.1016/j.ecolmodel.2015.06.006

Fjeld, L. (2012). Modeling MOR in Norway spruce (Picea abies (L.) karst.) structural lumber with stands and tree characteristics. *Holz Roh Werkst* 66, 219–227.

Fourcaud, T., Blaise, F., Lac, P., Cast E Ra, P., and de Reffye, P. (2003). Numerical modelling of shape regulation and growth stresses in trees. II. implementation in

the AMAPpara software and simulation of tree growth. *Trees* 17, 31–39. doi: 10.1007/s00468-002-0203-5

Fourcaud, T., and Lac, P. (2003). Numerical modelling of shape regulation and growth stresses in trees i. an incremental static finite element formulation. *Trees* 17, 23–30. doi: 10.1007/s00468-002-0202-6

Fourcaud, T., Zhang, X., Stokes, A., Lambers, H., and Körner, C. (2008). Plant growth modelling and applications: The increasing importance of plant architecture in growth models. *Ann. Bot.* 101, 1053–1063. doi: 10.1093/aob/mcm050

Fournier, M., Bailleres, H., and Chanson, B. (1994). Tree biomechanics: Growth, cumulative prestresses, and reorientations. *Biomimetics* 2, 229–251. Available at: https://agris.fao.org/agris-search/search.do?recordID=US9506632

Gardiner, B., Peltola, H., and Kellomäki, S. (2000). Comparison of two models for predicting the critical wind speeds required to damage coniferous trees. *Ecol. Model.* 129, 1–23. doi: 10.1016/S0304-3800(00)00220-9

Grotta, A. T., Leichti, R. J., Gartner, B. L., and Johnson, G. R. (2005). Effect of growth ring orientation and placement of earlywood and latewood on moe and mor of very-small clear Douglas-fir beams. *Wood Fiber.* 37, 207–212.

Guillon, T., Dumont, Y., and Fourcaud, T. (2012). A new mathematical framework for modelling the biomechanics of growing trees with rod theory. *Math. Comp. Model.* 55, 2061–2077. doi: 10.1016/j.mcm.2011.12.024

Guo, H., Lei, X., Courn E De, P.-H., and Letort, V. (2012). Characterization of the effects of inter-tree competition on source-sink balance in Chinese pine trees with the GreenLab model. *Trees* 26, 1057–1067. doi: 10.1007/s00468-012-0683-x

Hamant, O. (2013). Widespread mechanosensing controls the structure behind the architecture in plants. *Curr. Opin. Plant Biol.* 16, 654–660. doi: 10.1016/j.pbi.2013.06.006

Irvine, M. R., Gardiner, B. A., and Hill, M. K. (1997). The evolution of turbulence across a forest edge. *Boundary-Layer Meteorol.* 84, 467–496. doi: 10.1023/A:1000453031036

Jaffe, M. J. (1973). Thigmomorphogenesis: the response of plant growth and development to mechanical stimulation. *Planta* 114, 143–157. doi: 10.1007/BF00387472

Jaffe, M. J., and Forbes, S. (1993). Thigmomorphogenesis: the effect of mechanical perturbation on plants. *Plant Growth Regul.* 12, 313–324. doi: 10.1007/BF00027123

Jones, H. G. (1983). *Plants and microclimate. a quantitative approach to environmental plant physiology* (Cambridge: Cambridge University Press).

Kang, M. Z., Evers, J. B., Vos, J., and De Reffye, P. (2008). The derivation of sink functions of wheat organs using the GREENLAB model. *Ann. Bot.* 101, 1099–1108. doi: 10.1093/aob/mcm212

Kang, M. Z., Fan, X.-R., Hua, J., Wang, H. Y., Wang, X. J., and Wang, F.-Y. (2018). Managing traditional solar greenhouse with CPSS: A just-for-fit philosophy. *IEEE T. Cybernet.* 48, 3371–3380. doi: 10.1109/TCYB.2018.2858264

Kang, M. Z., Heuvelink, E., Carvalho, S. M. P., and de Reffye, P. (2012). A virtual plant that responds to the environment like a real one: the case for chrysanthemum. *New Phytol.* 195, 384–395. doi: 10.1111/j.1469-8137.2012.04177.x

Koizumi, A., J-i, M., Sawata, K., Sasaki, Y., and Hirai, T. (2010). Evaluation of drag coefficients of poplar-tree crowns by a field test method. *J. Wood Sci.* 56, 189–193. doi: 10.1007/s10086-009-1091-8

Kord, B., Kialashaki, A., and Kord, B. (2010). The within-tree variation in wood density and shrinkage, and their relationship in populus euramericana. *Turk. J. Agric. For.* 34, 121–126. doi: 10.3906/tar-0903-14

Letort, V., Cournède, P. H., Mathieu, A., de Reffye, P., and Constant, T. (2008). Parametric identification of a functional-structural tree growth model and application to beech trees fagus sylvatica. *Funct. Plant Biol.* 35, 951–963. doi: 10.1071/FP08065

Liu, B., Asseng, S., Liu, L., Tang, L., Cao, W., and Zhu, Y. (2016). Testing the responses of four wheat crop models to heat stress at anthesis and grain filling. *Global Change Biol.* 22, 1890–1903. doi: 10.1111/gcb.13212

Lopez, D., Michelin, S. E. B., and de Langre, E. (2011). Flow-induced pruning of branched systems and brittle reconfiguration. *J. Theor. Biol.* 284, 117–124. doi: 10.1016/j.jtbi.2011.06.027

Mayhead, G. J. (1973). Some drag coefficients for British forest trees derived from wind tunnel studies. *Agric. Meteorol.* 12, 123–130. doi: 10.1016/0002-1571 (73)90013-7

Moulia, B., Coutand, C., and Lenne, C. (2006). Posture control and skeletal mechanical acclimation in terrestrial plants: Implications for mechanical modeling of plant architecture. *Am. J. Bot.* 93, 1477–1489. doi: 10.3732/ajb.93.10.1477
Niklas, K. J., and Spatz, H.-C. (2000). Wind-induced stresses in cherry trees: evidence against the hypothesis of constant stress levels. *Trees* 14, 230–237. doi: 10.1007/s004680050008

Peltola, H., and Kellomäki, S. (1993). A mechanistic model for calculating windthrow and stem breakage of scots pines at stand age. *Silva Fenn.* 27, 5504. doi: 10.14214/sf.a15665

Perttunen, J., Sievänen, R., Nikinmaa, E., Salminen, H., Saarenmaa, H., and Väkevä, J. (1996). LIGNUM: A tree model based on simple structural units. *Ann. Bot.* 77, 87–98. doi: 10.1006/anbo.1996.0011

Qi, R., Letort, V., Kang, M., Cournede, P.-H., de Reffye, P., and Fourcaud, T. (2009). Application of the GreenLab model to simulate and optimize wood production and tree stability: a theoretical study. *Silva Fenn.* 43, 465–487. doi: 10.14214/sf.201

Sellier, D., Brunet, Y., and Fourcaud, T. (2008). A numerical model of tree aerodynamic response to a turbulent airflow. *Forestry* 81, 279–297. doi: 10.1093/forestry/cpn024

Takahara, Y., and Mesarovic, M. (2003). "Goal seeking system," in *Organization structure*. Eds. Y. Takahara and M. Mesarovic (Boston, MA: Springer US), 15–35.

Taylor-Hell, J., Lawrence, D., Harder, D., and Faramarz Samavati, D. (2005).
*Biomechanics in botanical trees* (Master’s thesis, The University of Calgary).

Telewski, F. W., and Pruyn, M. L. (1998). Thigmomorphogenesis: a dose response to flexing in ulmus americana seedlings. *Tree Physiol.* 18, 65–68. doi: 10.1093/treephys/18.1.65

Trewavas, A. (2017). The foundations of plant intelligence. *Interface Focus* 7, 20160098. doi: 10.1098/rsfs.2016.0098

Wang, H. Y., Kang, M. Z., and Hua, J. (2012). "Simulating plant plasticity under light environment: A source-sink approach," in *2012 IEEE 4th International Symposium on Plant Growth Modeling, Simulation, Visualization and Applications*. 431–438. doi: 10.1109/PMA.2012.6524869

Whitehead, F. H. (1962). Experimental studies of the effect of wind on plant growth and anatomy. II. helianthus annuus. *New Phytol.* 61, 59–62. doi: 10.1111/j.1469-8137.1962.tb06274.x

Whitehead, F. H., and Luti, R. (1962). Experimental studies of the effect of wind on plant growth and anatomy. i. zea mays. *New Phytol.* 61, 56–58. doi: 10.1111/j.1469-8137.1962.tb06273.x

Wu, L., Le Dimet, F.-X., de Reffye, P., Hug, B.-G., Cournede, P.-H., and Kang, M. Z. (2012). An optimal control methodology for plant growth-case study of a water supply problem of sunflower. *Math. Comput. Simul.* 82, 909–923. doi: 10.1016/j.matcom.2011.12.007

Yabuki, K. (2004). *Photosynthetic rate and dynamic environment* (Dordrecht: Springer).

Yang, G., Cao, W. Q., Liu, D. X., Kang, M. Z., Letort, V., Xing, M. J., et al. (2011). Growth and development simulation based on functional-structural model GreenLab for poplar (Salicaceae). *Afr. J. Agric. Res.* 6, 3071–3077. doi: 10.5897/AJAR10.855