Application of the GreenLab Model to Simulate and Optimize Wood Production and Tree Stability : A Theoretical Study

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.


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.mechanoperception 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 andHoullier 1997, Woodcock andShier 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.

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, Green-Lab 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 predefined 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 Green-Lab 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.

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).

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.
where Q 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: 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, 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; 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 Eq. 4 illustrates the interaction between organogenesis (defining the evolution of organ number) and physiological processes in GreenLab.

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.
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 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 λ.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.
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 where λ 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 λ takes null value, only mode D 1 takes effect.Total biomass for cambial growth is distributed equally among metamers.If λ is equal to one, only mode D 2 takes effect.With this last mode, older metamers get more biomass.

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 .Fig. 3 illustrates the corresponding variable definitions.
When the stem leans with an angle (ε) with regard to the vertical, the equilibrium of moments of the force at position s from the stem base is given by: ) , 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 1r 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; θ is the displacement angle with regard to the stem without load; F is the force applied on the stem top; ε 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, We multiply Eq. 12 by dθ / ds and integrate it with θ that varies from 0 to ω.The result is obtained as Eq. 13.
where I 1 is the second moment of the circular cross section at the stem base πr 1 4 / 4. Substituting Eq. 14 to Eq. 13, we get Eq.15, The final shape of the stem is given by integrating Eq. 15 When θ 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.
where I 1 , r 1 , r 2 are determined by the summation of q i and q 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 ε is applied to the tree.

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 ring ).Stem is assumed to be conic as a first approximation.The part of cone with solid line is the tree stem.
where Q trunk represents the trunk weight, Q 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.
Note that S 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 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 populationbased 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 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.
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 jth variable of the ith particle at iteration k; PI ij is the best position of variable j recorded by the ith 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 kth iteration; σ 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.

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 cambial growth (S ring ) and of the coefficient (λ) on stem wood production and tree stability has been investigated.S ring and λ have distinct effects: S ring determines the total biomass allocated to wood formation at the whole-plant level while λ 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.

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  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 slowdown 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- mass production Q to total demand of plant D.
In addition, the sink strength for cambial growth S 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 ring , relative growth rates tend to be negative as shown in Fig. 6(b), and it implies that trees are under suppression.

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 ring 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.

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 ring is set to be the optimal value 0.1.DBH decreases with decreasing in λ, due to the dominant effect of mode D 1 .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 ring , where λ is set to be 1.As shown in Fig. 4(a), when S ring 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 ring 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 ring .Nevertheless, Fig. 8

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 λ, 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 λ, while the rotation angle at the

Discussion
In this paper, a population-based optimization algorithm has been applied to the functionalstructural 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 ring 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 ring , 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 ring 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 λ drives the proportion of upward to downward propagation.The value of λ 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 λ 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 λ 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 λ 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 λ. Considering the contradictory effects of λ 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 λ 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.

Fig. 1 .
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.
is the secondary moment of the circular cross section at position s from the stem bottom,

Fig. 2 .
Fig. 2. Illustration of wood distribution into metamerswith two modes: left one corresponding to mode D 1 as shown in Eq. 8 and right one corresponding to mode D 2 as shown in Eq. 9.

Fig. 3 .
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.

Fig. 4 .
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.

Fig. 5 .
Fig. 5. Ratio of biomass for cambial growth and for piths Q wood to tree biomass Q with respect to S ring .

Fig. 6 .
Fig. 6.Comparison of relative growth rates between wood growth (pith and rings) and tree growth.(a) S ring is set to the optimal value 0.1 (b) S ring = 1.
(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.

Fig. 10 .
Fig. 10.(a) Difference between diameters at the bottom and at the top of the stem with respect to λ.(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 2 (approximately 1000 plants/ha).

Fig. 11 .
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 ring 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.

Table 1 .
Parameter values of the GreenLab model for tree.
of organ o of physiological age p, initiated at growth cycle n o symbol representing organ types, e.g. a represents a leaf and i represents a pith P m maximum physiological age ta leaf functioning duration slw specific leaf weight Q ring (n) total biomass allocated to cambial growth at the whole-plant level at growth cycle n D ring (n) total demand for cambial growth at growth cycle n, at the whole-plant level D 1 , D 2 demand of cambial growth for biomass allocation to metamer S ring sink strength for cambial growth l p (n) length of the metamer of physiological age p, initiated at growth cycle n allocated to cambial growth of the metamer of chronogical age j and physiological age p, when the plant age is n NL p (n) number of living leaves above the metamer of physiological age p, initiated at growth cycle n λ coefficient determining the way the biomass assigned to cambial growth is allocated to 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) number biomass