Estimating Sink Parameters of Stochastic Functional-Structural Plant Models Using Organic Series-Continuous and Rhythmic Development

Functional-structural plant models (FSPMs) generally simulate plant development and growth at the level of individual organs (leaves, flowers, internodes, etc.). Parameters that are not directly measurable, such as the sink strength of organs, can be estimated inversely by fitting the weights of organs along an axis (organic series) with the corresponding model output. To accommodate intracanopy variability among individual plants, stochastic FSPMs have been built by introducing the randomness in plant development; this presents a challenge in comparing model output and experimental data in parameter estimation since the plant axis contains individual organs with different amounts and weights. To achieve model calibration, the interaction between plant development and growth is disentangled by first computing the occurrence probabilities of each potential site of phytomer, as defined in the developmental model (potential structure). On this basis, the mean organic series is computed analytically to fit the organ-level target data. This process is applied for plants with continuous and rhythmic development simulated with different development parameter sets. The results are verified by Monte-Carlo simulation. Calibration tests are performed both in silico and on real plants. The analytical organic series are obtained for both continuous and rhythmic cases, and they match well with the results from Monte-Carlo simulation, and vice versa. This fitting process works well for both the simulated and real data sets; thus, the proposed method can solve the source-sink functions of stochastic plant architectures through a simplified approach to plant sampling. This work presents a generic method for estimating the sink parameters of a stochastic FSPM using statistical organ-level data, and it provides a method for sampling stems. The current work breaks a bottleneck in the application of FSPMs to real plants, creating the opportunity for broad applications.

For each DC under the notion of the development cycle, the event during which a meristem produces a phytomer is supposed to take place with a development probability bp, generating a typical Bernoulli process, where p is the physiological age of the corresponding axis. At a given DC, a failed Bernoulli event leads to a void entity in the time series, meaning that no phytomer is generated (de Reffye et al., 1988). The result is an alternation of random series of realized (1) and void (0) entities ( Figure S1). According to the above definition, at a given chronological age n without mortality, the number of phytomers in an axis shows a binomial distribution (n, b).
In plants, different development speeds can be observed between the axes of different branching orders. The ratio between development speeds is called the rhythm ratio w (w<1), taking the highest development speed as 1. In this discrete organogenesis model, development is simulated by a periodic series of real and void entities that describes the functioning of the meristem. For example, w=0.75 corresponds to the series 1110 1110 1110 ( Figure S1). Combined topological axis FIGURE S1 Illustration of the effects of development probability, rhythm ratio and reliability probability in the organogenesis model of an axis.

Reliability probability cp
An axis may end its development for many reasons (e.g., insufficient nutrition or light, insect attack), so at each DC, a reliability probability of a meristem is defined, denoted cp. Once the meristem dies at a given DC, no further phytomers will appear in future DCs, and void entities fill the entire axis from the current phytomer to the tip in the time series ( Figure S1). If b=1 and only the reliability probability c takes effect, the number of phytomers follows a geometrical distribution. The combination of both leads to a compound law, mixing dead axes with still-living axes.
The final simulation result and axis is a scalar product between the above three series. The combination leads to a stochastic axis where the chronological age of each phytomer no longer corresponds with its position on the axis ( Figure S1), which is the case of real plants.

Branching probability ap
Considering branching structure, it is supposed that there are a potential number of axillary branches on each phytomer, denoted m (typically being one or two). We consider that a ramification may occur with a branching probability, denoted ap. Accordingly, the probability that the axillary bud never develops a branch is 1-ap.

Supplement B: the dual-scale automaton (after Kang et al., 2008)
Parallel development of all terminal meristems results in a branching structure. Figure S2 gives an example with the Rauh model, in which each main-stem phytomer can potentially bear one branch. The potential number of phytomers that can bear branches is denoted upq, with p being the PA of the bearing axis and q being the PA of branches. In this example, u12=u13=1. The mean and the variance in the total number of phytomers in the full branching structure can be obtained through generating functions (Kang et al., 2008).

