Modelling the functional dependency between root and shoot compartments to predict the impact of the environment on the architecture of the whole plant. Methodology for model fitting on simulated data using Deep Learning techniques


 Tree structural and biomass growth studies mainly focus on the shoot compartment. Tree roots usually have to be taken apart due to the difficulties involved in measuring and observing this compartment, particularly root growth. In the context of climate change, the study of tree structural plasticity has become crucial and both shoot and root systems need to be considered simultaneously as they play a joint role in adapting traits to climate change (water availability for roots and light or carbon availability for shoots). We developed a botanically accurate whole-plant model and its simulator (RoCoCau) with a linkable external module (TOY) to represent shoot and root compartment dependencies and hence tree structural plasticity in different air and soil environments. This paper describes a new deep neural network calibration trained on simulated datasets computed from a set of more than 360 000 random TOY parameter values and random climate values. These datasets were used for training and for validation. For this purpose, we chose Voxnet, a convolutional neural network designed to classify 3D objects represented as a voxelized scene. We recommend further improvements for Voxnet inputs, outputs, and training. We were able to teach the network to predict the value of environment data well (mean error < 2%), and to predict the value of TOY parameters for plants under water stress conditions (mean error < 5% for all parameters), and for any environmental growing conditions (mean error < 20%).

as opposed to statistical-based or allometric-based modelling. This kind of modelling describes the processes that explain plant behaviour. It was first applied to cultivated agronomic and forest species and mostly focussed on biomass production. These studies included a very rough description of the structure (King 1991;Friend 2001;Buck-Sorlin 2013). On the other hand, many models represent plant structure. Concerning the aerial parts, de Reffye et al. (1988), Kurth (1994) and Smith et al. (2014) show the 3D shape of trees with varying degrees of botanical accuracy. Concerning roots, Pages (1999), Leitner et al. (2010) and Postma et al. (2017) provide models but these are mostly dedicated to monocotyledons. All these models make it possible to study interactions between one compartment and its environment (for instance, shoot vs. light or roots vs. soil). Studies on tree growth including both structure and biomass (functional-structural plant models, FSPMs) have mainly been carried out on shoot systems in agronomic and forest species (Vos et al. 2007(Vos et al. , 2010Sievanen et al. 2009;Reyes 2020;de Reffye et al. 2021). Particular attention has been paid to the influence of light in Buck-Sorlin et al. (2011) and hydraulics properties in Braghiere et al. (2020). The root system is usually taken apart due to the difficulties involved in observing and measuring it, especially root growth. However, climate change has made it crucial to study the structural plasticity of trees including root/shoot interactions. For this purpose, both shoot and root compartments must be considered simultaneously as they play collaborative roles with respect to climate traits (water availability for roots and light or carbon availability for shoots). Bornhofen and Lattaud (2007) proposed a theoretical model of whole trees to explain species evolution; however, it was only sketchy, far from the real plant architecture. LachenBruch and McCulloh (2014) and Ledo et al. (2018) also published wholes tree studies but these contained simplified description of tree structure. The question we aimed to answer with the present study is: what is plant plasticity with respect to biomass production and partitioning considering both shoot and root systems?
The production of biomass is the basis for the growth of both above-and below-ground tree axes. However, trees have the capacity to favour the growth of certain axes at the expense of others (e.g. trunk, stump), by partitioning biomass between the root and shoot compartments. Moreover, biomass partitioning has been shown to be driven by the resource that limits biomass production (Liu et al. 2012;Hertel et al. 2013). For example, under water stress, the tree allocates a larger share of the biomass produced to the root system, thus allowing more absorption of water. On the other hand, under light stress, the tree allocates a larger share of the biomass to the shoot system, thus allowing better absorption of light. The effects of drought are widely documented in the literature (Hommel et al. 2016;Ledo et al. 2018), light (Ballare and Pierik 2017) as are the effects of nutrient availability (Reynolds and D' Antonio 1996;Poorter et al. 2012) on biomass partitioning between root and shoot systems. For a review, see Poorter et al. (2012). More generally, here we consider the ability of the tree to adapt to resource availability by modulating biomass partitioning as defined in optimal partitioning theory (Niklas and Enquist 2002;Reich 2002;Colter Burkes et al. 2003): resources are allocated to both compartments in order to balance their contribution capabilities.
First we present a model and its simulator which merges the DigR root model (Barczi et al. 2018) and AmapSim Caulinar model (Barczi et al. 2008) to produce a whole-plant structural model we named RoCoCau (for Root, Collar, Caulinar). Since the architecture of both systems and the visible traits that express this partitioning are indispensable, RoCoCau was designed to mimic whole-plant growth with very good botanical accuracy.
Second, we present TOY, an external model and its simulator that is linked to RoCoCau. TOY represents shoot and root compartment dependencies and hence tree structural plasticity in different air and soil environments. The TOY model aims to describe links between the root and shoot systems for both the production and sharing of biomass when the plant is faced with a variable environment while simultaneously describing both compartments during growth with good botanical accuracy.
Before model validation, a recurrent question for any model is fitting. RoCoCau parameters are calibrated using traits directly measured in the field. The problem of model fitting was solved in previous studies (Barczi et al. 2008(Barczi et al. , 2018. On the other hand, most of the TOY model parameters are hidden, meaning they cannot be measured directly. In this paper, we calibrate 25 numerical values for the TOY model, which is a huge number for common fitting methods (least square, simulated annealing, genetic algorithm, etc.). The solution space is very wide and standard methods may be captured in local minima. Some model inversions have already been tried (Matsunaga et al. 2001;Higgins et al. 2010;de Reffye et al. 2016), but they targeted a small number of hidden parameters that cannot be directly measured on the plant. In this paper, we describe our experimental calibration using a convolutional neural network. This technique is mainly used for object detection or classification from numerical images, but is sometimes used in a generative way to predict parameter values for model calibration (Danson et al. 2003;Bacour et al. 2006;Su et al. 2015;Song and Xia 2016;Qi et al. 2016).
Another difficulty when fitting the TOY model is the lack of experimental data, especially on root system architecture. The architectural traits of a root system can only be measured through destructive methods involving the uprooting of the tree, which for logistic reasons, hinders the creation of data sets that are large enough to allow the training of a neural network using real data. To validate our claim that RoCoCau can produce simulated trees that are close to reality, we use standard traits (shoot length/root growth speed, branching density, geometry, etc.) to calibrate RoCoCau on particular species and we then calibrate the TOY model using simulated data sets to train a neural network.

Model
The aim of RoCoCau+TOY FSPM is to predict and simulate the impact of the environment, particularly resource availability, on wholeplant architectural development. The RoCoCau structural model enables accurate simulation of whole-plant architectural development, based on a few architectural measurements. However, RoCoCau does not explicitly model the functional processes that drive architecture development. These include both endogenous (genetic) and exogenous (environment) functional processes and prevent RoCoCau from simulating the impact of the environment on whole-plant architectural development. We designed the TOY model to (i) disentangle the effects of genotype and environment on whole-plant architectural Methodology for model fitting on simulated data • 3 development, (ii) predict the impact of environmental constraints, especially water and light availability, on biomass partitioning between the root and shoot compartments, (iii) plug in to RoCoCau, allowing RoCoCau+TOY FSPM to simulate the development of whole-plant architecture with respect to these functional constraints.

RoCoCau-a structural model to simulate whole-plant growth.
Field experiments and observations led us to develop a whole-plant model and its simulator named RoCoCau. This pure structural model is the result of merging two existing structural models; a model for the shoot system, AmapSim (Barczi et al. 2008) and a model for the root system, DigR (Barczi et al. 2018). Whole-plant architecture development is modelled as a set of virtual shoot and root meristems, whose development is driven by stochastic processes that represent apical growth (axis elongation), lateral growth (branching) and death. For the shoot system, the value of the parameters of these processes varies with the physiological age (partially expressed by the term 'vigour') of each virtual shoot apical or lateral meristem. The shoot system is described as a topology of shoots belonging to different categories. The physiological age of each shoot meristem evolves along a reference axis defining the different possible values of physiological age according to an oriented finite automaton whose law of occupation and transition follows a semi-Markovian function (Barczi et al. 2008;de Reffye et al. 2021). The root system is described as a topology of roots belonging to different categories in the same way as the shoot. The vigour of a root meristem depends on the type of root and on the current length of the root (Barczi et al. 2018). The geometry of the resulting axes is then modelled using simple descriptive functions. The stochasticity of the processes modelling the development and the differentiation of each axis implies the uniqueness of each simulated axis. Nevertheless, it is possible to classify axes in distinct structural categories according to their morphological specificity and to define the way these structural categories are combined into shoot and root branching systems (i.e. architectural units; Barthélémy and Caraglio 2007). Both models produce botanically accurate shapes obtained by calibrating their parameters using measurements taken on real plants. The models were tested and validated on many plant species (Barczi et al. 2008(Barczi et al. , 2018. The simulator of RoCoCau combines the AmapSim and DigR simulators in a single program in which all shoot and root meristems are synchronously simulated with a software architecture that allows external models that can interact with the default simulation kernel to be plugged in (Barczi et al. 2008).

RoCoCau+TOY, an FSPM to simulate the impact of environmental constraints on the development of whole-plant architecture.
RoCoCau is a purely structural model which incorporates both endogenous (i.e. growth and branching patterns according to structural categories of axes) and exogenous effects on the growth of the whole plant. To disentangle the two effects, we developed TOY, an extra model that represents (i) biomass production according to current plant structure and environmental constraints, (ii) biomass partitioning and its influence on plant growth and architectural development. Within TOY, functional dependency is modelled only at the compartment level, i.e. between the root and shoot compartments. Within each compartment, the role of the axes can be associated with three main functions: structuring, exploration and exploitation. This is the reason why, in TOY, we arbitrarily created six independent functional categories for the axis (three for the shoot compartment and three for the root compartment) according to the function to which they contribute the most. For instance, for the shoot system, the first functional category, grouping the axes that mostly contribute to the structuring function includes the trunk and main forking branches; the second functional category, grouping the axes that mostly contribute to exploration includes the other branches; the third functional category, grouping the axes that mostly contribute to exploitation includes the short shoots. It is important to note that the functional categories can include several structural categories. It is up to the modeller to define the architectural units of the shoot and root systems (i.e. structural categories and their relative organization). For RoCoCau settings and for TOY settings, it is also up to the modeller to define the content of each of the six functional categories in terms of structure. The flexibility of this approach makes it possible simulate a wide range of architectures considering both structure and function. Figure 1 shows a simple example of a shoot/root architectural unit and a shoot growth loop. In the rest of this paper, it is important to clearly distinguish between the structural axis categories (a classification defined in RoCoCau) and the functional axis categories (a classification defined in TOY).
The software implementation of RoCoCau enables dynamic linking of plug-ins that contain additional knowledge that will interact with the default RoCoCau algorithm. TOY plug-in processes take place in two different modes (Fig. 2). (i) At the end of each simulation time step, RoCoCau stops simulating growth and the TOY plug-in takes over. Knowing the actual shape of the tree and the quality of the environment (water and light), TOY computes biomass production, partitioning and the resulting growth attenuation factors. It then hands back to RoCoCau. (ii) In the following time step, RoCoCau uses the TOY actual attenuation factors to moderate apical growth, branching and death of meristems depending on which compartment and which functional category they belong to.

Description of the TOY model.
The growth of the whole plant is the result of endogenous processes, driven by the expression of the plant's genotype, and of exogenous processes, driven by the interaction between the plant and its environment. In simulations, each time step usually represents a year of growth. Exogenous processes are updated at the beginning of each time step, and endogenous processes are called on throughout the duration of the time step for the simulation of the growth of each axis. In the TOY model, the environment is defined as a set of three parameters, representing light (useable light energy, Light (w·m −2 )), inorganic carbon (concentration of atmospheric carbon, C air ) and water (soil relative humidity by volume, Hv soil) resources.
Exogenous growth processes, namely absorption of resources, sugar production and partitioning between the root and shoot compartments, are modelled by a set of equations using four parameters: k, Light max , R ‖ and RAL. The light extinction coefficient k, and the maximum useable light energy Light max (w·m −2 ) are parameters of the light absorption equation. As inspired by Nye (1973) and Jungk (2001), the root absorbing length (RAL) represents the length of the apical absorbing part of the roots and the R abs radius of the volume of soil investigated. We are aware that our representation of the soil and of water transfer are very rough compared to those of Doussan et al. (2006) and Garrigues et al. (2006), but our aim is to describe the different processes that occur in the plant at the same level of detail. Modelling exogenous processes also requires whole-plant variables that are provided by RoCoCau: leaf surface area (Sleaves), which represents the sum of the surface area of all leaves, leaf area index LAI, the radius of the root absorbing soil R abs and RAL.
Endogenous processes regulate the growth, branching and death of the axes. The processes are specific to each axis functional category and are modelled by a threshold function with three parameters: S m the mortality threshold, S r the branching threshold and A l the apical growth factor. Since the TOY model distinguishes six-axis functional categories, the total number of parameters used by the TOY model is 25 (7 + 6 * 3).
The TOY model functions in four successive stages: (i) estimating sugar production according to light, carbon and water resources, under actual and optimal environmental conditions, (ii) estimating actual and optimal sugar partitioning, i.e. estimating the quantity of sugar allocated to the root and shoot compartments, under actual and optimal environmental conditions, (iii) comparing the actual and optimal sugar partitioning in each compartment to estimate the distance between actual and optimal growth conditions (growth factor), (iv) for each axis, the functional category within each compartment, determining the response of the axes to the distance between actual and optimal growth conditions (growth factor), i.e. determining the impacts of a reduction in the quantity of sugar allocated, under actual growth conditions, on the death, branching and growth processes of the axis. We now detail each of the four steps.
(i) The production of sugar by photosynthesis is calculated as a function of the quantities of water, carbon and light absorbed. Carbon absorption is modelled as the product of leaf surface and the concentration of carbon in the air Carbon(t) = S leaves (t) × C air; water absorption is modelled as the product of the volume of the absorbing parts of the roots and soil relative humidity Water(t) = V(t) × Hv soil, with V(t) = π × R abs 2 × ∑ allroots min(Rootlength(t), RAL); light absorption is modelled as the product of available light energy Light(t) (w·m −2 ) and photosynthetic efficiency, which can be represented by a decreasing exponential as a function of the leaf area index LAI(t) and the light extinction coefficient k. Sugar production is assumed to be controlled by the limiting resource (the minimum between absorbed water and carbon quantities).
(1) Optimal environmental conditions are defined as soil relative humidity equal to 1, i.e. Water max (t) = 1, and light absorption equal to maximum useable light energy Light max (w.m −2 ). According to equation (1), optimal sugar production Sugar max (t) is then: The plant as a branched system of axes belonging to different structural categories (e.g. trunk, branches and short axes; taproot, lateral roots and fine roots), that can be grouped in six functional categories, three for each compartment, e.g. structuring category (green: shoot system, dark red: root system); blue, purple: exploration category; orange, pink: exploitation category. (B) Shoot axis development process implemented by RoCoCau (simplified). Φage stands for the physiological age of the virtual meristem.
Methodology for model fitting on simulated data • 5 (ii) According to the optimal partitioning theory, the partitioning of sugar production between the root and shoot compartments is assumed to be inverse to the contribution to production. Sugar partitioning between root and shoot is then calculated as follows: According to equation (2) optimal partitioning of sugar production is: (iii) For each compartment, the distance between actual and optimal growth conditions is estimated by calculating its growth factor; i.e. the ratio of the quantity of sugar allocated under actual growth conditions to the quantity of sugar allocated under optimal growth conditions.
(iv) For each axis functional category in each compartment, the response of an axis to the distance between actual and optimal growth conditions, i.e. the impacts of the value of the growth factor on axis death, branching and growth processes, are modelled by a threshold and attenuation function (TAF; Fig. 3). When the value of the growth factor is equal to 1, the growth of the axis is not attenuated. When the value of the growth factor is less than 1 and more than the branching threshold S r , the growth rate of the axis is attenuated by a coefficient that decreases linearly with the value of the growth factor, according to the values of S r and the apical growth factor A l. Then, when the value of the growth factor falls below the branching threshold S r , the axis stops branching, and the attenuation coefficient decreases linearly with the value of the growth factor, according to the values of S r , A l and the mortality threshold S m. Finally, when the value of the growth factor falls below the value of the mortality threshold S m, the axis dies.

Model calibration
The main difficulty in disentangling exogenous and endogenous growth processes is providing an adequate description of endogenous growth processes, and more specifically in estimating parameter values for endogenous growth equations. Indeed, endogenous effects on plant architecture are (i) inherently more difficult to observe and measure and (ii) perceivable only by monitoring the whole plant's architectural development over time. Standard calibration of the 18 TOY hidden parameter values for endogenous growth (Table 1), relying on simple relationships between parameter values and measurement of specific architectural traits is therefore impossible. In this section we present Figure 2. Schematic representation of the functioning of RoCoCau+TOY FSPM: the colours refer to the specificity of the modelling process; red for processes that occur at whole-plant level, blue for processes specific to each compartment (root and shoot) and green for processes specific to each development based on the filling of the zones of a segmented 3D image with the voxelized density maps of each successive plant structure.
a method for the calibration of the 18 TOY hidden parameters for endogenous growth.
Since the model is still in the early stages of development, the aim of this section is rather to address the main difficulties encountered in calibration than to provide complete calibration of the model. As the parameters for exogenous growth can be calibrated using standard calibration techniques, we do not carry out this task here, but this issue is addressed in the discussion section.

Model calibration on simulated data using Deep Learning techniques.
The most efficient way to calibrate TOY hidden parameters, as performed in Cournede et al. (2011) andde Reffye et al. (2018), is model inversion, so as to define a method enabling assessment of the input parameter values from the simulated plant architectural development. Generative deep neural networks can learn the intricate relationship between each parameter value and plant architecture development over time. Provided with an input consisting of a representation of plant architecture development over time, they can assess TOY model parameter values.
Generative deep neural networks can be trained on a simulated training data set, with each example of the training data set consisting of a simulation of the plant's architectural development by the RoCoCau+TOY FSPM, labelled with the values of TOY model parameters.
The two main choices to be made for this method of calibration are the choice of a suitable representation for the simulations of the plant's architectural development, and the choice of a suitable architecture for the deep generative neural network. Indeed, the simulations of a plant's architectural development must be in a format that can be processed by a deep neural network. The chosen representation must introduce as little bias as possible and must be obtained from a simulation as fast and easily as possible to be able to build a large training data set. The chosen network architecture must be able to learn the link between the representation of plant architectural development and TOY parameter values rapidly and with the best possible accuracy.  Methodology for model fitting on simulated data • 7

Choosing to represent plant architecture development over time as a voxelized density map that can be processed by the 3D deep convolutional neural network VoxNet.
One of the outputs of a plant's architectural development as simulated by the RoCoCau software is a time series of three-dimensional mock-ups. However, it is possible to get rid of the time dimension entirely by representing plant architecture development as a scene, by including the successive structures of the plant in the same 3D image. Point or voxel clouds are already used in LiDAR data processing (Wandinger 2005;Côté et al. 2009;Mak and Hu 2015;Ferrara et al. 2018;Fan et al. 2020;Estornell et al. 2021;Mayra et al. 2021), and are a well-documented and efficient way to represent a three-dimensional plant structure. Voxel clouds are equivalent to point clouds, except that the point coordinates are fixed and equidistant from one another, thus avoiding potential bias and making it possible to define the necessary neighbourhoods to process convolution in neural networks. Moreover, non-binary voxels can have a value. This value can be used to describe the amount of plant biomass in the voxel's volume. This makes it possible to represent plant structure as a distribution of biomass in a 3D space, i.e. as a voxelized biomass density map.
Since it is relatively easy and quick to transform RoCoCau simulations into voxelized density maps, this representation can create a 360 000 simulation training data set in less than 2 days using a standard PC. Moreover, even if the voxels are large, the information in each voxel concerning the amount of biomass it contains makes it possible to obtain a relatively accurate representation of plant structure.
The VoxNet network presented in Maturana and Scherer (2015) is a deep 3D convolutional neural network, and is promising for 3D object recognition tasks (Maturana and Scherer 2015;Liu et al. 2021). The use of convolution by VoxNet makes it possible to rapidly process 3D images such as voxelized density maps, which could be too big for non-convolutional networks. VoxNet was originally developed to predict the class label of a 3D object from its occupancy grid, but it can be adapted to generate parameter values from a voxelized density map representing plant architectural development.
The VoxNet network is also particularly well suited for processing point cloud (for instance, LiDAR) data (Maturana and Scherer 2015) after conversion to voxel data. This makes it possible to use it not only for calibration on simulated data, but also for calibration on real data acquired with LiDAR.

Pre-processing the simulation of architecture development by RoCoCau, to be used as inputs for VoxNet
Representing plant architectural development as a voxelized density map requires the definition of a unique pre-processing procedure that can transform a simulation by the RoCoCau software into a voxelized density map (Fig. 4). In this study, the simulations provided by the RoCoCau software are performed in 10 time steps, representing 10 successive plant structures. To be compatible with the VoxNet network input format, the size of the 3D image containing the voxelized density map is fixed at 30 × 30 × 30 voxels. For each simulation, the pre-processing procedure relies on (i) the voxelization of the 10 successive structures of the RoCoCau simulation, and (ii) the segmentation of the 3D image into 10 timespecific zones to host these structures at 10 successive ages.
Since the TOY model relies on absorbing zones (leaves and absorbing roots) for biomass production and on woody parts for growth control, TOY calibration takes advantage of separating these two traits to make calibration more efficient.
(i) Each successive structure of the plant is separated into two, providing two substructures for the absorbing parts (leaves and absorbing roots) and for the woody parts (branches and roots) of the plant. Each substructure is then placed in a 3D space made up of voxels of whose sides measure 20 cm. Each voxel is assigned a value according to the amount of plant biomass within its volume. Two voxels within the 3D space accommodating the woody substructure are attributed values representing water and light availability. This voxelization procedure makes it possible to obtain a voxelized density map of a single plant structure.
(ii) The 3D image is separated into two equal halves, on the left (resp. right) for the accommodation of absorbing (resp. woody) substructures. Each half is then divided into 10 time zones of appropriate size to accommodate each of the 20 voxelized substructures. Filling the time zones of this segmented image with the previously voxelized substructures, with respect to their chronological order, enables production of a voxelized density map of plant architectural development.

Training VoxNet to predict TOY hidden parameter values
The calibration of the TOY model is not a classification problem (Fig. 5A).
As it was originally developed to classify 3D objects, the architecture of VoxNet was adapted to predict TOY hidden parameter values. Adaptation concerned the fully connected layer and the loss function, which is set as the root mean square error (RMSE). To optimize speed and quality and thanks to the independence of the functional categories of the six axes defined in the TOY model, VoxNet was trained six times, each training run for one of the functional categories and aimed to predict their three parameter values. For this purpose, the size of the fully connected layer was reduced to three. The final result consists of six trained VoxNets, each having a set of weights computed to predict the three parameters of each of the functional categories.
According to the recommendation for the use of VoxNet (Maturana and Scherer 2015), the size of the training data set was set to 300 000 examples belonging to the same species which means that RoCoCau parameters were fixed to describe a particular species. The creation of the training data set followed a four-step procedure: (i) generating random values for TOY parameters; (ii) simulating plant architecture development over 10 time steps with RoCoCau+TOY FSPM using the TOY generated parameter values; (iii) pre-processing simulations of plant architectural development following the pre-processing procedure described above and labelling each voxelized density map with the three parameter values that VoxNet will have to predict; (iv) converting the resulting examples into a training-friendly format (.npy) (Maturana and Scherer 2015).
Since a single training run allows the prediction of three parameter values, the prediction of TOY 18 hidden parameter values requires six training runs. Beyond improving the quality and increasing the speed of each individual training run, this division also enables fine-tuning. Once a training has been performed, e.g. for the first shoot functional category, it is indeed possible to speed up the remaining trainings, e.g. for the second and third shoot functional categories, by using the weights of the previously trained VoxNet as initial values.
A priori, since TOY parameters depict plant response to environmental stress, their values should be exceedingly difficult to predict from a representation of a plant architecture that has developed in a stress-free environment. To avoid training failure, we isolated the subdata set SubDBStress containing only examples in which the plant architecture developed under light and water stresses (reduction of at least 25 % in light availability and of at least 25 % in water availability compared to the optimum). Once a training run has been performed on SubDBStress, it is possible to avoid the failure of the training on Methodology for model fitting on simulated data • 9 the whole data set by using the weights of the VoxNet trained on SubDBStress for initialization.
Finally, the training procedure for the calibration of TOY 18 hidden parameters with fine-tuning followed three steps: (i) training two VoxNets to predict the three parameter values for the first shoot (resp. root) functional category on SubDBStress, (ii) training four VoxNets to predict the three parameter values of the second and third shoot (resp. root) functional categories on SubDBStress fine-tuned with the weights of the training on the first shoot (resp. root) functional category on SubDBStress, (iii) training six VoxNets to predict the three parameter values of each functional category on the whole data set, fine-tuned with the weights of the training on SubDBStress.

Testing the quality of the prediction by trained VoxNet
In each training run, changes in the value of the training loss function (RMSE) were monitored to ensure the effective learning of VoxNet (Fig. 5B). Same monitoring was also carried out on the validation data set to avoid overfitting.
The quality of the training is tested on a test data set containing 400 test examples. For each example of the test data set, the input of the example, i.e. the voxelized density map is used as an input to the trained VoxNet. The trained VoxNet then delivers a prediction for the values of the chosen set of three hidden parameters it has been trained to predict, e.g. the death threshold, the branching threshold and the apical growth factor (see Table 1). Its prediction is then compared with the output of the example, containing the true values of the three hidden parameters. The comparison is performed through an estimation of the RMSE, and the relative absolute percentage error (RApE) between each hidden parameter true and predicted values.
To test the quality of VoxNet predictions under contrasted environmental conditions, the test data set was divided into four sub-data sets, according to the severity of light and water stress under which the plant architecture of the example developed. Each training run was tested on the whole test data set as well as on each of the four test sub-data sets (see Results)

RoCoCau+TOY FSPM use and analysis in silico
Depending on the modeller's knowledge of the tree species to be modelled, it is up to him/her to define the sensitivity of each of the six functional categories proposed by the TOY model. Combining the sensitivities of all six-axis functional categories defines the sensitivity of the whole plant. It will be recalled that we set the number of functional categories based on the three main functional roles of the shoot or root axes, i.e. structuring, exploration and exploitation. Once the RoCoCau model was set up to a particular species and the TOY model was set up and its simulator running, two kinds of tests were carried out to check its ability to mimic what is known about whole-plant functioning by: (i) testing the phenotypic plasticity of virtual plant species while varying the environment, (ii) testing the ability of the model to represent different plant strategies while varying the endogenous parameters of the six functional categories.

Predicting phenotypic plasticity under varying environmental conditions
By varying the TOY parameter values, we tested the capacity of the TOY model to express plant plasticity in different environments. For instance, we programmed a plant with a second functional branch category (lateral in this case) with highly sensitive apical growth (A l = 0) while all the other shoot functional categories are less sensitive (A l = 1). This results in strong apical dominance under stress. In the below-ground compartment, the third root category (fine roots in this case) is set sensitive to death (Sm = 0.3) while the other categories are not (Sm = 0.1). The first functional root category (taproot in this case) has highly sensitive apical growth (A l = 0) while other functional root categories are less sensitive (A l = 1). This results in the plant being sensitive to drought. We simulated the growth of this theoretical plant using combinations of two light conditions (maximum and half maximum) and two water conditions (maximum and half maximum). Figure 6 shows the result of the four simulations. At first sight the differences between the shape of the same plant in different environments are clearly visible. Taking the plant with no stress as reference (lower right), we see that the plant facing light and water stress (upper left) is globally weaker and did not even branch in Years 8 and 9. Both root and shoot systems are small. We see that plants under light stress or water stress are smaller than the reference plant. We also see that the root/shoot ratio of plants under water stress is opposite to the ratio of the plant facing light stress. If we look at the numbers, the impression we got at first sight is confirmed. Plant under both stresses have fewer and shorter axes than the plant in optimal environment. The three other plants have the same number of branches but the plant under water stress has twice as many fine roots as the plant under light stress, which has shorter side branches.

Generic modelling of genotype, defining plant sensitivity to environmental constraints
Three parameters of each of the six functional axis categories control its death (Sm), branching (Sr) and apical growth ( A l) depending on the current value of the growth factor of the corresponding compartment. Figure 7 shows four examples of the possible use of these three parameter values in the second category of branches (lateral branches). Parameter A l controls apical growth

Methodology for model fitting on simulated data • 11
attenuation between the branching and the non-branching zone of the growth factor. The higher the value of A l, the less sensitive it is to values of the growth factor greater than S r and the more sensitive it is to values of the growth factor between S m and S r (Fig. 7, upper two cases). The lower the value of A l, the more sensitive it is to values of the growth factor greater than S r and the less sensitive it is to values of the growth factor between S m and S r (Fig. 7, lower two cases). Parameter S r controls branching. The lower the value of S r , the larger the number of branches (Fig. 7, case in the upper left with many short axes). The higher the value of S r , the fewer branches (other cases with a few short axes). Parameter S m controls axis mortality. The higher the value of S m, the fewer the living axes in this category (Fig. 7, lower right case). The lower the value of S m, the larger the number of living axes (Fig.  7, other cases).
It is easy to grasp the huge number of possible combinations available to define different genotypes and plant sensitivity to environmental constraints by varying the 6 × 3 parameter values. The modeller can define them manually based on his/her knowledge. Another possible way would be to use a calibration tool able to predict the correct values of the parameters of the model based on measurements.

Calibration results
In this section, we present the results of implementing the calibration method (see section 2.2). The calibration method relies on training the VoxNet neural network on simulated data. The results of calibrating the TOY model are thus to be understood as predictions made by the trained VoxNet neural networks. We also assessed the quality of VoxNet predictions using simulated data (see section 2.2.4).

Calibration of the three hidden parameters of the TOY first functional shoot category
The aim of the calibration method was to assess TOY hidden parameter values. The TOY hidden parameters are the parameters of endogenous growth processes, i.e. the death threshold S m, the branching threshold S r and the apical growth factor A l for each functional axis category (see Table 1). Figure 8 shows  (b) The quality of the predictions by VoxNet trained on the whole training data set fine-tuned with the weights of the training on SubDBStress was significantly better, with a mean RApE on the whole

Methodology for model fitting on simulated data • 13
test data set of 13 % for A l, 14 % for S r and 22 % for S m. The predictions made by VoxNet trained on the whole data set for example, under water stress alone or under both water and light stress were still slightly more accurate than the predictions made for examples under light stress or no stress at all. However, they were significantly more accurate than the predictions made by the previous VoxNet trained on SubDBStress. The overall improvement in the quality of predictions given by the VoxNet trained on the whole data set was driven by the notable improvement in the quality of the predictions for the test examples under light stress alone and under no stress. Indeed, in comparison with the previous training, the median RApE for the test examples under light stress alone (resp. no stress) decreased significantly from 28.5 to 7.5 % (resp. from 37.6 to 10.0 %) for A l, from 19.2 to 9.0 % (resp. from 35.8 to 13.5 %) for S r and from 32.0 to 18.1 % (resp. from 46.4 to 17.5 %) for S m. Since the values of environmental parameters are included in some voxels of the voxelized representation of plant architectural development, the only explanation for this reduction in accuracy is that the information contained in those voxels is partially lost through convolution.
The results of both trainings share one similarity: the predictions of the values of the branching threshold S r and the apical growth factor A l were significantly more accurate than the prediction of the value of the mortality threshold S m. This may be due to the fact that the three parameters do not modulate plant responses to the same degree of severity of the stress; A l and S r modulate plant responses to mild environmental stresses, whereas S m modulates plant responses to more severe stresses. Since severe stresses, modelled by low values of the growth factor, are less frequent during and among the simulations that comprise the training data set, it is probably more difficult for VoxNet to learn how to predict the value of the mortality threshold S m.
A particularly interesting result is the prediction that VoxNet is able to provide on the examples under limited stress. Even though the prediction is of low quality, the ability of the network to predict plant responses to environmental stresses based on an example in which the plant architecture developed in the absence of stress is quite remarkable. This result requires more detailed investigation. In particular, the definition of stress used to classify the examples is rather vague (decrease of at least 25 % in water (resp. light) availability compared to the optimum). It would be interesting to refine this definition, to be able to quantify the extent of VoxNet's ability to predict plant responses to environmental constraints based on the simulation of a plant facing few or no environmental constraints.

Calibration method inversion.
Endogenous growth processes are assumed to be identical amongst all individuals of the same plant species. As such, the values of the 18 hidden parameters of the TOY model for endogenous growth are assumed to be associated with the plant species. For a given plant species, once these 18 values have been assessed, phenotypic plasticity between individuals is assumed to be modulated only by the values of the environmental parameters and exogenous growth processes.
Since the environment is measurable, the values of the environmental parameters can be assessed using standard calibration methods. However, since it is possible to construct a large data set containing architectural development simulations of an individual of a given species growing in different environmental constraints, it is also possible to train the VoxNet neural network to predict the environmental conditions under which the individual has developed.
For the purpose of inverting the calibration method, further training of VoxNet was undertaken. Eighteen hidden TOY model parameter values were fixed, as well as the three values of the parameters for exogenous growth processes. The fixed values of these parameters were carefully chosen, such that the resulting simulated individual expresses phenotypic plasticity, which, in turn, modifies its architectural development under different environmental constraints. Since most of the parameter values were fixed and only two values needed to be predicted, it was possible to reduce the size of the training data set to 20 000 examples. All the examples were constructed with the RoCoCau+TOY simulations of individuals under varying environmental constraints (for each simulation, the water and light availability values were chosen randomly between 0.5 and 1). The simulations were represented as voxelized density maps but without the two voxels in each structure representing water and light availability and labelled with the true values of the two parameters for water and light availability (see section 3). Finally, VoxNet architecture was slightly modified by reducing the size of the fully connected layer to 2, thus making it possible to predict two instead of three parameter values.

Results.
The predictions given by the trained VoxNet shown in Fig. 9 are fully accurate for every example of the test data set, with a mean RApE of 1.5 % for light availability, and of 2 % for water availability. The results of the inversion of the calibration method are therefore conclusive: they demonstrate VoxNet's ability to accurately predict the impact of water and light availability on the architectural development of a given individual of a given species.

DISCUSSION
In this section, we first discuss the main limitations and possible improvements of the TOY model. As the model is still in its early stages of development, the discussion is not intended to be exhaustive but rather to identify avenues for immediate improvement. We then discuss the method of calibration, mostly focusing on possible improvements to the calibration of current and future versions of the TOY model.
Throughout the conception of the calibration method for the TOY model, particular emphasis was placed on (i) minimizing the quantity of laborious tree architectural measurements to be performed for the calibration and (ii) separating the calibration of parameter values for the environment, and for exogenous and endogenous processes, so that RoCoCau+TOY FSPM can be calibrated to simulate phenotypic plasticity (Fig. 4). With respect to these two constraints, we propose a procedure for the calibration and use of RoCoCau+TOY FSPM based on measurements on a single in vivo tree.

Limitations and possible improvements of the TOY model
The main objective of this study was to build a modelling framework able to (i) accurately simulate whole tree growth, with particular attention paid to the interdependencies between the root and shoot compartments and (ii) to untangle genotype and environment effects on architectural development to predict phenotypic plasticity. Since such a modelling framework will be useful for a wide range of tree species when linked to RoCoCau, the TOY model can be applied to any type of tree and its set-up is therefore mainly focused on its genericity and precision rather than on its physiological realism. We also discuss the genericity of RoCoCau+TOY FSPM in this section. The first limitation of RoCoCau+TOY FSPM is the independence of the six-axis functional categories. It will be recalled that we chose six functional categories to express the structuring, exploration and exploitation processes that occur during plant growth. However, more structural categories may be needed to describe the architectural unit of a species. For each of the six functional categories the RoCoCau simulator modulates axes development independently, giving RoCoCau+TOY FSPM high genericity. However, such independence is not realistic, since the growth of all the axes is subject to biomass conservation, i.e. the cumulated biomass of all the new axes and of all the new portions of axes at a particular time step is equal to the amount of biomass allocated to tree growth. For RoCoCau+TOY FSPM to account for biomass conservation, the TOY model thus needs to be refined. Biomass sharing between functional axis categories can be explicitly written to express the strategy used by the plant species to face stresses over time (for instance, pioneer trees that mainly invest in trunk growth, whereas shade tolerant trees prefer to invest in lateral branches; Cai et al. 2013;Coll et al. 2008;Delagrange et al. 2008;Henry et al. 2010). It could also be improved by adding equations to estimate the amounts of biomass allocated for other essential functions of the tree, i.e. reserve, maintenance, reproduction and defence.
The second limitation concerns the branching process, which is mainly managed by the RoCoCau software. Indeed, axes of the first and the second categories can branch into new axes at each time step, the probabilities being fixed by the RoCoCau software. Within RoCoCau+TOY FSPM, the TOY model gives a binary modulation to those branching probabilities. It sets them to zero when the value of the growth factor is below the branching threshold. A first refinement of TOY model could consist in implementing continuous modulation of the branching probabilities, e.g. a linear decrease in the branching probability following the decrease in the growth factor. A more important refinement however, concerns the branching of the axes of the first functional category, e.g. the trunk and taproot. Indeed, axes of the first functional category can branch into axes of the second or the third category, e.g. branches or short shoots, lateral roots and fine roots, with different degrees of probability. It is likely that trees can modulate the branching process towards the second or the third category, towards exploration or exploitation, depending on the environment. Given that, the two probabilities are fixed by the RoCoCau software, and simultaneously set to zero when the growth factors value drops below the branching threshold, the RoCoCau+TOY FSPM is currently not able to account for new axis strategies, which reduces its genericity. The FSPM could be improved by adding new equations to the TOY model to allow it to modulate the two branching probabilities of the first category based on the environmental constraints.
The third limitation concerns modelling the environment, in particular the soil compartment, which is rudimentary in the current version of the FSPM. Indeed, the impact of the soil compartment on tree growth, and biomass partitioning is currently limited to water availability, which is itself considered homogeneous and is represented by the constant parameter Hv soil. The first refinement to the TOY model therefore concerns spatialization of soil water availability, e.g. by representing the soil compartment as a set of layers or a 3D grid of cells containing different water availabilities, as frequently observed in the field (Manoli et al. 2014;Braghiere et al. 2020).
Finally, as shown in Reynolds and D' Antonio (1996) and Poorter et al. (2012), nutrient availability is a determining parameter of the optimal partitioning theory (Reich 2002). It has a significant impact on tree biomass partitioning between the root and shoot compartments and on growth processes. Another refinement of the TOY model should therefore be including nutrient availability, as in McMurtrie and Landsberg (1992), D' Antonio (1996), andLacander et al. (2021).

Analysis
We chose a neural network suited to our specific task and improved it. Each plant growth simulation produces a 3D voxel image that contains 10 growth time steps split between the woody and absorbing parts. It is currently possible to predict 6 × 3 hidden TOY parameter values with acceptable accuracy after training on a simulated data set, which is a very promising result. As the input of the network consists in a voxelized representation of the plant, it is independent of the model to be inverted. It is thus possible to keep the same neural network to calibrate future improved versions of TOY. VoxNet was chosen for calibration because of its 3D input and because it is simple to adapt for generative purposes. It is clear that including time series in a single 3D image may not be the most accurate way to calibrate a dynamic model, not to mention including light and water availability in some voxels. It would thus be interesting to test recurrent networks whose structure natively includes this time feature. Another possible improvement would be to compute the loss according to the distance between the 3D images obtained from the target parameter values and from the predicted value rather than on the distance between the original and predicted parameter values. This will require calibrating all parameters in the same model, which was not possible for us due to hardware restrictions. Calibration will also be improved and tested further by our team, especially on yearly variations in the environment and using multiresolution 3D inputs rather than a single 3D input.

Future outlook
Although time series of environmental conditions can be used as inputs for the TOY model, thus accounting for inter-annual variability, the current architecture of the VoxNet neural network only allows the prediction of a single value for the parameters that influence light and water availability. As VoxNet is already able to predict the impact of a permanent stress on individual plant architectural development very accurately (see section 4.2), it would be interesting to introduce recurrence into its architecture (Medsker and Jain 2000;Liang and Hu 2015;Zuo et al. 2015), to enable it to predict time series of changing environmental parameter values, and also to test if it is able to predict the impact of episodic stresses. For instance, what would the difference in the final plant architecture be if it were affected by a single year of drought at the beginning or at the end of plant life (Backhaus et al. 2014).
Another improvement would be increasing the size of the input layer of the network to be able to process bigger or more complex trees. In this study, only the output of the network was adapted to make it generative. Also, the accuracy of the prediction was checked using the values of the predicted parameters. It would be interesting to compare the simulation run with the predicted values with the original tree. A more complex way would be to replace a single VoxNet with a single input 3D voxel image containing 10 growth time steps with a set of 10 VoxNets, each of which has an input of the corresponding shape at this time step, and each of the 10 VoxNets would be fully connected at the last layer. A big improvement in this strategy would be to adapt the resolution of each of the 10 3D images according to the current size of the growing plant, thus providing a more accurate input at each time step.
The next important step will be to test the method on real plants. This can be achieved through simulations of RoCoCau after calibration based on architectural measurements made in particular species. Predictions could then be made according to the reconstructed target. Yet another way would be to use LiDAR images converted into voxel 3D images as inputs. The second solution would avoid mock-up simulations and probably match reality better. This solution would only apply to the shoot system because no device is currently available to capture the root system without uprooting it. Some technologies including mini-rhizotrons, micro-x-ray chromatography CT, nuclear magnetic resonance (NMR imaging) and positron emission tomography (PET) are already capable of that (Huck and Taylor 1982;Cassidy Steven et al. 2020;Peruzzo et al. 2020) but are mainly used for small plants. Non-destructive observation of root growth of large plants at field scale still remains challenging despite recent attempts using radar (Hruska et al. 1999;Barton and Montagu 2004;Alani and Lantini 2020).
Another interesting use of RoCoCau+TOY with its calibration neural network could be to guess the characteristics of the root system with a known shoot system and known environmental conditions. Knowledge of the extension and depth of a root system is indispensable, especially for tree management in urban environments (Randrup et al. 2001;Bartens et al. 2008;Bartens et al. 2009). Calibration and use of RoCoCau+TOY FSPM based on measurements of a single phenotype.
The purpose of RoCoCau+TOY FSPM is to predict the architectural development of an individual of a given tree species under varying environmental conditions, in other words to predict phenotypic plasticity. To check the quality of RoCoCau+TOY FSPM predictions, measurements of several tree phenotypes are required, meaning measurements have to be taken on several tree stands growing in significantly different environments, as is the case in studies of shoot organogenesis and axis organization (Stecconi et al. 2010;Taugourdeau and Sabatier 2010;Sabatier et al. 2011;Taugourdeau et al. 2011Taugourdeau et al. , 2012Buissart et al. 2018). However, since the TOY model is based on the disentanglement of genome and environment (G×E) interactions, the calibration and use of the RoCoCau+TOY FSPM only require measurements on a single phenotype, hence on a single tree stand.

Calibrating the RoCoCau simulator
Calibrating the RoCoCau simulator for a given tree phenotype requires extensive measurements of architectural traits throughout the architectural development of the tree. As measuring root architectural traits is destructive, measurements have to be made on several trees belonging to the same species and sharing the same environment (at least one tree per architectural development stage). These measurements are then used to calibrate the RoCoCau simulator, allowing it to accurately simulate a botanical representation of the architectural development of a given tree phenotype. Measurement of whole-plant architecture is rare, some approaches with automatic devices are available but only for use on very simple herbaceous plants (Nagel et al. 2012).

Calibrating TOY parameters for exogenous growth on a single
phenotype TOY parameters for exogenous growth are (i) two parameters of the light absorption equation: the light extinction coefficient k, and the maximum useable light energy Light max (w·m −2 ), and (ii) the parameter of the water uptake equation: root absorbing length (RAL) and radius of absorbing soil R ∥ (see Table 1). Light extinction coefficient k and Light max parameters can be obtained by fitting the Beer-Lambert equation with measurements of photosynthetic active radiation (PAR) above, within and under the canopy, completed with the leaf area index (LAI). The root absorbing length (RAL) can be measured directly on the roots of uprooted trees during calibration of RoCoCau. The radius of absorbing soil R ∥ can be obtained from the literature or from measurements (Nadezhdina and Cermak 2003;Corners and Leuschner 2005).
More detailed equations for water uptake and light absorption are available in the literature (Reynolds and Thorney 1982;McMurtrie and Landsberg 1992) and could be used to improve the TOY model. However, it is important to point out that the TOY model relies on the disentanglement of genotype and environmental effects. The unambiguity of the TOY model is a precondition of its ability to predict phenotypic plasticity. When improving or replacing TOY equations for exogenous growth, it is therefore preferable to use only parameters whose calibration is not specific to the observed phenotype, parameters depicting either tree genotype or the environment.

Calibrating TOY parameters for endogenous growth of the tree species on a single phenotype
The 18 TOY hidden parameters for endogenous growth are the focus of the calibration method developed in this study (see section 3). Once the RoCoCau simulators and the exogenous growth parameters are calibrated, it is possible to implement RoCoCau+TOY FSPM with random endogenous growth parameters to build a simulated training data set and to train VoxNet neural networks to predict the values of the endogenous growth parameters.
Since the output formats of RoCoCau and RoCoCau+TOY FSPM are the same, it is possible to provide the trained VoxNet neural networks with a voxelized density map of the RoCoCau simulation of the tree phenotype measured in vivo (see section 5.1.1). The prediction given by the trained VoxNet neural networks for this input are the values of the endogenous growth parameters for the tree species measured in vivo.

Predicting in vivo tree architectural development under varying environmental conditions using RoCoCau+TOY FSPM
Once the RoCoCau simulator and the TOY model have been calibrated, the use of RoCoCau+TOY FSPM relies on the modeller's choice of environmental conditions. Provided with the calibrated values of exogenous growth parameters, the calibrated values of endogenous growth parameters and the chosen values for environment parameters, the RoCoCau+TOY FSPM can predict the architecture of tree measured in vivo if it had grown in the chosen conditions instead of in the actual environmental conditions. The RoCoCau+TOY FSPM can predict the phenotypic plasticity of the tree measured in vivo (Fig. 10).

CON CLUSION
RoCoCau is a new structural whole-plant growth model (describing shoot and root growth and their interactions at the same time). It includes an approach using the architectural unit concept and relies on structural axis categories. Basic functions that explain growth (apical growth, branching, death and self-pruning) are adapted to our specific field of knowledge and experimental measurements made in the shoot and root compartments (Taugourdeau et al. 2012). The RoCoCau simulator can produce 3D mock-ups that are botanically accurate. The functional TOY module is dynamically linked to RoCoCau to express shoot-root compartment dependencies for: -biomass production that requires carbon, water and light; -partitioning to balance shoot and root compartments absorbing contributions and modulating further growth.
The combination RoCoCau+TOY produces an FSPM that can express a variety of growth strategies and reveal plant plasticity under varying environmental conditions. It is possible to fit RoCoCau to a particular species thanks to measurements on both compartments but most of the 25 parameters of TOY are hidden parameters with values that can only be calibrated using numerical tools. It is possible but unlikely that such a large number of values can be discovered using standard methods, which is why we tested the capability of neural networks to accomplish this task. To this end, we adapted the well-known VoxNet convolutional neural network and tested the quality of its predictions using simulated plant data sets for training and validation. The input of the network consists in a 3D voxel representation of split woody and absorbing parts of the plant and the output consists in sets of three values of TOY hidden parameters. The method predicted 6 × 3 TOY hidden parameter values with a very good accuracy after training on a simulated data set and fine-tuning.
The next steps in this project will consist in improvements in two complementary directions. In the biological field, the TOY model will be improved to better match plant knowledge, which needs to be achieved while keeping TOY as simple and compact as possible. The TOY model should include new important features such as nitrogen contribution (Kobe et al. 2010;Liu et al. 2012); better management of the soil or biomass reserves. On the other hand, many actions need to be performed in the numerical and computer science field to improve calibration. The use of convolutional neural networks is promising, but the VoxNet was adapted for TOY calibration in a simple way while changing the last fully connected layer and using a single 3D image containing both a plant growth time sequence and climate values. We are aware that recurrent networks are considered to be more accurate for processing time series. The input content and the loss computing can also be improved. The last step will consist in using real species. The main idea is that, thanks to the botanical accuracy of RoCoCau, it will be possible to replicate plant architecture based on simplified measurements and to train the network on simulated data sets with random values of TOY and then to predict TOY model parameter values that match measured values.

DATA AVA IL A BILIT Y
The source code of the RoCoCau simulator with the plant parameter files used for this study is available on http://amapstudio.cirad.fr/. The source code of the neural network calibration tool can be downloaded from https://github.com/AbelMasson/VoxNet_Strat.