# Application of the GreenLab model to simulate and optimize wood production and tree stability: a theoretical study

Rui Qi, Véronique Letort, M. Z. Kang, Paul-Henry Cournède, Philippe de Reffye, Thierry Fourcaud

► To cite this version:

Rui Qi, Véronique Letort, M. Z. Kang, Paul-Henry Cournède, Philippe de Reffye, et al.. Application of the GreenLab model to simulate and optimize wood production and tree stability: a theoretical study. Silva Fennica, 2009, 43 (3), pp.465-487. halsde-00436297

HAL Id: halsde-00436297

https://hal.science/halsde-00436297

Submitted on 30 May 2020

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>
# Application of the GreenLab Model to Simulate and Optimize Wood Production and Tree Stability: A Theoretical Study

Rui Qi¹,², Véronique Letort¹, Mengzhen Kang², Paul-Henry Cournède ¹,³,
Philippe de Reffye ³,⁴ and Thierry Fourcaud⁴

Qi, R., Letort, V., Kang, M., Cournède, P-H., de Reffye, P. & Fourcaud, T. 2009. Application of the GreenLab model to simulate and optimize wood production and tree stability: a theoretical study. Silva Fennica 43(3): 465–487.

The GreenLab model was used to study the interaction between source-sink dynamics at the whole tree level, wood production and distribution within the stem, and tree mechanical stability through simulation and optimization. In this first promising numerical attempt, two GreenLab parameters were considered in order to maximize wood production: the sink strength for cambial growth and a coefficient that determines the way the biomass assigned to cambial growth is allocated to each metamer, through optimization and simulation respectively. The optimization procedure that has been used is based on a heuristic optimization algorithm called Particle Swarm Optimization (PSO). In the first part of the paper, wood production was maximized without considering the effect of wood distribution on tree mechanical stability. Contrary to common idea that increasing sink strength for cambial growth leads to increasing wood production, an optimal value can be found. The optimization results implied that an optimal source and sink balance should be considered to optimize wood production. In a further step, the mechanical stability of trees submitted to their self weight was taken into account based on simplified mechanical assumptions. Simulation results revealed that the allocation of wood at the stem base strongly influenced its global deformation. Such basic mechanical criterion can be an indicator of wood quality if we consider further the active biomechanical processes involved in tree gravitropic responses, e.g. formation of reaction wood.

**Keywords** biomechanics, FSPM, Particle Swarm Optimization, source-sink dynamics, optimization, wood quality, biomass allocation

**Addresses** ¹Ecole Centrale Paris, Laboratory of Applied Mathematics, Grande Voie des Vignes, 92295 Chatenay-Malabry, France; ²Institute of Automation, Chinese Academy of Sciences, LIAMA/NLPR, P.O.Box 2728, Beijing, China, 100190; ³INRIA saclay Ile-de-France, EPI Digiplant, Parc Orsay Université, 91893 Orsay cedex, France; ⁴CIRAD, UMR AMAP, Montpellier, F-34000 France **E-mail** qiruitree@gmail.com

**Received** 2 June 2008 **Revised** 12 February 2009 **Accepted** 14 May 2009

**Available at** http://www.metla.fi/silvafennica/full/sf43/sf433465.pdf
# 1 Introduction

Nowadays, the economical importance of wood production is enhanced by an increasing demand from South-East Asia (Zhang et al. 2006), by increasing custom duties on wood export from Russia (Russia to go... 2007) and by extending fields of applications: paper, pulp, joinery, furniture, biofuel, etc (Nepveu 1993). In terms of environmental protection and resource utilization, wood as a renewable source attracts more and more attention (Wang et al. 2004, Forestry commission in England). For most application fields, not only wood production but also wood quality is important. These two criteria are of primordial importance to determine the economical value of logs.

Woodcock and Shier (2003) observed that the factors influencing wood specific gravity were tree diameter, height and growth rate. Furthermore, Wiemann and Williamson (1989ab) pointed out the correlation between the wood specific gravity and the structural Young's modulus that defines wood rigidity. Wood quality is thus closely linked with the need in mechanical stability of trees through both structural and material effects.

Tree response to mechanical forces, e.g. self weight and wind loads, depends on tree structure and wood mechanical properties, which can be both modified by insects and diseases or silvicultural practices. In addition to external loads, trees can develop internal biomechanical stresses due to wood cells maturation (Fournier et al. 2006). This active response aims at controlling the shape of tree axes and their position in space. In the particular case of trunks, this active mechanism usually corresponds to a negative gravitropism that allows trees to remain vertical. Even though these biomechanical phenomena have been widely studied in the past, the underlying biological mechanisms, e.g. mechanoperation and signal transduction, are still poorly known (Moulia et al. 2006). Modelling and simulating tree biomechanics is a helpful approach that can be used to support or test biological hypotheses (Fourcaud et al. 2003). Nevertheless such biomechanical models are not classical as the mechanical solicitations apply on a growing structure (Ancelin et al. 2004a, Fourcaud and Lac 2003, Fourcaud et al. 2003). However, the goal of this study was not to

consider these complex biomechanical processes, but to introduce a simple mechanical criterion of stem stability in a functional-structural plant growth model (FSPM) in order to assess the influence of model parameters on wood production and wood quality.

Wood quality can be defined with regard to several criteria including for instance ring and stem profiles, knot distribution and sizes, wood mechanical and chemical properties, heterogeneity of these properties within the stem and rings, stem straightness and/or cross section shapes. So far, most of the existing wood quality models include one or more of these properties (Nepveu 1993, Denne et al. 1994, Houllier et al. 1995, Deleuze and Houllier 1997, Woodcock and Shier 2003, Lasserre et al. 2005). Jayawickrama (2001) pointed out that among these wood properties, wood stiffness and stability are of major importance in fast growth radiate pines. Moreover, Lasserre et al. (2005) found that there is a negative relationship between stand density and wood stiffness (an increase in stand density being associated to a decrease in stem DBH). Therefore, in this paper, wood quality was evaluated according to tree mechanical stability, which involved the coupling influences between tree structure, load distribution and wood mechanical properties (Ancelin et al. 2004a). For this purpose, a mechanical stability criterion has been used, which was based on stem buckling equations. We assume that this criterion summarizes both wood mechanical properties and the mechanical state of the tree, excluding the effects due to biomechanical responses as defined above, i.e. formation of reaction wood and induction of growth stresses. The final displacement angle of the top stem was considered as an indicator of tree stability.

Wood production and wood quality are the results of processes involving both topological development and physiological processes of plants. Studying these process interactions through simulation requires specific properties for plant growth models. Plant growth models can be classified into four categories: geometrical models, empirical models, process-based models and functional-structural models. Geometrical models (Greene 1989, Deussen et al. 1998, Dornbusch et al. 2007) aim at generating plant shapes looking realistic for their morphology.
Organ sizes are determined by computer graphic techniques, without considering any physiological information. Empirical models (Peng 2000, Fox et al. 2008, Repola 2008) use linear or nonlinear regression techniques to establish the relationship among some particular features of plants such as diameter at breast height and height. They provide ways to analyse plant growth by extracting statistical correlations between variables of interest, however, they are not based on a mechanistic approach. Hence, their range of application is often limited. On the contrary, process-based models offer a mechanistic representation of the plant physiological mechanisms (Mäkelä et al. 2000, Peng 2000), yet the plant is represented as a set of compartments, i.e. compartment weight of each type of organs at a certain stage. The topological description is so simplified that detailed information of each organ at different positions (weight, dimension) is impossible to obtain. For this work on tree mechanical stability, the final displacement angle is related to morphological characteristics, such as stem diameter and variation of ring width (Nepveu 1993, Denne et al. 1994), and tree height (Deleuze and Houllier 1997, Woodcock and Shier 2003). These properties are influenced by growth physiological processes (Sievänen et al. 2000, Woodcock and Shier 2003, Mathieu et al. 2008). As a result, the widely used process-based plant growth models are not sufficient to model interactively the plant functioning and structural development. This is the role of functional-structural models (Sievänen et al. 2000, Prusinkiewicz 2004, Fourcaud et al. 2008), such as AMAPpara (de Reffye et al. 1997), L-Peach (Allen et al. 2005), LIGNUM (Perttunen et al. 1998), and GreenLab (Yan et al. 2004, Cournède et al. 2006).

In this study, the functional-structural plant growth model GreenLab is chosen to simulate wood production and tree mechanical stability. Its ability to describe the interaction between plant architecture and plant physiological processes (Mathieu et al. 2008) is important to carry out studies on wood production and wood quality with regard to plant development. The plant topological structure is generated using predefined rules, which can be formulated by a single equation (Yan et al. 2004, Cournède et al. 2006). Biomass production and allocation is modelled using