Biomass production
An equation similar to the widely used photosynthesis crop model is used in GreenLab. In the case of a plant with stochastic development, we have the biomass production for a random sample plant s in DC n in (1): where s is the identification of a stochastic plant sample, being an integer value between 1 and sample size; Q s (n) is the biomass increment of plant s at the n-th DC (g plant -1 DC -1 ); LUE is the light use efficiency expressed in gMJ -1 ; PAR is the photosynthetically active radiation above the canopy; k is a cultural coefficient issued from leaf orientation; SP is the production surface, which is a hidden parameter, of the root of (1); LA s (n) is the total plant leaf area of plant s at the n-th DC. Note that LA s (n) is a stochastic variable that differs from one plant to another.

Sink functions
In the GreenLab model, the sink value of an individual organ at the i-th cycle of its expansion is defined as the product of its sink strength and its sink variation function: where po p is the dimensionless relative sink strength of the organ o (leaf, internode and fruit) of physiological age p; the sink strength of a blade of physiological age 1 (main stem) is set to 1 as a reference; to is the duration of the expansion of organ o, assumed to be the same for a given physiological age; Fo is an empirical normalized function describing the evolution of sink strength, modeled here by a Beta law (bell shaped) whose parameters are estimated by model inversion (Guo et al., 2006). Sink function is regarded as the same for both the deterministic or stochastic model.

Plant demand
The plant demand at cycle n, i.e., its overall biomass request, is the sum of all the active sinks in the growing structure.
The plant demand is a convolution on all physiological ages and organ types. No p,s (n-i+1) is the number of organs o at cycle n resulting from the stochastic organogenesis model, with i expansion cycles, born at cycle n-i+1. Another component of plant demand is that of secondary growth (Kang et al., 2012), which is omitted here for clarity.

Biomass allocation
At each DC n, the incoming biomass, ∆qo p,s (i, n), allocated to organ o of age i is proportional to its sink value, Po p (i), and the ratio between the plant production Q s (n-1) and plant demand D s (n): According to this equation, if Q remains the same, fewer created organs will result in smaller plant demand and then larger organs, as can be seen in simulation.
The weight of an organ is then the sum of its increments:

Simulated chronological series
Considering organs of different ages, equation (S5) can be written for age n in a matrix form: where [ , ( * , )] is an n×1 matrix standing for a chronological organic series in an axis of physiological age p; ] is a triangular matrix containing the sink evolution of the organs according to their chronological ages; and [ −1 ] is an n×1 matrix representing the evolution of the supply to demand ratio along the n cycles. The organic series [ , ( * , )] contains the memory of the growth of stochastic sample plant s. From multiple stochastic chronological organic series (indexed with ch), one can obtain the average of all organic series ( * , ) ̅̅̅̅̅̅̅̅̅̅ ℎ . Choosing multiple growth stages, it is possible to determine the evolution of the organ weights along the ranks.

Leaf area
Organ shape dimension can be retrieved from allometric rules. If the specific leaf weight is e, and the leaf weight is , ( , ), then its surface is: This enables the individual plant leaf area, LA s (n), to be computed: where ta is the functioning time of a leaf; , ( − + 1) is the number of functioning leaves of chronological age i; and , ( , ) is the individual leaf weight.

Recursive form of biomass production
Replacing , ( , ) with its value in (S5) finally gives the incoming biomass into the plant at cycle n: The system begins at n=1 from Q0 provided by the seed. The analytical biomass production ( ) shares the same form.

Supplement D: estimating parameters for development -plant crown analysis (after Diao et al., 2012)
The upper part of a plant, or plant crown, is chosen for sampling and provides the data source, as this part is usually the most intact and thus contains the best information about meristem activities. The crown analysis retrieves the development parameters (a, b and c) from the plant architecture measurements (Diao et al., 2012). This approach has been used for numerous species such as coffee (de Reffye et al., 1990), bamboo, cotton, and eucalyptus (Diao et al., 2012). The principles are restated here.