a source-sink approach. Several methods have been considered to represent the plant organogenesis depending on the tree species. One of the methods is to consider a stochastic development (Kang et al. 2008), where the number of metamers appearing at each time depends on certain probabilities. A more mechanistic technique to simulate the indeterminate topological structure of plants is to correlate it with physiological processes (Mathieu 2006), i.e. whether the metamers emerge depends on the ratio of biomass supply to demand (Mathieu et al. 2009). The main goal of this work is to present an optimization procedure for wood production and to study the influence on tree mechanical stability, which was considered as a criterion to assess wood quality, based on a functional-structural growth model. To our knowledge, such a numerical approach is novel. Therefore we choose in this preliminary study to consider a deterministic plant topological development. However, let note that although the number of metamers appearing at each time is predefined, the sizes of organs depend on the plant biomass production that is determined by the sizes of leaves (leaf surface area) through photosynthesis.

The GreenLab model that has been used is a generic model, not devoted to a particular species. Carson et al. (2006) pointed out that due to the complex system of plant growth, the challenge on plant growth models for parameter estimation and model validation must be overcome for model applications. The parameters of GreenLab have been well estimated from real data for kinds of species. Besides being applied on a dozen of crop species (de Reffye et al. 2008), it has been applied to several tree species including chinese pine tree (Guo et al. 2006), as well as beech tree (Letort et al. 2008a). In all cases, GreenLab proved to be an efficient tool to describe the source-sink dynamics in plant growth. Furthermore, the mathematical formulation of GreenLab under the form of a discrete dynamic system enables the efficient use of mathematical algorithms for solving optimization problems. An application is illustrated in the present study where the total stem wood production was maximized considering the sink strength for cambial growth as an optimization variable in GreenLab. For this purpose we used a heuristic optimization algorithm, namely the
Particle Swarm Optimization (PSO) (He et al. 2004). In a second step, stem straightness and stem profile were investigated numerically with regard to another GreenLab parameter that defines the distribution of biomass to each metamer along the stem. Initial simulation analyses showed that the trunk tends to be straighter, i.e. more stable, with an increase in the parameter. These previous results indicated that the use of the complex optimization algorithm to optimize stem stability with regard to the parameter is not required.

# 2 Methods

## 2.1 The GreenLab Model

GreenLab is a functional-structural, organ-based plant growth model. Organs are identified by two properties: the *physiological age* that is a botanical variable characterizing the morphological differentiation of organs (Barthélémy et al. 1997) and the *chronological age* that is equivalent to the duration from the emergence of organs. Organs (leaves, internodes, fruits if they exist) with the same physiological and chronological ages have the same physiological behaviour. GreenLab simulates the response of individual plants to their local environmental conditions. The model takes into account the competition between individuals at a given population density (Cournède et al. 2008). From a mathematical point of view, GreenLab is a discrete dynamic system. The choices of spatial and temporal discretization scales relied on botanical observations (Barthélémy and Caraglio 2007): spatial units are metamers and the model time steps are growth cycles. The growth cycle duration is the time span needed to generate a new growth unit containing one or several metamers, e.g. for trees in temperate zones, the growth cycle is generally one year. Organogenesis (production of new metamers) is modelled according to pre-defined rules using a dual-scale automaton (Zhao et al. 2003). Biomass production and allocation are based on source-sink relationships and calculated using recurrence mathematical equations within a growth cycle, which is synchronous with organogenesis mechanism. The biomass produced through photosynthesis is supposed to be first

stored in a common pool at the beginning of each growth cycle, and then allocated to each organ, including cambial growth. Emergence and elongation of organs (e.g. length increase of metamers) occur during a phase of primary growth, and cambial growth along the stem and branches (e.g. diameter increase of the stem and branches) is called secondary growth. The principles of GreenLab were extensively described in detail in several papers (de Reffye and Hu 2003, Yan et al. 2004, Cournède et al. 2006, Letort et al. 2008a). Here, we recall its basic principles and mostly focus on the specificities used in this work.

### 2.1.1 Organogenesis: Production of Metamers and Growth Units

In the tree development, metamers are characterized by two indices: physiological age of the organs and that of their potential buds. The metamers belonging to a same growth unit have the same physiological age of the organs, but their axillary buds can have different physiological ages. In this study, we will not consider the case of reiteration, i.e. the case where the physiological age of the axillary bud of a metamer is the same as the one of its main axis. The emergence sequence of metamers follows a predefined rule, modelled by a dual-scale automaton (Zhao et al. 2003). In Fig. 1, the different potential states of the automaton are represented with the rules defining their succession. These rules determine the number of growth cycles before transitions from current metamer or growth unit to the following ones (represented by arrows). The corresponding skeleton of a two-year-old tree and the corresponding topological structure of a 40-year-old tree are shown at the right-bottom of Fig. 1. The topology of the tree in this work is generated by the automaton illustrated in Fig. 1. In addition to axillary buds, a metamer can possibly bear leaves and fruits. According to the dual-scale automaton, the number of each kind of metamers at any time can be calculated. Hence, the number of each organ type is known. Owing to the use of substructure splitting and factorization techniques, the computing performance of growth simulations was significantly improved (Cournède et al. 2006).
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Presentation of GreenLab organogenesis mechanism based on a dual-scale automaton: the different physiological ages (PAs) of buds (or metemers) are distinguished by different colours; circles represent buds; rectangles represent metamers; the number above each small ellipse represents the number of metamers with the corresponding PA inside the growth cycle; the number above each big ellipse represents the number of growth units with the corresponding PA in the whole tree; The arrow represents the transition from one metamer or growth unit to another one. At the right-bottom of the figure, according to the organogenesis mechanism, the skeleton of a two-year-old tree and the topological structure of a 40-year-old tree without leaves are shown, where dead branches are not drawn.</figcaption></figure>

### 2.1.2 Biomass Production and Allocation to Primary Growth

Biomass production and allocation processes in GreenLab are modelled using a source-sink approach. The biomass increment of an individual plant at growth cycle $n$, denoted by $Q(n)$, is calculated as Eq. 1, which describes the light interception by foliage.

$$Q(n) = E(n)Sp \left( 1 - \exp \left( -k \cdot \frac{S(n)}{Sp} \right) \right) \quad (1)$$

$$Q(0) = Q_{\text{seed}}$$

where $Q_{\text{seed}}$ is the seed biomass; $E(n)$ is a variable representing the plant local environment at growth cycle $n$, which includes the incident photosynthetically active radiation, radiation used efficiency and the hydraulic resistance to transpiration; $Sp$ is the total ground projection area available of the crown for plant modulated by the effects of self-shading and neighbour competition

that is related to plant density; $k$ is a light extinction coefficient; $S(n)$ is the total green leaf surface area at growth cycle $n$; hence, the ratio of $S(n)$ to $Sp$ represents the leaf area index (LAI) adapted to individual plants.

For trees in temperate zones, organs finish elongation within one year that corresponds to one growth cycle in GreenLab. Therefore, the quantity of biomass allocated to an organ $o$ of physiological age $p$ when the plant age is $n$, denoted by $q_p^o(n)$, is calculated as follows:

$$q_p^o(n) = P_p^o \cdot \frac{Q(n)}{D(n)} \quad (2)$$

where subscript $o$ represents a leaf (a) or a pith (i); $P_p^o$ is the sink strength of organ $o$ of physiological age $p$; $D(n)$ is the total plant demand at growth cycle $n$, which is the sum of the sinks of all organs that need elongation of the plant (see more details in a following section).

According to Eq. 2, the total green leaf weight of the plant at growth cycle $n$, denoted by $QB(n)$,
is given by Eq. 3,

$$QB(n) = \sum_{p=1}^{P_m} \sum_{k=1}^{ta} \left[ N_p^a(n-k+1) \cdot P_p^a \cdot \frac{Q(n-k+1)}{D(n-k+1)} \right] \quad (3)$$

where $N_p^a(n)$ is the number of leaves of physiological age $p$, initiated at growth cycle $n$, which is derived from the organogenesis mechanism;

$$Q(n) = E(n)Sp \left( 1 - \exp \left( - \frac{k}{Sp \cdot slw} \sum_{p=1}^{P_m} N_p^a(n-1) \cdot P_p^a \cdot \frac{Q(n-1)}{D(n-1)} \right) \right) \quad (4)$$

Eq. 4 illustrates the interaction between organogenesis (defining the evolution of organ number) and physiological processes in GreenLab.

### 2.1.3 Biomass Allocation to Secondary Growth

At each growth cycle, a new ring is formed along the stem or branches. In GreenLab, this phenomenon is modelled in two steps. First, the total amount of biomass allocated to cambial growth at growth cycle $n$, denoted by $Q_{ring}(n)$, is calculated at the whole-plant scale, as Eqs. 5–7.

$$Q_{ring}(n) = \frac{D_{ring}(n)}{D(n)} Q(n) \quad (5)$$

with

$$D_{ring}(n) = S_{ring} \sum_{p=1}^{P_m} N_p^a(n-1) \quad (6)$$

and

$$D(n) = \sum_{p=1}^{P_m} \left[ N_p^a(n) \cdot P_p^a + N_p^i(n) \cdot P_p^i \right] + D_{ring}(n) \quad (7)$$

where $D_{ring}(n)$ is the total demand for cambial growth at growth cycle $n$, which is proportional to the total number of leaves alive at the previous growth cycle; $S_{ring}$ is the sink strength for cambial growth; the total plant demand $D(n)$ corresponds to demands of new organs of all physiological ages and to the demand of the cambial growth compartment.

This total amount of biomass for cambial growth $Q_{ring}(n)$ is then allocated to each metamer depending on its topological position. According

$ta$ is the leaf functioning duration, i.e. the time during which leaves are photosynthetically active; $P_m$ is the maximum physiological age.

The leaf surface area is equal to the leaf weight divided by the specific leaf weight (slw). Considering that the leaf functioning duration $ta$ is equal to one growth cycle for the tree presented in this work (deciduous trees), the increment of plant biomass at growth cycle $n$ becomes

to the Pressler law and to the pipe model theory (Shinozaki et al. 1964), biomass increment for cambial growth of a metamer depends on the number of active leaves above it in the plant architecture. However, some restrictions to this rule were observed (Deleuze and Houllier 2002). To overcome these limitations, a mixed approach was proposed, which allowed allocation of the biomass for cambial growth both along upward and downward pathways, using a coefficient $\lambda$. In the first mode ($D_1$), metamer demand for cambial growth is determined as a mere sink sub-model, while in the second mode ($D_2$) the number of active leaves is taken into account, as Eqs. 8–9.

$$D_1(n) = \sum_{p=1}^{P_m} \sum_{j=1}^{n} N_p^i(n-j) \cdot l_p(n-j) \quad (8)$$

$$D_2(n) = \sum_{p=1}^{P_m} \sum_{j=1}^{n} NL_p(n-j) \cdot l_p(n-j) \quad (9)$$

where $N_p^i(n-j)$ is the number of metamers (piths) of physiological age $p$, chronological age $j$, at growth cycle $n$, which initiated at growth cycle $n-j$; $NL_p(n-j)$ is the number of living leaves above the metamer of physiological age $p$, chronological age $j$, at growth cycle $n$; $l_p(n-j)$ is the corresponding metamer length, which is determined by primary growth and related to $q_p^i(n-j)$.

In the mixed approach, the biomass allocated to cambial growth of a metamer of physiological age $p$, chronological age $i$ at growth cycle $n$ is given by
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Illustration of wood distribution into metamers with two modes: left one corresponding to mode D₁ as shown in Eq. 8 and right one corresponding to mode D₂ as shown in Eq. 9.</figcaption></figure>

$$q_p^{\text{ring}}(i,n) = \left( \frac{1-\lambda}{D_1(n)} + \frac{\lambda \cdot NL_p(n-i)}{D_2(n)} \right) \cdot l_p(n-i) \cdot Q^{\text{ring}}(n) \quad (10)$$

where $\lambda$ is a coefficient between 0 and 1.

Fig. 2 illustrates this modelling approach under the assumption that only one metamer is generated at each growth cycle with a constant length equal to one. According to Eq. 10, if $\lambda$ takes null value, only mode $D_1$ takes effect. Total biomass for cambial growth is distributed equally among metamers. If $\lambda$ is equal to one, only mode $D_2$ takes effect. With this last mode, older metamers get more biomass.

## 2.2 Calculation of Stem Mechanical Stability

Due to growth phenomena, modelling of tree biomechanical responses to external loads is a complex problem. Stem and branches are inhomogeneous elements composed of central piths surrounded by stacks of rings that were added at each time step. Moreover, mechanical stresses in tree axes can originate from external, e.g. wind or gravity, and internal, i.e. maturation stresses, actions. In addition, the variable wood density and mechanical properties are factors that enhance the complexity to model tree biomechanical responses to external loads. In our study, we aim at assessing the influence of growth model parameters on wood production and stem mechanical stability by introducing a simple mechanical criterion in a functional-structural model. For our primary work, we calculated the top stem deflection at the