Estimating development parameters b1 and b2 and rhythm ratio w
Take T crowns for analysis with each defined by an axis and the set of first-order branches on it. Note that only living branch axes are chosen for observation. For each axis, the position (phytomer rank) of the branch from the top (noted as K) is recorded, and the number of phytomers in the branch axis (noted as Xs,K, sϵ [1, T] being the identification of the sampled plant crown). The aim is to estimate the development probability for the axis (b1) and the branch axis (b2) as well as the rhythm ratio between them (w) (that of the main axis is set to 1 as a reference).
These probabilities can be estimated from the mean and variance of the number of phytomers in branches located at different positions. In case there are branches of higher order, this analysis can be applied to a pair of orders two and three and so on to obtain the corresponding parameters (Diao et al., 2012).
Suppose there are T crowns to be analyzed, each defined by a stem and the set of branches on it. At rank K from the top of the stem, there is a global mean and a global variance for the number of phytomers of the corresponding branches: The observed global mean is: While the analytical value is: The observed global variance is: The analytical global variance is: The local variance between two branches on the same stem separated by L phytomers is: The analytical local variance is: ( ) = • 2 (1 − 2 ) + 2 1 + ( • 2 ) 2 (1 − 1 ) 2 1 2 + ( • 1 ) 2 2 2 2 2 (S14) Solving the system, we obtain b1, b2 and w: The computation of b1, b2 and w can be done for each selected rank k, and when more k are selected, their average values are more accurate. For two branches belonging to the same whorl, we have L=0, while we have L=1 for two adjacent branches.

Estimating the branching probability a
The branching rate at rank K, a(K), is estimated as the ratio between the number of observed branches and the number of available axillary meristems.

Estimating the reliability probability c
The abortion rate Fc1(K) at the phytomer of rank K is the ratio between number of the dead branches and the total number of branches observed. This proportion depends on the probability of bud survival ci at each DC i from 1 to n. Data assimilation is performed using the mean of target files to be fitted by the model that describe the status of the phytomers on the stem at rank K from the axis tip, i.e., how many phytomers get branches and whether they are alive or dead. Due to the Bernoulli process, the average age at rank K is: Rewriting the rank-governed reliability function Fc1(K) on the thermal time basis Fc2(n) using equation, the reliability, cn, of the meristem, i.e., the probability to survive at cycle n, can be deduced from the following recurrence system: If parameter c is constant, we have: Generally, Fc2 is a sigmoid that can be assessed by two parameters, Ba and Bb, in a process of n duration from meristem birth to death.

Supplement E: Example of the experiments of simulation and computation using the software 'Gloups'
The experiments of simulation and computation was performed in the software 'Gloups' developed under Windows environment using MATLAB. The main procedures include one for simulation and another for calibration.

Simulation of stochastic development
Monte-Carlo method is applied for the simulation of stochastic development. In Gloups, 'ZSub_GL5_2S' is used for the simulation. The growth cycle of plants can be set as needed in 'ZSub_GL5_2S'.
To obtain the simulated chronological organic series for a single stem, the parameter values in Line 4 needs to be set to 1 for chronological and 0 for topological mode, and that in Line178 needs to be set to 0.7 both for the main stem and branches, as follows:

Calibration of stochastic development
In Gloups, 'ZSub_GL5_2S' is used for the parameter estimation. In calibration, it calls the analytical computation of organic series according to initial parameters values, the result is to be fit with the target data, either from simulation, or from real plants.

Calibration from the virtual target data
For Figure 6, as we have obtained the virtual target data 'targ.m' renamed as 'Fig6_tar.m', we can estimate the parameters using 'ZSub_GL5_2A'. Run 'ZSub_GL5_2A', choose the parameter file 'Fig 6_par.m' and 'Fig 6_tar.m', the parameter values can be obtained, and the estimated weights of each organ according to their positions, the source-sink ratio, etc., will be given. Generally speaking, the source-sink parameters can be estimated simultaneously or separately. In the parameter files, we can choose which parameters we want to compute.

Calibration from the simulated target data
For Figure 9, as we have obtained the measured data (weights of blade and internode, number of phytomers on the main stem and branches), first, we compute the probabilities of a, b and w with the crown analysis; Second, the values of these probabilities can be set in the parameter file 'Fig9_par.m'; Third, run 'ZSub_GL5_2S', choose parameter file 'GLOUP\Target\Fig 9_par.m', thus, the virtual target data file can be obtained in 'targ.m', then replaced the target data by the measured data, named as 'Fig9_tar.m'; Finally, we can estimate the parameters with the method mentioned above.