end of the tree growth simulation, considering that the total self weight was applied in one stage only on a slightly leaning stem, i.e. buckling. Such a mechanical criterion, which does not account for progressive growth and gravitropism phenomena due to the formation of reaction wood, is commonly used to assess tree stability (see the use of Euler's buckling criteria by Spatz and Bruechert (2000) for instance), and it was considered as sufficient for the theoretical purpose of this paper.

To simplify the calculations, all the loads are supposed to be integrated and applied at the stem top. Furthermore tree stems were considered slender enough to allow using the beam theory. We used the equilibrium equations adapted to conic beams and proposed by de Reffye (1979) in order to calculate the deflection, i.e. displacement angle (w), of the stem tip (Eqs. 11–19). Fig. 3 illustrates the corresponding variable definitions.

When the stem leans with an angle ($\epsilon$) with regard to the vertical, the equilibrium of moments of the force at position $s$ from the stem base is given by:

$$\begin{align}
E_y I_s \frac{d\theta}{ds} &= F \cos \epsilon \cdot (e - y) + F \sin \epsilon \cdot (L - x) \tag{11} \\
dx &= \cos \theta \cdot ds \\
dy &= \sin \theta \cdot ds
\end{align}$$

where $I_s$ is the secondary moment of the circular cross section at position $s$ from the stem bottom,

$I_s = \pi r_s^4 / 4$, $r_s = r_1(1 - s/H)$, describing the property of the shape and being used to predict the resistance to bending and deflection; $r_1$ is the radius at the bottom of the cone, and $r_s$ is the radius at position $s$; $H$ is the length of the cone, which equals to $L(r_1/(r_1-r_2))$; $L$ is the length of the stem, which is determined by $q^i$; $r_2$ is the radius at the top of the stem; $E_y$ is structural Young's modulus, a measure of the stiffness of an elastic material, less Young's modulus, less stiffness; $\theta$ is the displacement angle with regard to the stem without load; $F$ is the force applied on the stem top; $\epsilon$ is an initial displacement angle at the stem bottom in the vertical plane; $x$ is the projection distance to the plane of the stem without load; $y$ is the deflection distance which is vertical to the stem without load; $e$ is the final deflection distance which is vertical to the stem without load.

We differentiate Eq. 11 with respect to $s$,
$$E_s I_s \frac{d^2\theta}{ds^2} = F \cos\varepsilon \cdot (-\sin\theta) + F \sin\varepsilon \cdot (-\cos\theta) \quad (12)$$

We multiply Eq. 12 by $d\theta/ds$ and integrate it with $\theta$ that varies from 0 to $\omega$. The result is obtained as Eq. 13.

$$\frac{1}{2} E_s I_s \left(\frac{d\theta}{ds}\right)^2 = F \cos\varepsilon(\cos\theta - \cos\omega) - F \sin\varepsilon(\sin\theta - \sin\omega) \quad (13)$$

Suppose

$$K_s = \sqrt{2 \frac{F}{E_y I_s}} = \sqrt{2 \frac{F}{E_y \pi \frac{r_1^4 \left(1 - \frac{s}{H}\right)^4}{4}}} = \sqrt{\frac{2F}{E_y I_1} \frac{1}{\left(1 - \frac{s}{H}\right)^2}} = K_0 \frac{1}{\left(1 - \frac{s}{H}\right)^2} \quad (14)$$

where $I_1$ is the second moment of the circular cross section at the stem base $\pi r_1^4/4$.

Substituting Eq. 14 to Eq. 13, we get Eq. 15,

$$K_0 \frac{ds}{\left(1 - \frac{s}{H}\right)^2} = \frac{d\theta}{\sqrt{\cos(\theta + \varepsilon) - \cos(\omega + \varepsilon)}} \quad (15)$$

The final shape of the stem is given by integrating Eq. 15

$$\int \frac{1}{\left(1 - \frac{s}{H}\right)^2} ds = L = \int_0^\omega \frac{d\theta}{K_0 \sqrt{\cos(\theta + \varepsilon) - \cos(\omega + \varepsilon)}} \quad (16)$$

When $\theta$ reaches the value $w$, the sum of the $ds$ equals to the length of the beam $L$.

Integrating Eq. 16, the final displacement angle of the stem $w$ is as Eqs. 17–19.

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Illustration of variable definitions in Eqs. 11–19. Stem is assumed to be conic as a first approximation. The part of cone with solid line is the tree stem.</figcaption></figure>

$$w = \frac{\sin(\varepsilon) \cdot (1 - \cos(g))}{\cos(\varepsilon) \cdot \cos(g)} \quad (17)$$

with

$$g = \sqrt{\frac{2F}{E_y \cdot I_1}} \cdot H \cdot \sqrt{|\cos(\varepsilon)|} \quad (18)$$

and

$$F = 9.8 \cdot \left( \pi \cdot d \cdot L \cdot \frac{r_1^2 + r_2(r_1 + r_2)}{3} + \sum_{t=1}^{n} Q(t) \right) \quad (19)$$

where $I_1$, $r_1$, $r_2$ are determined by the summation of $q^i$ and $q_{\text{ring}}$ shown in Eq. 2 and Eq. 10; $d$ is the density of the stem; $w$ is the final rotation angle at the top of the deflected stem when a small initial leaning angle $\varepsilon$ is applied to the tree.

## 2.3 Optimization Formula and Algorithm

The first part of our work is to maximize wood production (trunk weight). The formula of this maximization problem $J$ is given by Eq. 20. The optimization variable is the sink strength for cambial growth ($S_{\text{ring}}$).

$$\begin{aligned} J &= Q_{\text{trunk}} = Q_{\text{pith}} + Q_r \\ &= \sum_{k=0}^{n} N_1^i(k) \cdot q_1^i(k) + \sum_{k=0}^{n} \sum_{i=1}^{k} q_1^{\text{ring}}(i, k) \end{aligned} \quad (20)$$
where $Q_{\text{trunk}}$ represents the trunk weight, $Q_{\text{pith}}$ represents the total pith weight of the trunk and $Q_r$ represents the total weight of the rings surrounding the pith on the trunk; the physiological age of the metamers on the trunk is 1.

Substituting Eq. 2, Eqs. 5–7, and Eq. 10 in Eq. 20, the formula of the maximization problem $J$ is rewritten as Eq. 21.

$$J = \sum_{k=0}^{n} N_{1}^{i}(k) P_{1}^{i} \frac{Q(k)}{D(k)} + \sum_{k=0}^{n} \left( S_{\text{ring}} \sum_{p=1}^{P_{m}} N_{p}^{a}(k-1) \frac{Q(k)}{D(k)} \right) \sum_{i=1}^{k} \left[ \left( \frac{1-\lambda}{D_{1}(k)} + \frac{\lambda \cdot NL_{1}(k-i)}{D_{2}(k)} \right) I_{1}(k-i) \right] \quad (21)$$

Note that $S_{\text{ring}}$ is not only a multiplier of the second item of Eq. 21 explicitly, but also implicitly in the equation for plant total demand $D$ as shown in Eq. 7. As plant biomass increment $Q$ is a function of $D$, $S_{\text{ring}}$ appears also implicitly in the equation of $Q$. Hence, it is difficult to get the differentiation information of the objective function of the maximization problem as shown in Eq. 21.

Because heuristic optimization algorithms do not depend on initial values of optimized variables and do not require the differentiation information of objective functions, we turned to heuristic optimization algorithms. Our preference led us to consider population-based methods, which usually allow obtaining quickly acceptable solutions at each iteration. Compared to other population-based heuristic optimization algorithms, PSO has high convergence rate for a wide range of optimization problems (Kennedy and Eberhart 2001). Inspired by these features, we use the PSO method in order to find the optimal solution of $S_{\text{ring}}$ to get the maximal trunk weight. Among the variations of PSO proposed by researchers to improve the performance of the initial PSO (Kennedy and Eberhart 1995), we chose passive congregation PSO which is proposed by He et al. (2004), thanks to its simplicity, generalization of parameters and high accuracy.

The population of PSO is called *Swarm*, and each individual in the population is called *Particle*. There are two properties to describe particles: position and velocity. The position of each particle represents a possible solution of the optimization problem, and the velocity is used to change the position. The velocity calculation and the position modification are as Eq. 22. The last item of the equation for calculating particle velocities is the distinction from the initial PSO (Kennedy and Eberhart 1995), which is used to avoid trapping in local extremum.

$$\begin{gathered} v_{ij}^{k+1} = \sigma \cdot v_{ij}^{k} + c_1 \text{rand}_1 (PI_{ij} - xp_{ij}^{k}) \\ + c_2 \text{rand}_2 (PG_j - xp_{ij}^{k}) + c_3 \text{rand}_3 (PR_j - xp_{ij}^{k}) \\ xp_{ij}^{k+1} = xp_{ij}^{k} + v_{ij}^{k+1} \end{gathered} \quad (22)$$

$i = 1, 2, \dots, \text{Num}$

$j = 1, 2, \dots, Nd$

where *Num* is the number of particles in the population; *Nd* is the number of variables of the problem; $v_{ij}^k$ is the velocity of the $j$th variable of the $i$th particle at iteration $k$; $PI_{ij}$ is the best position of variable $j$ recorded by the $i$th particle during the previous iterations; $PG_j$ is the global best position of variable $j$ in the swarm; $PR_j$ is the best position of variable $j$ recorded by a particle which is randomly selected in the swarm; $xp_{ij}^k$ is the current position of variable $j$ of particle $i$ at the $k$th iteration; $\sigma$ is inertia weight that varies linearly from 0.9 to 0.7; $c_1$, $c_2$ and $c_3$ are acceleration coefficients, where the values are 0.5, 0.5 and 0.6, respectively; $rand_1$, $rand_2$ and $rand_3$ are uniformly distributed random values between 0 and 1.

## 3 Results

The GreenLab parameters used in this study (see Table 1) have been chosen based on the previous published research works carried out on trees (Guo et al. 2006, Letort et al. 2008a). The constitutive material of a tree stem is supposed to be isotropic. It was observed on several tree species that the structural Young's modulus is only slightly different as well along the stem (Niklas 1997ab) as within the stand (Brüchert et al. 2000, Alméras et al. 2004). Therefore the structural Young's modulus is set to be constant in our study. The effect of the sink strength for
<table><caption>**Table 1.** Parameter values of the GreenLab model for tree.</caption><thead><tr><th>Parameter</th><th>Definition</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td><i>E</i></td><td>Environmental factor</td><td>0.1</td><td>g/cm<sup>2</sup></td></tr><tr><td><i>Sp</i></td><td>Total ground projection area available of the crown for plant</td><td>1E+6 (isolated)<br>1E+5 (high density)</td><td>cm<sup>2</sup></td></tr><tr><td><i>K</i></td><td>Light extinction coefficient</td><td>1</td><td>--</td></tr><tr><td><i>slw</i></td><td>Specific leaf weight</td><td>0.04</td><td>g/cm<sup>2</sup></td></tr><tr><td><i>ta</i></td><td>Leaf functioning duration</td><td>1</td><td>growth<br>cycle</td></tr><tr><td><i>P</i><sub>m</sub></td><td>Maximum physiological age</td><td>5</td><td>--</td></tr><tr><td><i>P</i><sub>p</sub><sup>a</sup></td><td>Sink strength of leaf with physiological age (PA) <i>p</i></td><td>1, 0.8, 0.6, 0.4, 0.2<br>from PA 1 to PA <i>P</i><sub>m</sub></td><td>--</td></tr><tr><td><i>P</i><sub>p</sub><sup>i</sup></td><td>Sink strength of pith with physiological age (PA) <i>p</i></td><td>0.4, 0.32, 0.2, 0.1, 0.08<br>from PA 1 to PA <i>P</i><sub>m</sub></td><td>--</td></tr><tr><td><i>E</i><sub>y</sub></td><td>Structural Young's modulus</td><td>10</td><td>GPa</td></tr></tbody></table>

“--” means no unit is associated to the parameter.

cambial growth ($S_{ring}$) and of the coefficient ($\lambda$) on stem wood production and tree stability has been investigated. $S_{ring}$ and $\lambda$ have distinct effects: $S_{ring}$ determines the total biomass allocated to wood formation at the whole-plant level while $\lambda$ drives the way this total amount of biomass is partitioned in the tree architecture. The effects of these two parameters on tree growth being distinct, they were considered independently in the optimization process.

## 3.1 Influence of the Sink Strength for Cambial Growth ($S_{ring}$) on Wood Production

Simulated trunk weight and tree height with regard to $S_{ring}$ are shown in Fig. 4(a). Tree height decreases monotonously with increasing in $S_{ring}$. This result is straightforward if we consider the GreenLab equations describing biomass allocation in the tree. Tree height is indeed deduced

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** (a) Variations of trunk weight and tree height with regard to $S_{ring}$ (y-axis is in logarithm scale.). (b), (c) and (d) 3D image of tree: in (b) $S_{ring}=0$, trunk weight is 2.03 kg, height is 19.79 m; in (c) $S_{ring}=0.1$, trunk weight is 187 kg, height is 13.54 m; in (d) $S_{ring}=0.5$, trunk weight is 0.11 kg, height is 0.39 m.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Ratio of biomass for cambial growth and for piths $Q_{wood}$ to tree biomass $Q$ with respect to $S_{ring}$.</figcaption></figure>

from the summation of trunk metamer lengths. Metamer length depends on pith weight and is therefore a function of the amount of biomass allocated to primary growth. Due to mechanisms of sink competition, this quantity is reduced when the demand of the biomass compartment for cambial growth increases with increasing in $S_{ring}$. According to Eq. 2 and Eqs. 5–7, increasing in $S_{ring}$ leads to increasing plant demand $D$, which leads to the decrease of the amount of biomass allocated to leaves and piths. Hence increase in $S_{ring}$ results in lower pith weights which lead to shorter metamers.

Trunk weight is determined by summing the pith weights of all its constitutive metamers and the total weight of rings surrounding these piths. Increasing $S_{ring}$ has two contradictory effects: a negative effect on pith weights and a positive effect on ring weights. With increasing in $S_{ring}$, a larger amount of the tree total biomass is allocated to cambial growth at the expense of primary growth and in particular of leaf growth. Hence leaf area index (LAI) decreases, which reduces in turn the plant photosynthetic potential. This means that an optimal $S_{ring}$ value can be found that maximizes the trunk weight (Fig. 4(a)). Under the condition that the other parameters of the GreenLab model are fixed, the PSO algorithm with 20 particles in

the swarm provides after 50 iterations the optimal value of $S_{ring}$, which is 0.1, that maximizes the trunk weight. Fig. 4(b–d) represent the 3D shapes of three simulated trees with different values of the sink strength for cambial growth.

Fig. 5 shows that when sink strength for cambial growth takes a value exceeding 50, 95% of the biomass produced by the tree is allocated to cambial growth and piths. As a result, not enough biomass remains for leaves, thus inducing a slow-down of tree growth. This result coincides with the previous results of Mäkelä (1986). Therefore, far from generating higher yield, increasing tree cambial sink strength ($S_{ring}$) has a negative effect on the tree growth rate. It leads to shorter tree height, branches and sizes of leaves which can reduce the tree performances in situations of high competition with neighbour trees, due to less light amount intercepted by the leaves.

Fig. 6 enables to compare relative growth rates of wood (pith and rings) to relative tree growth rate. Relative growth rate is defined as the ratio $(V_n - V_{n-1}) / V_{n-1}$, $V_n$ being the biomass allocated to wood or the tree biomass production at growth cycle $n$. As shown in Eqs. 4–7, the biomass production of plant $Q$ and the total biomass for cambial growth are both related to the product between numbers of organs and the ratio of bio-
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Comparison of relative growth rates between wood growth (pith and rings) and tree growth. (a) $S_{\text{ring}}$ is set to the optimal value 0.1 (b) $S_{\text{ring}}=1$.</figcaption></figure>

mass production $Q$ to total demand of plant $D$. In addition, the sink strength for cambial growth $S_{\text{ring}}$ is set to be constant. Therefore, the relative growth rates of wood and of photosynthetic potential are similar. At the first growth cycle, as plant biomass production and biomass allocated to wood increase from 0, their growth rates are both 100%. After that, the growth rates decrease due to the competition of cambial growth against

the growth of leaves for photosynthesis. For the particular tree with the maximal trunk weight and with the optimal balance between cambial growth and photosynthetic potential, the relative growth rates for wood growth and tree growth are zero at the end of plant growth, as shown in Fig. 6(a). However, for larger $S_{\text{ring}}$, relative growth rates tend to be negative as shown in Fig. 6(b), and it implies that trees are under suppression.
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** (a) Variation of the trunk weight with respect to λ. (b) and (c) 3D image of tree with the optimal value of sink strength for cambial growth: S<sub>ring</sub>=0.1: in (b) λ=0, trunk weight is 1.63 kg, tree height is 13.54 m; in (c) λ=1, trunk weight is 187 kg, tree height is 13.54 m.</figcaption></figure>

## 3.2 Effect of λ on Wood Production

The partition coefficient λ in Eq. 10 determines the distribution of total biomass assigned to cambial growth to each metamer. It does not change the tree height which is controlled by primary growth. Eq. 10 implies that λ should be equal to 1 in order to allocate the maximal proportion of wood to the stem. In that case, the proportion of biomass allocated to the stem is higher than that allocated to branches. Fig. 7 shows the variation of the trunk weight as λ increases through simulation and gives two examples to illustrate the impact of λ, all the other parameters being fixed and S<sub>ring</sub> being set to be the optimal value. As the simulation results showed that the trunk weight increases monotonously with increasing in λ, the use of the optimization algorithm to optimize trunk weight with respect to λ is not required.

## 3.3 Stem Buckling Criteria

As described in section 2, the mechanical criterion of stem stability is represented by the final rotation angle at the stem tip. It depends on

tree weight, tree height and diameters, and stem conicity. The variation of this final rotation angle versus stem diameter at breast height (DBH) with increasing in λ is shown in Fig. 8(a), where S<sub>ring</sub> is set to be the optimal value 0.1. DBH decreases with decreasing in λ, due to the dominant effect of mode D<sub>1</sub>. Meanwhile, the stem tip tends to bend more.

Fig. 8(b) shows the variations of stem diameter versus the final rotation angle at the stem tip with increasing in S<sub>ring</sub>, where λ is set to be 1. As shown in Fig. 4(a), when S<sub>ring</sub> is larger than 0.3, tree height is less than 1.30 m. In this particular case, the diameter at the stem tip is considered instead of DBH. Contrary to λ, S<sub>ring</sub> influences the photosynthesis process and tree height. Consequently the load (self-weight) applied to the stem as well as the tree height vary with increasing in S<sub>ring</sub>. Nevertheless, Fig. 8(a) and Fig. 8(b) reveal that the stem DBH and the final rotation angle at the stem tip always have opposite slopes, i.e. increasing bending angle always coincides with decreasing DBH and reciprocally. Two examples are given to illustrate the impact of the coefficient λ on stem stability as shown in Fig. 9.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Relationship between stem diameter and bending angle at the stem tip with respect to (a) lambda (λ) (b) sink strength for cambial growth (S<sub>ring</sub>). (y-axis is in logarithm scale.)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** 3D image of bending stems with respect to λ, where S<sub>ring</sub>=0.1, tree weight is 1120 kg. (a) λ=0, trunk weight is 1.63 kg, bending angle at the stem tip is 176.2°; (b) λ=1, trunk weight is 187 kg, bending angle at the stem tip is 0°.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** (a) Difference between diameters at the bottom and at the top of the stem with respect to $\lambda$. (b) Simulation results of bending angle at the stem tip and stem diameter at breast height for old trees growing in high stand density (y-axis is in logarithm scale). The maximal total ground projection area available for an old tree in high density is 10 m² (approximately 1000 plants/ha).</figcaption></figure>

## 3.4 Impact of Stand Density on Tree Stability

Fig. 10(a) shows differences in old tree diameters considering different densities. In the GreenLab model, density is related to total ground projection area available of the crown for plant *Sp*. From Fig. 10(a), we see that with respect to $\lambda$, the dif-

ferences between diameters at the stem base and at the stem tip for an old tree in high stand density (approximately 1000 plants/ha in Fig. 11) are all less than the ones for an old isolated tree. It implies that as stand density increases, the stem profile tends to be more cylindrical. Furthermore, for any density value, DBH monotonously increases with respect to $\lambda$, while the rotation angle at the
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Comparison of 3D images of an isolated tree on the left and a tree growing in high density (approximately 1 000 plants/ha) on the right, where the other parameter values are the same, especially S<sub>ring</sub> is 0.1 and λ is 1. The heights of an isolated tree and a tree growing in high density are 13.54 m and 5.60 m, respectively; DBHs are 16.18 cm and 6.56 cm, respectively.</figcaption></figure>

stem tip decreases, as shown in Fig. 8(a) and Fig. 10(b). This result is in agreement with the beam theory in mechanics where the beam stiffness is proportional to beam diameter to the fourth power. Two 3D images of an old isolated tree and an old tree growing in high density are shown in Fig. 11. These two trees share the same parameter values except one: the total ground projection area available of the crown for plant which is related to stand density.

# 4 Discussion

In this paper, a population-based optimization algorithm has been applied to the functional-structural plant growth model GreenLab, which can be written under the mathematical form of a discrete dynamic system. The aim was to optimize a particular problem in the domain of forestry, namely maximization of wood production and

tree stability. Two parameters of GreenLab were considered: S<sub>ring</sub> and λ that drive respectively biomass allocation and partitioning for cambial growth. We first discuss the influence of these parameters on the variables of interest in this work, i.e. wood production, tree stability.

Contrary to the common idea that increasing sink strength for cambial growth would increase the final wood production, we showed in this paper that there is an optimal value of the sink strength. This sink strength, S<sub>ring</sub>, determines the quantity of total biomass allocated to the rings of trunk and of branches surrounding the piths. It controls the global allocation of the tree biomass into cambial growth. The level of competition of cambial growth against the other organs is reflected by the ratio of biomass production to biomass demand as shown in Eq. 2 and Eqs. 4–7 (see also Mathieu et al. 2009). In terms of source-sink dynamics, an increase in S<sub>ring</sub> has two contradictory effects. On one hand, it induces a slowdown of the plant photosynthetic potential.
It results in less biomass allocated to leaves, and consequently in a decrease in leaf area index. Meanwhile, it reduces biomass allocated to piths (pith weight). On the other hand, it increases the total weight of rings surrounding piths. Therefore, to maximize wood production (pith and rings), a balance had to be found between cambial growth and primary growth.

In the literature, for a branching point on the trunk, one part of the assimilates produced by a leaf goes downward towards the root system and the other part goes upward to the crown. However, the Pressler law only reproduces the downward propagation of assimilates in the plant hydraulic system. In our model, we introduced an original way of modelling more general cases: the coefficient $\lambda$ drives the proportion of upward to downward propagation. The value of $\lambda$ can be estimated from real data for a given tree species, by fitting the stem profile (Guo et al. 2006, Letort et al. 2008a). According to our simulation results, in order to maximize the trunk weight, the optimal value of $\lambda$ is 1, which leads to minimizing the load of branches. Hence, this preferential allocation of biomass to the trunk penalizes the cambial growth of branches and thus reduces their hydraulic conductivity. The case where $\lambda$ is equal to 0 is favourable to water conduction in branches, whereas it has negative effect on stem stability; and vice versa. Thus, the value of this parameter $\lambda$ has also indirect effect on water conductivity, which should be considered as additional constraints when applying the methods presented in this paper to real case studies. This can be done through multi-objective optimization, i.e. maximization of stem stability and water conductivity of branches simultaneously with respect to $\lambda$. Considering the contradictory effects of $\lambda$ on these two objectives can be an interesting extension of this work.

The analysis of the simulated results reveals emergent properties of GreenLab that are in accordance with the literature. No matter considering variations of $\lambda$ or $S_{ring}$, as DBH increases, the stem bends less and less. The same relationship between DBH and stem stiffness is obtained for trees both in isolated and dense states, as illustrated in Fig. 8(a) and Fig. 10(b), respectively. With the simplified mechanical assumptions, because the bending angle of a stem resulting

from a load is negatively proportional to the second moment of area which is biquadratic of diameter, the stiffness of a stem is mainly influenced by stem dimensions. Therefore, a bigger diameter coincides with a smaller bending angle. This phenomenon is observed for most kinds of trees (Wiemann and Williamson 1989ab, Woodcock and Shier 2003, Ancelin et al. 2004b). The results coincide with the hypothesis in Wiemann and Williamson (1989ab): radial increases are related to mechanical support; stiffness increases as diameter increases. The increase in stiffness conferred with an increase in diameter was observed, no matter whether a tree is isolated or it is in high density. Compared to isolated trees, the stem profile of trees grown in high density tends to be cylinder (Fig. 10(a)). It is in agreement with common observations in forestry (Assmann 1970).

The optimization results that we presented open interesting perspectives for applications in breeding programs. In the present work, it was assumed that the model parameters were genotypic, i.e. they only depend on the plant genes and are independent of environmental influences. A preliminary study to validate this assumption on real tree species is presented in Letort et al. (2008a) where two beech trees (*Fagus sylvatica*) growing in different local environments were fitted with the same set of genotypic parameter values. Recently, efficient methods have been developed in the field of quantitative genetics to identify particular loci on the plant chromosomes that contribute to phenotypic traits (de Vienne 1998). These methods rely on establishing statistical correlations between quantitative traits that can be measured on plants (e.g. yield, duration, height) and the values of a set of particular genes (markers). Several authors have emphasized the potential benefits of coupling these methods to models of plant growth in order to unravel the complex genotype × environment interactions and thus to improve genetic selection (Hammer et al. 2002, Tardieu 2003, Yin et al. 2004, Hammer et al. 2006). Models such as GreenLab can play the role of intermediate link between genotype and phenotype. It requires identifying the genotype loci that have influence on the parameters of GreenLab instead of the classical breeding traits. A simulation study of these methods applied to
the case of GreenLab can be found in Letort et al. (2008b): a simple genetic model was built to determine the model parameters as functions of the genotype of the plant considered. In this context, defining a set of optimized parameter values under given objectives is equivalent to defining ideotypes. Ideotype is the set of desirable traits that a plant should present to enhance yield or any other objective trait under specified conditions (Donald 1968). To benefit from the recent advances in plant growth modelling, this set of traits should include the values of model parameters (Dingkuhn et al. 2007, Letort et al. 2008c). It implies that modellers should be able to provide a set of optimal parameter values in response to specified objectives, which is precisely what was done in this paper. This is of crucial importance in the domain of plant breeding as it provides new criteria to improve genetic selection. In this theoretical study, we present a generic method without considering any particular species. To be as general as possible, it was chosen to set a wide range of variations for each parameter value. Therefore, depending on the species considered in application, the optimum found may not be biologically relevant. Even though the perfect ideotype may not be possible to reach, it can be considered as ideal values that could guide the breeding efforts.

This paper aimed at exploring new applications of functional-structural growth models, in particular in the framework of optimizing wood production in conjunction with tree stability. This last criterion was considered as an indicator of wood quality and is known to be affected by tree biomechanical responses, e.g. formation of reaction wood involving tree gravitropic movements. It is also well known that besides internal factors (internal biomechanical stress), the shape of trees is influenced by environmental factors (wind, storm, etc). James et al. (2006) stated that wind was the main external force influencing stem shape. In our theoretical work, we introduced a simplified biomechanical model where neither the impact of wind nor of maturation stress was considered. Furthermore, we hypothesized that i) the number of metamers appearing at each growth cycle is fixed, ii) the tree is only subjected to its self-loading, and iii) the stem has a smooth conic profile. Despite these hypotheses, the numerical

approach showed promising perspectives that can be used further in breeding programs. This method is generic and could be applied to any tree species. To conclude, comparing qualitatively the knowledge in the literature, studying the ideal case of a tree with simplified interactions between tree growth and biomechanics provided reasonable and interesting results that could be extended to more complex cases, including for instance environmental factors and retroactive influences of biomass production on topology in further studies.

## Acknowledgement

We are grateful to two anonymous reviewers for their suggestions. We thank the French Embassy in China, the National Natural Science Foundation of China (NSFC, # 60723005, 60703043), and Chinese 863 program (2008AA10Z218) for their support. AMAP (Botany and Computational Plant Architecture) is a joint research unit which associates CIRAD (UMR51), CNRS (UMR5120), INRA (UMR931), IRD (2M123), and Montpellier 2 University (UM27); http://amap.cirad.fr/.

## References

* Allen, M., Prusinkiewicz, P. & DeJong, T. 2005. Using L-systems for modeling source-sink interactions, architecture and physiology of growing trees, the L-peach model. New Phytologist 166: 869–880.

* Alméras, T., Costes, E. & Salles, J.C. 2004. Identification of biomechanical factors involved in stem shape variability between Apricot Tree Varieties. Annals of Botany 93: 455–468.

* Ancelin, P., Fourcaud, T. & Lac, P. 2004a. Modelling the biomechanical behaviour of growing trees at the forest stand scale. Part 1: Development of an incremental Transfer Matrix Method and application to simplified tree structures. Annals of Forest Science 61: 263–275.

* —, Courbaud, B. & Fourcaud, T. 2004b. Development of an individual tree-based mechanical model to predict wind damage within forest stands. Forest Ecology and Management 203: 101–121.
Assmann, E. 1970. The principles of forest yield study. Pergamon Press, Oxford.

Barthélémy, D., Caraglio, Y. & Costes, E. 1997. Architecture, gradients morphogénétiques et âge physiologique chez les végétaux. In: Bouchon, J. (ed.). Modélisation et simulation de l'architecture des végétaux, Science Update INRA. p. 89-136.

— & Caraglio, Y. 2007. Plant architecture : a dynamic, multilevel and comprehensive approach to plant form, structure and ontogeny. Annals of Botany 99: 375-407.

Brüchert, F., Becker, G. & Speck, T. 2000. The mechanics of Norway spruce [Picea abies (L.) Karst]: mechanical properties of standing trees from different thinning regimes. Forest Ecology and Management 135: 45-62.

Carson, E., Feng, D.D., Pons, M.N., Soncini-Sessa, R. & van Straten, G. 2006. Dealing with bio- and ecological complexity: Challenges and opportunities. Annual Reviews in Control 30: 91-101.

Cournède, P.H., Kang, M.Z., Mathieu, A., Barczi, J.F., Yan, H.P., Hu, B.G. & de Reffye, P. 2006. Structural factorization of plants to compute their functional and architectural growth. Simulation 82: 427-438.

—, Mathieu, A., Houllier, F., Barthélémy, D. & 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. Annals of Botany 101: 1207-1219.

Deleuze, C. & Houllier, F. 1997. A transport model for tree ring width. Silva Fennica 31(3): 239-250.

— & Houllier, F. 2002. A flexible radial increment taper equation derived from a process-based carbon partitioning model. Annals of Forest Science 59: 141-154.

Denne, M.P., Atkinson, C.J. & Dodd, R.S. 1994. Quantification of Trends in Wood Production within Trees of Silver Birch (Betula pendula Roth.). Annals of Botany 73: 655-664.

Deussen, O., Hanrahan, P., Lintermann, B., Mech, R., Pharr, M. & Prusinkiewicz, P. 1998. Realistic Modeling and Rendering of Plant Ecosystems. In: 25th annual conference on Computer graphics and interactive techniques. Association for Computing Machinery SIGGRAPH, New York, NY, USA. p. 275-286.

Dingkuhn, M., Luquet, D., Clément-Vidal, A., Tambour, L., Kim, H.K. & Song, Y.H. 2007. Is plant

growth driven by sink regulation? In: Spiertz, J.H.J., Struik, P.C., Van Laar, H.H., (eds.). Scale and complexity in plant systems research: gene-plant-crop relations. p. 157-170.

Donald, C.M. 1968. The breeding of crop ideotypes. Euphytica 17: 385-403.

Dornbusch, T., Wernecke, P. & Diepenbrock, W. 2007. A method to extract morphological traits of plant organs from 3D point clouds as a database for an architectural plant model. Ecological Modelling 200: 119-129.

Forestry commission in England. http://www.forestry.gov.uk/england-woodfuel.

Fourcaud, T. & Lac, P. 2003. Numerical modelling of shape regulation and growth stresses in trees I. An incremental static finite element formulation. Trees-Structure and Function 17: 23-30.

—, Blaise, F., Lac, P., Castera, P. & 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-Structure and Function 17: 31-39.

—, Zhang, X.P., Stokes, A., Lambers, H. & Koerner, C. 2008. Plant growth modelling and applications: The increasing importance of plant architecture in growth models. Annals of Botany 101: 1053-1063.

Fournier, M., Stokes, A., Coutand, C., Fourcaud, T. & Moulia, B. 2006. Tree biomechanics and growth strategies in the context of forest functional ecology. In: Herrel, A., Speck, T. & Rowe, N. (eds.). Ecology and biomechanics - a mechanical approach to the ecology of animals and plants. Boca Raton: CRC Press Taylor & Francis. p. 1-33.

Fox, J.C., Bi, H. & Ades, P.K. 2008. Modelling spatial dependence in an irregular natural forest. Silva Fennica 42(1): 35-48.

Greene, N. 1989. Voxel space automata: modeling with stochastic growth processes in voxel space. In: 16th annual conference on Computer graphics and interactive techniques. Association for Computing Machinery SIGGRAPH, New York, NY, USA. p. 175-184.

Guo, H., Letort, V., Hong, L., Fourcaud, T., Cournède, P.H., Lu, Y. & de Reffye, P. 2006. Adaptation of the GreenLab model for analyzing sink-source relationships in Chinese pine saplings. In: Fourcaud, T. & Zhang, X.P. (eds.). Plant growth modeling, and applications. IEEE Computer Society, Los Alamitos, California. p. 236-243.
Hammer, G., Kropff, M.J., Sinclair, T.R. & Porter, J.R. 2002. Future contributions of crop modelling—from heuristics and supporting decision making to understanding genetic regulation and aiding crop improvement. European Journal of Agronomy 18: 15–31.

—, Cooper, M., Tardieu, F., Welch, S., Walsh, B., Van Eeuwijk, F., Chapman, S. & Podlich, D. 2006. Models for navigating biological complexity in breeding improved crop plants. Trends in Plant Science 11: 587–593.

He, S., Wu, Q.H., Wen, J.Y., Saunders, J.R. & Paton, R.C. 2004. A particle swarm optimizer with passive congregation. BioSystems 78: 135–147.

Houllier, F., Leban, J.M. & Colin, F. 1995. Linking growth modelling to timber quality assessment for Norway spruce. Forest Ecology and Management 74: 91–102.

James, K.R., Haritos, N. & Ades, P.K. 2006. Mechanical stability of trees under dynamic loads. American Journal of Botany 93: 1522–1530.

Jayawickrama, K.J.S. 2001. Breeding radiata pine for wood stiffness: review and analysis. Aust. Forest. 64: 51–56.

Kang, M.Z., Cournède, P.H., de Reffye, P., Auclair, D. & Hu, B.G. 2008. Analytical study of a stochastic plant growth model: Application to the GreenLab model. Mathematics and Computers in Simulation 78: 57–75.

Kennedy, J. & Eberhart, R.C. 1995. Particle swarm optimization. In: Proc. IEEE Conf. on Neural Networks, IV, Piscataway, NJ. p. 1942–1948.

— & Eberhart, R.C. 2001. Swarm intelligence. Morgan Kaufmann Publishers.

Lasserre, J.P., Mason, E.G. & Watt, M.S. 2005. The effects of genotype and spacing on Pinus radiata [D. Don] corewood stiffness in an 11-year old experiment. Forest Ecology and Management 205: 375–383.

Letort, V., Cournède, P.H., Mathieu, A., de Reffye, P. & Constant, T. 2008a. Parametric identification of a functional-structural tree growth model and application to beech trees (Fagus sylvatica). Functional Plant Biology 35: 951–963.

—, Mahe, P., Cournède, P.H., de Reffye, P. & Cournois, B. 2008b. Quantitative Genetics and Functional-Structural Plant Growth Models: Simulation of Quantitative Trait Loci Detection for Model Parameters and Application to Potential Yield Optimization. Annals of Botany 101: 1243–1254.

—, Mahe, P., Cournède, P.H., de Reffye, P. & Cournois, B. 2008c. Optimizing plant growth model parameters for genetic selection based on QTL mapping. In: Fourcaud, T., Zhang, X.P. (eds.). Plant growth Modeling, and Applications, IEEE Computer Society (Los Alamitos, California). p. 16–21.

Mäkelä, A. 1986. Implications of the pipe model theory on dry matter partitioning and height growth in trees. Journal of Theoretical Biology 123: 103–120.

—, Landsberg, J., Ek, A.R., Burk, T.E., Ter-Mikaelian, M., Agren, G.I., Oliver, C.D. & Puttonen, P. 2000. Process-based models for forest ecosystem management: current state of the art and challenges for practical implementation. Tree Physiology 20: 289–298.

Mathieu, A. 2006. Essai sur la modélisation des interactions entre la croissance et le développement d'une plante – Cas du modèle GreenLab. PhD Thesis. Ecole Centrale de Paris.

—, Cournède, P.H., Barthélémy, D. & de Reffye, P. 2008. Rhythms and Alternating Patterns in Plants as Emergent Properties of a Model of Interaction between Development and Functioning. Annals of Botany 101: 1233–1242.

—, Cournède, P.H., Letort, V., Barthélémy, D. & de Reffye, P. 2009. A dynamic model of plant growth with interactions between development and functional mechanisms to study plant structural plasticity related to trophic competition. Annals of Botany 103: 1173–1186.

Moulia, B., Coutand, C. & Lenne, C. 2006. Posture control and skeletal mechanical acclimation in terrestrial plants: implications for mechanical modeling of plant architecture. American Journal of Botany 93: 1477–1489.

Nepveu, G. 1993. The possible status of wood quality in oak breeding programs (Quercus Petreae Liebl and Quercus robur L). Annals of Forest Science 50 Suppl 1: 338–394.

Niklas, K.J. 1997a. Mechanical properties of Black locust (Robinia pseudoacacia L.) wood. Size- and age-dependent variations in sap- and heartwood. Annals of Botany 79: 265–272.

— 1997b. Size- and age-dependent variation in the properties of sap- and heartwood in Black locust (Robinia pseudoacacia L.). Annals of Botany 79: 473–478.

Peng, C. 2000. Growth and yield models for uneven-
aged stands: past, present and future. Forest Ecol-
ogy and Management 132: 259-279.

Perttunen, J., Sievänen, R. & Nikinmaa, E. 1998. LIGNUM: a model combining the structure and the functioning of trees. Ecological Modelling 108: 189-198.

Prusinkiewicz, P. 2004. Modeling plant growth and development. Current Opinion in Plant Biology 7: 79-83.

Repola, J. 2008. Biomass equations for birch in Finland. Silva Fennica 42(4): 605-624.

de Reffye, P. 1979. Modélisation de l'architecture des arbres par des processus stochastiques. Simulation spatiale des modeles tropicaux sous l'effet de la pesanteur. Application au coffea robusta. Thèse à l'université de Paris sud Centre d'Orsay.

—, Fourcaud, T., Blaise, F., Barthélémy, D. & Houllier, F. 1997. A Functional Model of Tree Growth and Tree Architecture. Silva Fennica 31(3): 297-311.

— & Hu, B.G. 2003. Relevant qualitative and quantitative choices for building an efficient dynamic plant growth model: GreenLab Case. In: Hu, B.G. & Jaeger, M. (eds.). Plant growth modeling and applications proceedings - PMA03. Tsinghua University Press and Springer, Beijing, China. p. 87-107.

—, Heuvelink, E., Barthélémy, D. & Cournède, P.H. 2008. Plant growth models. In: Jøgensen, S.E., Fath, B.D. (eds.). Ecological models. Vol. [4] of Encyclopedia of Ecology, 5 vols. p. 2824-2837.

Russia to go ahead with export tax increase on wood:
Fradkov. Available at: http://www.eubusiness.com/
news.live/1180468812.5.

Shinozaki, K., Yoda, K., Hozumi, K. & Kira, T. 1964. A quantitative analysis of plant form – the pipe model theory. I Basic analysis. Japanese Journal of Ecology 14: 97-105.

Sievänen, R., Nikinmaa, E., Nygren, P., Ozier-Lafontaine, H., Perttunen, J. & Hakula, H. 2000. Components of functional-structural tree models. Annals of Forest Science 57: 399-412.

Spatz, H.C. & Bruechert, F. 2000. Basic biomechanics of self-supporting plants: wind loads and gravitational loads on a Norway spruce tree. Forest Ecology and Management 135: 33-44.

Tardieu, F. 2003. Virtual plants: modelling as a tool for the genomics of tolerance to water deficit. Trends in Plant Science 8: 9-14.

de Vienne, D. 1998. Les marqueurs moléculaires en génétique et biotechnologies végétales. Paris: INRA Editions.

Wang, X.P., Divos, F., Pilon, C., Brashaw, B.K., Ross, R.J. & Pellerin, R.F. 2004. Assessment of Decay in Standing Timber Using Stress Wave Timing Nondestructive Evaluation Tools: A Guide for Use and Interpretation. [Online report]. Available at: http://www.fpl.fs.fed.us/documents/fplgtr/fpl_gtr147.pdf.

Wiemann, M.C. & Williamson, G.B. 1989a. Radial gradients in the specific gravity of wood in some tropical and temperate trees. Forest Science 35: 197-210.

— & Williamson, G.B. 1989b. Wood specific gravity gradients in tropical dry and montane rain forest trees. American Journal of Botany 76: 924-928.

Woodcock, D.W. & Shier, A.D. 2003. Does canopy position affect wood specific gravity in temperate forest trees? Annals of Botany 91: 529-537.

Yan, H.P., Kang, M.Z., de Reffye, P. & Dingkuhn, M. 2004. A Dynamic, Architectural Plant Model Simulating Resource-dependent Growth. Annals of Botany 93: 591-602.

Yin, X., Struik, P. & Kropff, M. 2004. Role of crop physiology in predicting gene-to-phenotype relationships. Trends in Plant Science 9: 426-432.

Zhang, G., Jiang, H., Niu, G.D., Liu, X.W. & Peng, S.L. 2006. Simulating the dynamics of carbon and nitrogen in litter-removed pine forest. Ecological Modelling 195: 363-376.

Zhao, X., de Reffye, P., Barthelemy, D. & Hu, B.G. 2003. Interactive Simulation of Plant Architecture Based on a Dual-Scale Automaton Model. In: Hu, B.G., Jaeger, M. (eds.). Plant Growth Modeling and Applications Proceedings - PMA03. Tsinghua University Press and Springer, Beijing, China. p. 144-153.

Total of 70 references
**Appendix. Symbol list.**

<table><tr><td>Q(n)</td><td>biomass increment of an individual plant at growth cycle n</td></tr><tr><td>Q<sub>seed</sub></td><td>seed biomass</td></tr><tr><td>D(n)</td><td>total demand of plant at growth cycle n</td></tr><tr><td>E(n)</td><td>environmental factor at growth cycle n</td></tr><tr><td>Sp</td><td>total ground projection area available of the crown for plant</td></tr><tr><td>k</td><td>light extinction coefficient</td></tr><tr><td>S(n)</td><td>total green leaf surface area at growth cycle n</td></tr><tr><td>q<sup>o</sup><sub>p</sub>(n)</td><td>biomass of organ o of physiological age p, when the plant age is n</td></tr><tr><td>P<sup>o</sup><sub>p</sub></td><td>sink strength of organ o of physiological age p</td></tr><tr><td>QB(n)</td><td>total green leaf weight at growth cycle n</td></tr><tr><td>N<sup>o</sup><sub>p</sub>(n)</td><td>number of organ o of physiological age p, initiated at growth cycle n</td></tr><tr><td>o</td><td>symbol representing organ types, e.g. a represents a leaf and i represents a pith</td></tr><tr><td>P<sub>m</sub></td><td>maximum physiological age</td></tr><tr><td>ta</td><td>leaf functioning duration</td></tr><tr><td>slw</td><td>specific leaf weight</td></tr><tr><td>Q<sub>ring</sub>(n)</td><td>total biomass allocated to cambial growth at the whole-plant level at growth cycle n</td></tr><tr><td>D<sub>ring</sub>(n)</td><td>total demand for cambial growth at growth cycle n, at the whole-plant level</td></tr><tr><td>D<sub>1</sub>, D<sub>2</sub></td><td>demand of cambial growth for biomass allocation to metamer</td></tr><tr><td>S<sub>ring</sub></td><td>sink strength for cambial growth</td></tr><tr><td>l<sub>p</sub>(n)</td><td>length of the metamer of physiological age p, initiated at growth cycle n</td></tr><tr><td>q<sup>ring</sup><sub>p</sub>(j,n)</td><td>biomass allocated to cambial growth of the metamer of chronological age j and physiological age p, when the plant age is n</td></tr><tr><td>NL<sub>p</sub>(n)</td><td>number of living leaves above the metamer of physiological age p, initiated at growth cycle n</td></tr><tr><td>λ</td><td>coefficient determining the way the biomass assigned to cambial growth is allocated to each metamer</td></tr><tr><td>E<sub>y</sub></td><td>structural Young's modulus</td></tr><tr><td>I<sub>s</sub></td><td>secondary moment of the circular cross section at position s from the stem bottom</td></tr><tr><td>I<sub>1</sub></td><td>secondary moment of the circular cross section at the stem base</td></tr><tr><td>r<sub>s</sub></td><td>radius at position s from the stem bottom</td></tr><tr><td>r<sub>1</sub></td><td>radius at the stem base</td></tr><tr><td>r<sub>2</sub></td><td>radius at the stem tip</td></tr><tr><td>s</td><td>position index</td></tr><tr><td>H</td><td>cone length</td></tr><tr><td>L</td><td>stem length</td></tr><tr><td>θ</td><td>displacement angle</td></tr><tr><td>F</td><td>force applied on the stem tip</td></tr><tr><td>ε</td><td>initial displacement angle at the stem bottom in the vertical plane</td></tr><tr><td>x</td><td>projection distance to the plane of the stem</td></tr><tr><td>y</td><td>deflection distance, vertical to the stem</td></tr><tr><td>e</td><td>final deflection distance, vertical to the stem</td></tr><tr><td>w</td><td>final displacement angle of the stem</td></tr><tr><td>d</td><td>density of the stem</td></tr><tr><td>J</td><td>objective function of the optimization problem</td></tr><tr><td>Q<sub>trunk</sub></td><td>total biomass of the trunk (trunk weight)</td></tr><tr><td>Q<sub>pith</sub></td><td>total biomass of the piths on the trunk</td></tr><tr><td>Q<sub>r</sub></td><td>total biomass of the rings surrounding the piths on the trunk</td></tr><tr><td>v</td><td>particle velocity</td></tr><tr><td>σ</td><td>inertia weight</td></tr><tr><td>c<sub>1</sub>, c<sub>2</sub>, c<sub>3</sub></td><td>acceleration coefficient</td></tr></table>
rand<sub>1</sub>, rand<sub>2</sub>, rand<sub>3</sub> uniformly distributed random value between 0 and 1

PI best position of the particle itself, recorded during the previous iterations

PG global best position in the population

PR best position of the particle selected in the population randomly

xp current position of the particle

Num number of particles in the population

Nd number of variables (dimension of the problem)