Antonin Della

Della Noce

Véronique Letort

Amandine Hansart

Charlotte Baey

Gautier Viaud

Sébastien Barot

Jean-Christophe Lata

Xavier Raynaud

Paul-Henry Cournède

Jacques Gignoux

Antonin Della Noce

Modeling the inter-individual variability of single-stemmed plant development

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The concept of thermal time and its use as time scale to model plant development has a long history since its first introduction by [START_REF] Allee | Principles of animal ecology[END_REF], [START_REF] Trudgill | Thermal time-concepts and utility[END_REF]. It is usually the first and mandatory step for building a functional-structural plant model (FSPM): relying on the assumption that organogenesis is mainly monitored by thermal time [START_REF] Gramig | Leaf Appearance Base Temperature and Phyllochron for Common Grass and Broadleaf Weed Species[END_REF], it consists, for single-stemmed plants (e.g. crops or trees in their early stages i.e. before the appearance of secondary axes), in expressing the number of phytomers along its main axis as a predefined function of thermal time.

Thermal time (expressed in degree-day) can be defined as:

τ (t) = t t0 max(0, T (s) -T b )ds ( 1
)
where t is the time expressed in days, T (t) is the temperature in • C, t 0 is the emergence time of the plant, i.e. the date at which its first phytomer is observed, and T b is the base temperature, below which the organogenesis process stops.

If the function of thermal time is chosen deterministic, plant organogenesis should be the same from an individual to another, as long as they share the same environmental conditions. However, in many cases (see Figure 1 for an illustration), a strong inter-individual variability is encountered, that cannot be straightforwardly associated to any environmental or microenvironmental factors. It has nevertheless to be accounted for in the models and cannot be neglected since it is an important characteristic that can influence the model behavior and interpretations. Besides, by integrating random effects into the model, parametric estimations at population level and at individual level are intrinsically coupled : the predicted development of an individual is consistent with the overall trend of the phenomenon in the population, with individual characteristics matching the observed population variability.

Heterogeneity was studied in a sugar beet population in [START_REF] Baey | Modelling the interindividual variability of organogenesis in sugar beet populations using a hierarchical segmented model[END_REF] for the organogenesis and in [START_REF] Baey | A non linear mixed effects model of plant growth and estimation via stochastic variants of the em algorithm[END_REF] for a full FSPM model. The approach is based on the identification of non-linear mixed effect models with the Stochastic Approximation Expectation Maximization (SAEM) algorithm as proposed by [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF]. In such framework, the parameters of the model describing the individual behavior, such as phyllochrons, are assumed to be randomly distributed over the population. In this article, we follow such a population approach using hierarchical models for plant development, developing in details each step of the identification procedure.

The method will be presented using a specific case-study that will be used all along the paper to illustrate the method and interpretation of the results. This dataset was collected on young Acacia erioloba plants, grown in CEREEP greenhouse. Acacia Erioloba is a key species of savanna ecosystems. In some areas, it is a dominant component of the tree cover, providing food (specifically nutritive pods during the dry season [START_REF] Barnes | Seed predation, germination and seedling establishment of Acacia erioloba in northern Botswana[END_REF]) and shadow to many browser herbivores and possibly enriching the soil through nitrogen fixation. A noticeable plasticity between individuals was observed in [START_REF] Seymour | Grass, rainfall and herbivores as determinants of Acacia erioloba (Meyer) recruitment in an African savanna[END_REF] and explained by the plant adaptability to the multiple constraints it is subject to: fires, drought, herbivory. The design of mechanistic models to predict the dynamics of savana tree cover according to environmental conditions is of prior importance for the man-agement of the ecosystem resources. A first attempt to capture the morphological variability induced by browser aggressions was undertaken for a closely related species, Acacia tortillis, growing in similar biotopes, using a GreenLab FSPM in [START_REF] Mathieu | Simulation of morphological plasticity of Acacia tortilis in response to herbivore attacks[END_REF]. It was able to reproduce some observed architectural patterns such as the ball-shaped crown resulting from intense grazing but was lacking proper follow-up of the individual growth trajectories to be validated, especially for the organogenesis part. This is why new experiments were conducted, with the objective of a better understanding of the growth of Acacia in controlled conditions. The study of organogenesis appears a first step for the design of a complete GreenLab FSPM, as architecture monitors biomass allocation.

Section 2 sets the notations of the problem and develops in details a step of the iterative Expectation-Maximization procedure ( [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF]). According to the model formulation, different versions of the EM algorithm are derived. In section 3, the method is applied to the identification of different models of organogenesis.

II. DATA AND HIERARCHICAL MODELS

A. Experimental conditions and observed variability

A large measurement campaign has been carried out since 2013 on Acacia erioloba specimens in the greenhouse of CEREEP-Ecotron IleDeFrance (http://www.foljuif.ens.fr/) to analyze the dynamics of development of trees according to their environmental conditions. Plants were growing in regularly spaced pots, so that the effect of competition for light resources can be neglected. The study focuses on single axis plants, which germinated in July 2013. Weekly measurements of phytomers number, along with daily recordings of the temperature inside the greenhouse, were gathered in a dataset, totalling 567 observation points (temperatures, numbers of phytomers), corresponding to the first five months of development of 122 plants. Figure 1 gives a visualization of the number of phytomers according to the cumulated temperature, i.e. if germination date of a plant is denoted by t 0 , then the associated cumulated temperature at day t is t s=t0 T s .

A global pattern can be observed on Figure 1 : the growth has essentially two stages, characterized by two distinct phyllochrons. Around this overall behavior, important variations of growth rate can be noticed from one individual to another. A piecewise linear growth can be a good approximation of the dynamics, similarly to the model adopted in [START_REF] Baey | Modelling the interindividual variability of organogenesis in sugar beet populations using a hierarchical segmented model[END_REF]. The objective is to design a model that reproduces this overall dynamics and takes into account this inter-individual variability.

B. Design of two mixed-effect models

For simplicity sake, a non generic definition of mixed-effect model is given in this section. A general formulation can be found in [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF].

Let y = (y i ) 1≤i≤N be the experimental observations for all individuals of the population, i.e. for all i ∈ 1; N , y i is a vector in R ni , and for all j ∈ 1; n i , y ij is the number of phytomers of plant i at the observation time t ij . t i0 denotes the observed emergence date of plant i. Let T b be the base temperature specific to Acacia erioloba, which is unknown a priori, its estimation being an interesting aspect of this study. Then the thermal time associated to the plant i at time t ij > t i0 is given by a discretization of equation [START_REF] Allee | Principles of animal ecology[END_REF].

τ (t ij , T b ) = τ ij (T b ) = tij -1 s=ti0 max(0, T (s) -T b ) (2)
with the initialization τ (t i0 ) = 0. In the previous equation, T (s) is the averaged temperature at day s ∈ t i0 ; t ij . A model of organogenesis can therefore be formulated at the individual level as follows, for all i ∈ 1; N and for all j ∈ 1; n i

y ij = 1 + f (τ ij , φ, α) + ij ( 3
)
where f is a function of the thermal time, parameterized by φ ∈ R d and α ∈ R m , and ij ∼ i.i.d.

N (0, σ 2 ) are additive centered Gaussian noise processes. As the germination date t i0 is the reference date for the thermal time, and as the number of phytomers in the plant at that date is 1, we must have f (0, φ, α) = 0.

To model the inter-individual variability, it can be assumed that for each plant i a parameter vector φ i is sampled from a multivariate normal distribution in the population

φ i ∼ i.i.d. N (φ p , Σ p ) (4)
with (φ i ) 1≤i≤N , the set of individual parameters, or random effects, ξ = (φ p , Σ p ), the population parameters. α are the fixed effects, that is to say that these parameters are constant within the population. Conditionally to φ i , the distribution of the observation model is given by equation ( 3)

y ij |φ i ∼ N (1 + f (τ ij , φ i , α), σ 2 ) (5)
In the following, in order to alleviate the derivations, the notation ỹij = y ij -1 will be used. For convenience, the dependence of the model with respect to the thermal time may be occasionally written implicitly by the notation

f i (φ i , α) = (f (τ ij , φ i , α)) 1≤j≤ni ∈ R ni .
We have chosen two candidate models f to reproduce the organogenesis dynamics of Acacia. Both of them represent the organogenesis as a piecewise linear function, following the guidelines given in [START_REF] Baey | Modelling the interindividual variability of organogenesis in sugar beet populations using a hierarchical segmented model[END_REF], but the first model is linear with respect to the individual parameters, whereas the other one is non-linear. a) Model 1 : Linear with respect to the individual parameters: For simplicity sake, let us first consider the case of the linear model w.r.t. φ. The model is formulated as follows,

for all i ∈ 1; N f (τ ij , φ i , α) = a i (τ ij -τ 0 )I{τ ij > τ 0 }+b i (τ ij -τ c )I{τ ij > τ c } (6) Here, the individual parameter is φ i = (a i , b i ) T ∈ R 2 are
the phyllochrons associated with the two stages of organogenesis: the first one delimited between initiation thermal time τ 0 and the thermal time of the change point τ c and has a rate a i ; the second one takes place after τ c and has a rate a i + b i .

The full parameter vector for our model is denoted θ and θ = (φ p , Σ p , T b , τ 0 , τ c , σ 2 ). The constant parameters that are connected to the notion of thermal time are called the thermal parameters and are denoted α = (T b , τ 0 , τ c ). To underline the linear aspect with respect to the individual parameters, we will use the notation f i (φ i , α) = f i (α)φ i . A descriptive model which is linear with respect to the random effects and non-linear with respect to the fixed effects was studied in the case of orange tree growth in [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS[END_REF]. b) Model 2 : non-linear mixed effects: The assumptions of the previous model are relaxed concerning the thermal times. The equation is the same as in [START_REF] Kuhn | Maximum likelihood estimation in nonlinear mixed effects models[END_REF], except that here τ 0 and τ c are also considered variable in the population and

φ i = (a i , b i , τ i 0 , τ i c ) T ∈ R 4 .
The full parameter vector is now θ = (φ p , Σ p , T b , σ 2 ) with φ p ∈ R 4 , and the only fixed effect parameter being the base temperature α = T b .

In a frequentist approach, we want to estimate the parameters θ by maximizing the likelihood given the observations. The likelihood is given by:

p θ (ỹ) = N i=1 R d N (φ i ; φ p , Σ p )   ni j=1 N (ỹ ij ; f (τ ij , φ i , α), σ 2 )   dφ i
where we denote N (x; µ, Σ) the probability density function of N (µ, Σ) evaluated at x. The evaluation and/or the maximization of such likelihood, in the case of both model 1 and model 2, is analytically intractable. Expectation-Maximization (EM) algorithms, designed for the estimation in incomplete data problems, is the intrinsic method to identify hierarchical models, by generating a sequence (θ n ) such that (p θ n (y)) is increasing. The next subsection specifies how such algorithm is applied to the identification of our two models.

III. AN ITERATION OF THE EM ALGORITHM

We recall the principles of the EM algorithm from [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF], and how it can be adapted to different model formulations. Let us suppose that we have an estimate θ n of the parameter θ from the initialization or from a previous iteration. The observation log-likelihood associated to individual i i (θ) = log p θ (ỹ i ) can then be decomposed as follows

i (θ) = R d i (θ)p θ n (φ i |ỹ i )dφ i (7) = R d log p θ (ỹ i , φ i )p θ n (φ i |ỹ i )dφ i (8) - R d log p θ (φ i |ỹ i )p θ n (φ i |ỹ i )dφ i (9) = Q i (θ|θ n ) + H i (θ|θ n ) (10)
where the derivation of the lines ( 8) and ( 9) is obtained by applying Bayes' rule. The first term ( 8) is the expected joint log-likelihood

Q i (θ|θ n ) = E θ n [log p θ (ỹ i , φ i )|ỹ i ]. The second term (9) H i (θ|θ n ) = -E θ n [log p θ (φ i |ỹ i )|y i ]
is the cross entropy of the two posterior distributions. Using Jensen inequality, it can be proved that for all θ, H i (θ|θ n ) ≥ H i (θ n |θ n ). To maximize i (θ), we are therefore interested in maximizing the total expected joint likelihood

θ n+1 = argmax θ N i=1 Q i (θ|θ n ) = argmax θ Q(θ|θ n ) (11)
A qualitative description of EM can be formulated as follows : at each iteration of the algorithm, we use our current estimate θ n to generate the latent variables φ i , then we maximize the joint log-likelihood Q(θ|θ n ).

The next paragraphs specify the expectation and the maximization steps for models 1 and 2.

A. Expectation step

The objective of this step is to find a close form of Q(θ|θ n ). For all i ∈ 1; N , let φ i be a realization of an individual parameter, then for all θ

log p θ (y i , φ i ) = log p θ (ỹ i |φ i ) + log p θ (φ i ) log p θ (ỹ i , φ i ) = - n i + d 2 log(2π) -n i log(σ) - ỹi -f i (φ i , α) 2 σ 2 -log det(Σ p ) - 1 2 (φ i -φ p ) T Σ -1 p (φ i -φ p )
The expectation step therefore requires the evaluation of the following functions

MSE n i (α) = E θ n ỹi -f i (φ i , α) 2 |ỹ i (12)
D n i (φ p , Σ p ) = E θ n (φ i -φ p ) T Σ -1 p (φ i -φ p )|ỹ i (13)
a) Model 1: For model 1, thanks to the linearity with respect to φ i , the posterior distribution p θ n (φ i |y i ) can be analytically derived, giving close form expressions for the previous functions with respect to θ. Let us derive those expressions. p θ n (φ i |y i ) is proportional to the following distribution

p θ n (φ i |y i ) ∝ p θ n (y i |φ i )p θ n (φ i ) (14)
∝ exp - ỹi -f i (α n )φ i 2 2(σ n ) 2 (15)
- 1 2 (φ i -φ n p ) T (Σ n p ) -1 (φ i -φ n p ) (16)
The previous equation shows that φ i |y i is normally distributed, and its covariance matrix and mean are respectively given by

Σ n i = (Σ n p ) -1 + f i (α n ) T f i (α n ) (σ n ) 2 -1 (17)
µ n i = Σ n i (Σ n p ) -1 φ n p + f i (α n ) T ỹi (σ n ) 2 (18)
The expectation ( 13) is therefore

D n i (φ p , Σ p ) = E θ n (φ i -µ n i ) T Σ -1 p (φ i -µ n i )|ỹ i + (µ n i -φ p ) T Σ -1 p (µ n i -φ p ) D n i (φ p , Σ p ) = Tr Σ -1 p Σ n i + (µ n i -φ p ) T Σ -1 p (µ n i -φ p ) (19)
Similarly, for expectation (12)

MSE n i (α) = Tr f i (α) T f i (α)Σ n i + ỹi -f i (α)µ n i 2 (20)
The expected joint log-likelihood is obtained by combining all the individual terms

Q(θ|θ n ) = - M + N d 2 log(2π) -M log(σ) - 1 2σ 2 N i=1 MSE n i (α) (21)
- N 2 log det(Σ p ) - 1 2 N i=1 D n i (φ p , Σ p )
where the total number of observations is M = N i=1 n i . b) Model 2: : The situation in the non-linear case is more complex since the first and second order moments of p θ n (φ i |y i ) cannot be derived analytically as in the linear case. The idea is then to estimate the quantities of interest MSE and D n i by generating samples from the distribution. It is only known up to a multiplicative constant

p θ n (φ|ỹ i ) ∝ exp - ỹi -f i (φ i , T n b ) 2 2(σ n ) 2 - 1 2 (φ i -φ n p ) T (Σ n p ) -1 (φ i -φ n p ) (22) p θ n (φ|ỹ i ) ∝ exp(-J n i (φ i )) (23)
MCMC techniques, such as Metropolis-Hasting, enables the generation of a Markov Chain, whose stationary distribution is the target distribution (22), and the moments of interest can be computed by the ergodic theorem. A MCMC simulation within the Expectation step of p θ n (φ i |ỹ i ) is therefore adopted here (for a a more detailed description on MCMC-EM methods, see [START_REF] Fort | Convergence of the monte carlo expectation maximization for curved exponential families[END_REF]). An adaptative Metropolis-Hasting algorithm, introduced by [START_REF] Haario | Adaptive proposal distribution for random walk Metropolis algorithm[END_REF], is used in this paper. The procedure is described in the following frame Algorithm 1 Adaptative Metropolis-Hasting algorithm [START_REF] Chen | Inférence bayésienne dans les modèles de croissance de plantes pour la prévision et la caractérisation des incertitudes[END_REF] Let L be the length of the chain, λ 0 > 0, β * ∈]0, 1[ and (γ k ) k∈N a sequence of learning rates.

for i ∈ 1; N do Initialization find good values for the starting point of the chain φ 0 i , the initial covariance of the proposal distribution Σ 0 , and the initial mean µ 0 for covariance update Iterations

for k ∈ 1; L do 1. sample a new candidate z ∼ N (φ k-1 , λ k-1 Σ k-1 )
2. select z with probability

β(φ k-1 i , z) = min p θ n (z|ỹ i ) p θ n (φ k-1 i |ỹ i ) , 1
3. select z with probability β :

φ k i = z or reject it φ k i = φ k-1 i . 4. update the mean µ k = µ k-1 + γ k (φ k i -µ k-1 ) 5. update the covariance Σ k = Σ k-1 + γ k (φ k i -µ k )(φ k i -µ k ) T -Σ k-1 6. Update the scale of the covariance λ k = exp log(λ k-1 ) + γ k (β(φ k-1 i , z) -β * ) end for Add the chain (φ k i ) to the stack end for return the stack of all chains (φ k i ) 1≤k≤L 1≤i≤N
Efficiency of the MCMC algorithm is sensitive to its initialization. We propose to choose φ 0 i maximizing p θ n (φ|ỹ i ), which is equivalent to minimizing the quantity:

J n i (φ i ) = 1 (σ n ) 2 ỹi -f i (φ i , T n b ) 2 +(φ i -φ n p ) T (Σ n p ) -1 (φ i -φ n p ) (24)
As the function f i is non-smooth, we have chosen to use the Nelder-Mead algorithm, embedded in the function optimize in the Optim package (JULIA language), to compute φ 0 i . At this point, the posterior density has high values, and the acceptance rate α is likely to be non zero (unless this maximum is too sharp). Therefore, the generated chain will potentially not get stuck in an area where the target density is low. The covariance of the proposal is initialized at the population covariance, so that the update of the parameter is adapted to the scale of each component φ i .

We use the generated samples (φ k i ) 1≤k≤N for the Monte-Carlo estimation of MSE n i and D n i for each individual i

MSE n i (T b ) ≈ 1 L -K + 1 L k=K ỹi -f i (φ k i , T b ) 2 D n i (φ p , Σ p ) ≈ Tr(Σ -1 p Σn i ) + (μ n i -φ p ) T Σ -1 p (μ n i -φ p )
where K is a burn-in period,

μn i = 1 L -K + 1 L k=K φ k i is the empirical mean and Σn i = 1 L -K + 1 L k=K (φ k i -μn i )(φ k i - μn i )
T is the empirical covariance. Injecting those expressions in 21, we obtain the joint-likelihood of the expectation step. c) Alternative E-step for Model 2: An alternative to MCMC simulation of the posterior distribution appears when considering its shape. As it is shown in the Results section, p θ n (φ|ỹ i ) are unimodal distributions, very sharp around their respective modes. Such situation motivated [START_REF] Lindstrom | Nonlinear Mixed Effects Models for Repeated Measures Data[END_REF] and [START_REF] Vonesh | A note on the use of Laplace's approximation for nonlinear mixed-effects models[END_REF] to approximate the target by a multivariate Gaussian, either by a linearization in the neighborhood of the fixed effects, either by using the Laplace approximation.

The qualitative principle of such approximation is as follows: near the argmax of p θ n (φ i |y i ) (that is to say near the argmin of J n i ), denoted by µ n i , the distribution can be approximated by the second order development

J n i (φ i ) = J n i (µ n i ) + 1 2 (φ i -µ n i ) T ∂ 2 φi J n i (µ n i )(φ i -µ n i ) (25)
as the gradient is zero, leading to

p θ n (φ i |ỹ i ) ∝ exp 1 2 (φ i -µ n i ) T ∂ 2 φi J n i (µ n i )(φ i -µ n i ) (26) If φi is an isolated maximum, the hessian ∂ 2 φi J n i ( φi ) is positive definite, and p θ n (φ i |ỹ i ) can be approximated by a Gaussian distribution N µ n i , -∂ 2 φi J n i (µ n i ) -1
.

A problem specific to model 2 is the fact that it is nonsmooth (because of the indicator function) and taking second order derivatives of f i (φ i , T b ) might be numerically unstable. A C ∞ approximation of the model 2 was considered, using a sharp sigmoid function σ as a surrogate of the indicator function I{. > 0}1 .

σ (x) = 1 1 + exp -x L 2 (R) -→ →0 I{x > 0} (27)
Similar differentiable approximations of piecewise-model were used in [START_REF] Tishler | A Maximum Likelihood Method for Piecewise Regression Models with a Continuous Dependent Variable[END_REF]. This enables to have analytic expressions of the derivative in every points of the parameter space.

The minimum µ n i of J n i can be obtained numerically, once again by using Nelder-Mead or Gauss-Newton algorithms for instance. The expression of the covariance matrix Σ n i as a function of the model derivatives is

Σ n i =   (Σ n p ) -1 + 1 (σ n ) 2 ni j=1 ∂ φi f (τ ij , µ n i )∂ φi f (τ ij , µ n i ) T (28)
+ 1 (σ n ) 2 ni j=1 (f (τ ij , µ n i ) -ỹij )∂ 2 φi f (τ ij , µ n i )   -1 (29)
The previous equation is consistently equal to (17) if f is linear with respect to the individual parameters. The positivity of the matrix holds since µ n i is a minimizer of J n i . The deviation D n i (φ p , Σ p ) has with these notations exactly the same expression as in [START_REF] Mathieu | A dynamic model of plant growth with interactions between development and functional mechanisms to study plant structural plasticity related to trophic competition[END_REF]. However, because of the nonlinearity, the mean square error MSE n i has to be estimated by an empirical mean, using an i.i.d sample

(φ k i ) 1≤k≤L ∼ N (µ n i , Σ n i ) MSE n i (T b ) ≈ 1 L L k=1 ỹi -f i (φ k i , T b ) 2 (30)

B. Maximization step

This section derives the method to find the maximum argument of the function Q(θ|θ n ) for model 1 and model 2.

a) Model 1: In the case of the linear model, the update of the parameters φ p , Σ p , σ is analytic, while the update of the fixed effects α (base temperature T b and thermal times τ 0 , τ c ) are performed by the numerical minimization of the

functional MSE n = N i=1 MSE n i (α).
Let us start by the updates of φ p and Σ p which are obtained by solving the respective equations

∂ φp Q = 0 and ∂ Σp Q = 0. φ n+1 p = 1 N N i=1 µ n i (31) Σ n+1 p = 1 N N i=1 Σ n i + (µ n i -φ n+1 p )(µ n i -φ n+1 p ) T (32)
As for the fixed effects α, it is interesting to notice that the quantity to minimize can be written as

MSE n (α) = N i=1 ỹi -f i (α)µ n i 2 + N i=1 Tr(f i (α) T f i (α)Σ n i )
(33) which appears as a bias and variance decomposition. The update of the fixed -effect estimate is obtained by minimizing this mean square error term.

α n+1 = argmin α MSE n (α) (34)
This minimization is performed using a simplex algorithm as for equation 24. This procedure estimates at the same time the optimal change points τ 0 and τ c , and also the base temperature T b . As a matter of fact, the procedure does not take into account the fact the model is piecewise-linear. A methodology specific to the identification of change points in piecewise-linear models were derived in [START_REF] Hawkins | Point Estimation of the Parameters of Piecewise Regression Models[END_REF]. Unfortunately, such procedure cannot be straightforwardly applied here as it requires to estimate the slopes and the change points at the same time, whereas in our case the slopes are computed during the expectation stage, prior to the estimation of the change points. Finally the update of the observation noise parameter σ is given by the equation ∂ σ Q = 0, which leads to

σ n+1 = MSE n (α n+1 ) M (35)
b) Model 2: The maximization step in the case of the non-linear model follows exactly the same pattern as previously, except that the characteristics of the posterior distribution µ n i and Σ n i have to be replaced by their Monte-Carlo estimates μn i and Σn i or by the characteristics of the Gaussian posterior derived by Laplace approximation. Besides, the MSE function is only dependent on the base temperature T b , such function can be efficiently minimized using univariate optimization methods like the Brent method.

IV. RESULTS

Parameter identification of hierarchical models enables multi-level predictions : at the population level, the distribution of the observations depends on the estimated parameters φ p , Σ p and on σ 2 ; at the individual level, a distribution of observations is also identified through the individual estimates.

A. Identification of Model 1

In the case of the linear model, the EM algorithm was run for 200 iterations, starting from initialization values that are presented in Table IV-A. After 150 iterations, the sequence of (Q(θ n+1 |θ n )) n was quite stable, so as the standard deviations of the observation (σ n ) n . The increase of Q(θ n+1 |θ n ) ensures that the parameters update is going in a good direction, just like the decrease of the standard deviation. The following table gives the estimates of the model parameters after convergence. At the population level (Figure 2), the observations made at thermal time y(τ ) = 1 + f (τ ) T φ + are normally distributed:

y(τ ) ∼ N 1 + f (τ ) T φ p , σ 2 + Tr f (τ )f (τ ) T Σ p (36)
which gives the confidence intervals of the prediction distribution. We can also check in Figure 3 whether the shape of the empirical distribution of the individual means (µ n i ) 1≤i≤N is close to a Gaussian. If the number of individuals in the population is too small, this graphical representation might be inappropriate. In conclusion, the linear model shows a good adequacy with the observations. As model 2 is a more generic formulation of model 1, in the sense that one can go from one to another by setting the variance of τ 0 and τ c to zero, the parameter values estimated in the linear case can be used to initialize the EM algorithm performed in the non-linear case.

B. Identification of Model 2

For the non-linear model, a MCMC-EM was first performed for 100 iterations starting from the values estimated in the previous model. The length of the Monte-Carlo Markov chain was set to 3000, with a burn-in period of 1000. The parameter λ of Metropolis-Hasting algorithm was first initialized to 1 and the optimal acceptance rate was set to β * = 0.234, following advice in [START_REF] Haario | Adaptive proposal distribution for random walk Metropolis algorithm[END_REF]. Such model did not improve significantly the mean squared error compared to the first one, in terms of mean square error and observation noise: the initial value of σ was 0.697 (see previous model) and it ends up being equal to 0.67. The base temperature was also very stable, the final estimate being 14.4. Nevertheless it was observed, by considering the sample generated by Metropolis-Hasting algorithm, that the Gaussian approximation could be a quite good estimate of the posterior distribution in itself, leading us to implement a Laplace approximation method.

After convergence of the EM with the Laplace approximation, we have noticed that the graph of the predicted population distributions and individual distributions were quite similar to the ones obtained in the linear case (Model 1). The standard deviations of parameters τ 0 and τ c were found to be critically low, respectively 9×10 -3 • C.d and 1×10 -2 • C.d. At this point, there are two possible explanations of such results :

• either modeling initiation thermal time and thermal time of the change point as random effects is inappropriate and the linear assumption is justified • or the fact to initialize the EM algorithm in the region of parameter space corresponding to the configuration where those two parameters were constant at the population level maintains them in this disposition. In that case, a new initialization procedure needs to be investigated.

V. DISCUSSION AND PERSPECTIVES

The organogenesis of young Acacia erioloba provides a new and interesting case of application of hierarchical models. Such population approach seems relevant: the heterogeneity within the population is quite well represented by simple linear models, whose identification with EM algorithm is analytic. The non-linear version brought only minor improvements, due to the fact that the thermal times of change point and initiation could be appropriately modeled by fixed effects in our case.

The methodology developed could be easily extended to other non-linear mixed effects models. Indeed, the first method only requires linearity with respect to the individual parameters, and such assumption can be justified in a wide range of dynamics. In our case, the model is a linear combination of piecewise linear functions, but other families of non-linear functions, such as sinusoids, Gaussian kernels would have lead to the same derivations. Identifying a linearized version of the model provides relevant initialization values to EM algorithms applied on a generic form of the model, with non-linearity with respect to the individual parameters.

Although the presentation and illustration of our methodology was the core of this paper, this study is also a step towards a more complete modeling of Acacia erioloba development using a complete GreenLab FSPM. For further stages of development, with branchings and reiterations, modeling the rhythm of appearance of phytomers as a function of the thermal time is no longer sufficient as current architecture plays a role. The development of the branching patterns involves complex processes that can be modeled through the retroaction of plant functioning on architecture (see [START_REF] Mathieu | A dynamic model of plant growth with interactions between development and functional mechanisms to study plant structural plasticity related to trophic competition[END_REF] within the GreenLab framework) or with stochasticity ( [START_REF] Loi | Generating functions of stochastic Lsystems and application to models of plant development[END_REF]).

Finally, the real objective is to derive a full population model for Acacia erioloba based on an individual FSPM in order to account for the inter-individual variability. Likewise, and as introduced by [START_REF] Baey | A non linear mixed effects model of plant growth and estimation via stochastic variants of the em algorithm[END_REF], both organogenesis and functional parameters should be modelled as mixed effects. The main difficulty will rely on the joint estimation of organogenesis and functioning, especially in the case of interaction between both.

Figure 1 .

 1 Figure 1. Number of phytomers along plant main axes with respect to the cumulated sum of temperature, i.e

Figure 2 .

 2 Figure 2. Predicted distribution of the observations at the population level

Figure 3 .

 3 Figure 3. Empirical and theoretical distributions of the individual means

Figure 4 .

 4 Figure 4. Prediction of the observations distribution at the individual level

This C ∞ approximation of the indicator function is obtained by taking its convolution with the approximate Dirac distribution 1 σ (x)(1 -σ (x)) → →0 δ 0 in the sense of the distributions.


P Reinhardt
email: peter.reinhardt@upmc.fr

A possible valence-bond approach to symmetry-adapted perturbation theory

The multi-configurational ansatz of valence-bond theory may serve as basis for calculating intermolecular interaction energies in a non-orthogonal basis. We look in the present contribution at the possibility to obtain the 1st-order electrostatic interactions from breathing-orbital valencebond densities, and the 2nd-order dispersion energy from dipole-dipole interactions. The discussion is based on numerical results for the interaction of two N 2 O molecules.

Introduction

Since the beginnings of describing electronic structure in terms of orbitals, the two pictures of completely localized bonds or completely delocalized molecular orbitals compete [START_REF] Heitler | Wechselwirkung neutraler Atome une homöopolare Bindung nach der Quantenmechanik[END_REF][START_REF] Mulliken | The Assignment of Quantum Numbers for Electrons in Molecules. I[END_REF][START_REF] Brush | Dynamics of Theory Change in Chemistry: part I -the benzene problem 1865 -1945[END_REF]. Both pictures have their advantages and shortcomings, and obvious applications. From molecular orbitals we have Hartree-Fock theory and beyond all the methods for describing dynamical electron correlation. In cases of bond breaking, multi-reference methods have been developed, adding non-dynamical electron correlation in situations where relying on a single reference is not any more justified. This machinery becomes very demanding on computer resources for handling immense amounts of data due to the large number of Slater determinants or Configuration State Functions (CSF) involved, and their interactions.

In a chemist's language this apparently brute-force approach is difficult to treat -bonds, bond formation, bond breaking and electron rearrangement need only few electrons involved in a picture, despite the indiscernibility of elementary particles and long-range Coulomb interactions between them.

Valence Bond theory tries to describe these essentials of a chemist's view on molecules and matter. It relies on one or several Lewis structures (i.e. electron pairs assembled in bonds or lone pairs, or left alone in singly occupied orbitals). Each of these structures may be developed in a common (localized) orbital set, or in a more modern formulation [START_REF] Hiberty | Compact valence bond functions with breathing orbitals: Application to the bond dissociation energies of F 2 and FH[END_REF], in orbital sets for each involved structure, the Breathing Orbitals Valence Bond (BOVB) approach. The main obstacle for routine calculations in this type of development is the inherent non-orthogonality of the orbital sets. All expressions have to carry overlap matrix elements throughout, rendering calculations much heavier to be executed than expansions in determinants in common orthogonal orbitals, where today one may handle dimensions of millions or even several billions of determinants.

Ground-state energies, geometries and thermochemistry may be obtained nowadays with routine calculations at experimental precision, using gradients and higher derivatives, and sophisticated optimization algorithms like for instance meta-dynamics [START_REF] Laio | Escaping free-energy minima[END_REF]. A continuation of this is to include a chemical environment, via explicit solvatation or by implicit models of a solvatation sphere [START_REF] Mo | Block-localized wavefunction (BLW) based two-state approach for charge transfers between phenyl rings[END_REF]. The first step toward this is the study of intermolecular interactions for dimer systems. Ideally, the interaction may be constructed from the separate monomer wavefunctions, which gave rise to a series of developments known as symmetry-adapted perturbation theory (SAPT) [START_REF] Jeziorski | Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes[END_REF]. The total energy of the dimer system is constructed from a triple perturbation series, starting from the separate Hartree-Fock orbitals (wavefunctions) of the respective monomers, including then intermolecular interaction as one perturbation, and the dynamical, intra-molecular correlation of the monomers as a second and third one. Recently, the latter perturbation series could be replaced with good precision via density-functional theory [START_REF] Hesselmann | Density-functional theory-symmetry-adapted intermolecular perturbation theory with density fitting: A new efficient method to study intermolecular interaction energies[END_REF][START_REF] Vissers | New CO-CO interaction potential tested by rovibrational calculations[END_REF]. If, however, the monomer wavefunction is only poorly described by a single reference wavefunction, this approach becomes again difficult to be carried out. On the other hand, treating the interaction of several molecules with significant non-dynamical correlation in a super-molecular approach (∆E = E dimer -E monomers ) leads immediately to very large active spaces, and thus to calculations of immense computational effort.

Here, the VB approach may be useful in order to reduce this effort by replacing the multireference expansion by only a few configurations for each monomer, and by calculating directly the contribution of each pair of configurations to the overall interaction energy, without the need for a super-molecular set-up. In the SAPT development, non-orthogonality of orbitals is as well dealt with, as the monomer wavefunctions are obtained completely independently. As well antisymmetrization is added in a second step in the SAPT expansion ("perturb first -anti-symmetrize later" ansatz [START_REF] Jeziorski | Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes[END_REF]). In practice, however, a dimer-centred basis set is generally used to expand the molecular orbitals of the monomers, and in the hybrid SAPT approach the infinite series of induction contributions is obtained via a Hartree-Fock calculation on the complete dimer system, abandoning partially the approach of completely independent monomer calculations.

In the present manuscript we would like to explore two terms of the SAPT development: the electrostatic energy E 100 pol and the second-order dispersion energy E 200 disp . Of course, for a complete description one needs to approximate at the same level of precision the other SAPT terms as well, notably the first-order exchange and induction. This we will leave for further investigations.

Having set the frame, we will use the dimer of nitrogen protoxide (N 2 O) as example. As for most of the nitrogen oxides the electron distribution is rather complicated due to the relatively high local density -in a chemist's view the concurrence of triple and double bonds around a central nitrogen atom [START_REF] Harcourt | Valence Bond Structures for N 2 O and HCNO with Apparently Pentavalent Nitrogen Atoms[END_REF], to be cast into orbitals constructed from the four hybrid lobes of the nitrogen 2sp valence space.

In the following we will first develop the theoretical approach, then give some relevant data for the N 2 O single molecule before discussing the results on the dimer. All purely technical details will be collected in an appendix.

Theoretical aspects

The VB development of the wavefunction

Modern Valence Bond theory seemed to be reserved to a few experts in the field. However, with the inclusion of VBSCF in larger program packages and notably the development of the Xiamen VB program [START_REF] Song | XMVB 2.1: An ab initio Non-Orthogonal Valence Bond Program[END_REF], calculations can be carried out more routinely nowadays.

The basic idea is that the wavefunction is developed in several contributions, distinguished by orbital occupations and Rumer-like spin couplings [START_REF] Rumer | Zur Theorie der Spinvalenz[END_REF] of orbital pairs leading to an overall singlet wavefunction:

Ψ VB = I cI ΦI ΦI = A {(Closed)(Active)}
where (Closed) designates the doubly occupied closed-shell orbitals φ c i φc i and (Active) the singletcoupled pairs of singly occupied orbitals φ 1 φ2 -φ1 φ 2 / √ 2. The anti-symmetrizer A generates from this a completely anti-symmetric sum of 2 n determinants for n unpaired electrons in a configuration. Even if the VB wavefunction is written in configurations Φ, we decompose it into simple determinants Φ I . In the same way decompose the coefficients cI into coefficients of determinants c I .

The wavefunction is normalized

1 = Ψ|Ψ = IJ c I c J Φ I |Φ J (1)
needing thus explicitly the overlap matrix elements between determinants -as well as Rumer diagrams themselves are not orthogonal.

From the expansion in determinants a first-order (or one-particle) density matrix in atomic orbitals can be obtained, and, of course, a second-order matrix as well. From the first-order density matrix of either monomer in the intermolecular set-up the electrostatic interaction energy is obtained as

E 100,V B pol = I∈1,J∈B Z I Z J | R I -R J | - α,β∈A P A,V B αβ χ α | I∈B 1 | r -R I | |χ β - γ,δ∈B P B,V B γδ χ γ | I∈A 1 | r -R I | |χ δ + α,β∈A P A,V B αβ γδ∈B P B,V B γδ (αβ|γδ) (2)
with the density matrices P V B and mono-and bielectronic integrals in atomic orbitals, in the same way as for instance in the SAPT approaches based on a density-functional description of the monomers [START_REF] Vissers | New CO-CO interaction potential tested by rovibrational calculations[END_REF][START_REF] Hesselmann | Density-functional theory-symmetry-adapted intermolecular perturbation theory with density fitting: A new efficient method to study intermolecular interaction energies[END_REF]. As the VB expansion includes already some of the intramolecular correlation energy, the calculated term should lie in between the bare E 100 pol SAPT term with Hartree-Fock density matrices and the infinite sum of including the intramolecular perturbation series ǫ (1) pol = i j≥0,i+ j>0 E 1i j pol .

When addressing the dispersion energy, we will rather use each determinant separately instead of the density matrix. This means that we should carry along the expansion coefficients c I or pairs of expansion coefficients from the normalization of the wavefunction ψ as in equation [START_REF] Heitler | Wechselwirkung neutraler Atome une homöopolare Bindung nach der Quantenmechanik[END_REF]. This is all but convenient; we resort therefore to (Coulson-Chirgwin) weights instead of products of coefficients, in order to linearize the normalization condition, in the spirit of a Mulliken analysis [START_REF] Wu | Classical Valence Bond Approach by Modern Methods[END_REF]:

1 = I c I J c I Φ I |Φ J w I = I w I (3)
A product of contributions from the determinants of monomer A (like Φ I |M|Φ I ) and the determinants of monomer B (like Φ J |M ′ |Φ J ), as we will need it later, will reduce then to a double sum

I∈A J∈B w I w J M I M ′ J .
In this context it is irrelevant whether the molecular orbitals spanning each determinant are common to all determinants of one monomer or not. A first VBSCF calculation yields one set of reasonable, localized orbitals and furnishes an estimation of relative weights of the different contributions to the wavefunction.

In second step, the BOVB step, the orbitals are varied for each configuration separately, minimizing the total energy of the overall wavefunction. From the BOVB calculation we obtain a global one-particle density matrix, including all determinants, and a list of individual determinants with weights and individual sets of orbitals with double and with single occupation.

The dispersion term

The key element of the present paper is the estimation of the dispersion interaction, the most eminent contribution to intermolecular interactions with a purely quantum mechanical origin. In textbooks one may find interacting, fluctuating dipoles for these interactions, via the simple MP2 formulation in spin orbitals

E MP2 = 1 4 ia jb ((ia| jb) -(ib| ja)) 2 ǫ i + ǫ j -ǫ a -ǫ b (4)
or, for a closed-shell case in spatial orbitals

E MP2 = ia jb (ia| jb) (2(ia| jb) -(ib| ja)) ǫ i + ǫ j -ǫ a -ǫ b = 4 i< j,a<b (ia| jb) 2 + (ib| ja) 2 -(ia| jb)(ib| ja) ǫ i + ǫ j -ǫ a -ǫ b ( 5
)
where we have the dispersion terms when having orbitals i and a on monomer A, and orbitals j and b on monomer B. In localized orbitals this sum can be divided into six classes, as pointed out already a while ago [START_REF] Runeberg | The aurophilic attraction as interpreted by local correlation methods[END_REF][START_REF] Reinhardt | The decomposition of intermolecular interaction energies in localized orbitals -critical analysis and an invariance[END_REF] and as displayed in Figure 1. The class of intra-molecular excitations is not meant to contribute to the interaction energy, and classes describing charge polarization, spin polarization and the BSSE term will become exponentially small with increasing intermolecular distance. At the end only the contribution

E inter MP2 ≈ E disp MP2 = 4 i,a∈A j,b∈B (ia| jb) 2 ǫ i + ǫ j -ǫ a -ǫ b (6)
will remain with a significant contribution to the interaction energy. If we deal as well with singly occupied orbitals, the factor of four becomes a sum over the different contributions for different spin settings, as and sums run over spin-orbitals, in principle different for each spin. Developing the interaction of the two overlap densities ρ

E disp MP2 = σ A i σ A ,a σ A ∈A σ B j σ B ,b σ B ∈B (i σ A a σ A | j σ B b σ B ) 2 ǫ i σ A + ǫ j σ B -ǫ a σ A -ǫ b σ B (7) B A B A B A B A BSSE B A B A
A ( r) = φ i σ A ( r) φ a σ A ( r) and ρ b ( r) = φ j σ B ( r) φ b σ B ( r) into interactions of transition multipoles at respective centres R of A and B M R i jk = ρ( r) (x -X) i (y -Y) j (z -Z) k d 3 r = ρ( r + R) x i y j z k d 3 r
we may use the fact that i and a are orthogonal within one Slater determinant (the "transition charge"). Thus, the first non-zero multipole is the (transition) dipole µ i→a = φ i | r | φ a , which we may situate for the multipolar development on the average of the positions of the occupied and the virtual orbital r

µ i→a = φ i | r | φ i + φ a | r | φ a /2.
This leaves us then with the classical dipole-dipole interactions

E dip-dip = µ A • µ B R 3 AB -3 ( µ A • R AB )( µ B • R AB ) R 5 AB ( 8
)
as a model for the bielectronic integral of the dispersion contribution to the intermolecular interaction energy. By replacing the bielectronic integral in Eq. 7 the global formula to evaluate the dispersion interaction energy between two molecules in he VB context becomes thus a double sum over all determinants with respective weights, orbitals and spins as

E disp VB ≈ I∈A w I σ A i σ A ,a σ A ∈A J∈B w J σ B j σ B ,b σ B ∈B (E dip-dip (i σ A , a σ A , j σ B , b σ B )) 2 ǫ i σ A + ǫ j σ B -ǫ a σ A -ǫ b σ B ( 9
)
All quantities are evaluated strictly for each monomer separately, maintaining thus the spirit of SAPT and VB calculations. Once the dipoles and their respective positions determined, the only intermolecular quantity to be evaluated are distance vectors and distances, allowing for a rapid summation over the different terms. The result for the H 2 dimer in a parallel setup, in the minimal STO-9G basis set. The dipole approximation reproduces exactly the SAPT or MP2 results. For short distances orthogonality and penetration give the deviation between dipole/SAPT and the MP2 energy. Note that the energy scale is logarithmic. Taking as simplest possible example two hydrogen molecules in a minimal STO-9G basis set and a parallel set-up, we may compare the MP2 dispersion -4

(σ 1 σ * 1 |σ 1 σ * 1 ) 2 2(ǫ σ * -ǫ σ )
and the SAPT

E 200

Disp term with our dipole model, at different distances. As expected we recover the SAPT result at not-too-small intermolecular distances, see Figure 2.

We may try thus our model for a more non-trivial intermolecular system as presented by the N 2 O dimer.

The N 2 O dimer

For the interaction of two linear molecules the relative geometries may be defined in several ways. We will use a general set-up, in two local coordinate systems the origins of which being placed at the respective barycenters aligned with the x axis. The overall barycentre is situated at the global origin. With respect to the global frame, the local coordinate systems may be rotated around the three Euler angles α, θ and φ, for successive rotations around the z, the y and again the z axis [START_REF] Sherman | Transformation matrices fo the rotation of real p, d and f atomic orbitals[END_REF].

The N 2 O dimer has been well studied experimentally with matrix spectroscopy, see for instance Hecker et al. [START_REF] Hecker | IR-High resolution spectroscopy of the N 2 O dimer: The torsional mode[END_REF]. The geometry of the ground state shows a slipped anti-parallel set-up with a distance of 3.4225 Å between the respective barycenters, and an angle θ 1 of -30.0 degrees (and θ 2 = 120 degrees, see Figure 3) [START_REF] Hecker | IR-High resolution spectroscopy of the N 2 O dimer: The torsional mode[END_REF][START_REF] Huang | The structure of the nitrous oxide dimer from sub-Doppler resolution infrared spectroscopy[END_REF][START_REF] Valdés | The N 2 O•N 2 O, N 2 O•SO 2 , and (N 2 O) 2 •SO 2 van der Waals Complexes: An ab Initio Theoretical Analysis[END_REF][START_REF] Glendening | Ab initio calculations of nitrogen oxide reactions: Formation of N 2 O 2 , N 2 O 3 , N 2 O 4 , N 2 O 5 , and N 4 O 2 from NO, NO 2 , NO 3[END_REF][START_REF] Zheng | Theoretical studies of the N 2 O van der Waals dimer: Ab initio potential energy surface, intermolecular vibrations and rotational transition frequencies[END_REF][START_REF] Berner | Nitrous oxide dimer: An ab initio coupled-cluster study of isomers, interconversions, and infrared fundamental bands, and experimental observation of a new fundamental for the polar isomer[END_REF], with the two oxygen atoms lying close to each other. Corroborated by calculation [START_REF] Glendening | Ab initio calculations of nitrogen oxide reactions: Formation of N 2 O 2 , N 2 O 3 , N 2 O 4 , N 2 O 5 , and N 4 O 2 from NO, NO 2 , NO 3[END_REF][START_REF] Zheng | Theoretical studies of the N 2 O van der Waals dimer: Ab initio potential energy surface, intermolecular vibrations and rotational transition frequencies[END_REF], this is consistent with an interaction governed by dipole (favouring anti-parallel) and quadrupole (favouring slipped) electrostatic interactions.

The theoretical study of Glendening and Halpern [START_REF] Glendening | Ab initio calculations of nitrogen oxide reactions: Formation of N 2 O 2 , N 2 O 3 , N 2 O 4 , N 2 O 5 , and N 4 O 2 from NO, NO 2 , NO 3[END_REF] documents the difficulties to obtain reliable potential surfaces with single-reference methods [START_REF] Valdés | The N 2 O•N 2 O, N 2 O•SO 2 , and (N 2 O) 2 •SO 2 van der Waals Complexes: An ab Initio Theoretical Analysis[END_REF], and the need for a multi-reference treatment right from the beginning. Looking at the T1 term in a Coupled-Cluster calculation (CCSD), as suggested by Lee and Taylor [START_REF] Lee | A diagnostic for determining the quality of single-reference electron correlation methods[END_REF], leads to a value of 0.02, close to the borderline indicating the possible inadequacy of a single-reference calculation.

With a planar ground-state geometry obtained by all cited reports on the N 2 O dimer, we will place all atoms in one single plane which is described by the distance of the barycenters and the two rotation angles around the respective y axes (see Figure 3) and restrict the geometry variations to three possibilities: (1) rotating one monomer only whilst the other is held fixed along the z axis (θ 1 = 0, 0 ≤ θ 2 ≤ 2π), (2) rotating both monomers in parallel (0 ≤ θ 1 = θ 2 ≤ π/2), and (3) rotating both monomers anti-parallel (π/2 ≤ θ 2 = π + θ 1 ≤ π/2). We will denote these three variants as "T", "P" and "AP", respectively.

Actual calculations

We are interested whether the proposed idea may work out to be useful, not for quantitative results. For being able to carry out the necessary calculations in a reasonable time, we restrict ourselves to a "standard" aug-cc-pvdz basis set. For the reference SAPT calculations we use a dimer-centred basis set with 5 d orbitals, whereas for VB calculations, which are conducted only for one isolated monomer, the orbitals are expanded in a cartesian basis with 6 d orbitals, including the spherically symmetric x 2 + y 2 + z 2 ), employing MP2 in an aug-cc-pvdz or aug-cc-pvtz basis set, we find N-N bond distances of 1.174 and 1.155 Å, and longer N-O bond distances of 1.190 and 1.181 Å, respectively. The N-N bond length is obtained significantly longer than the experimentally documented one [START_REF] Herzberg | Electronic spectra and electronic structure of polyatomic molecules[END_REF] of 1.128 Å, whereas the other one is close to the experimental data of 1.184 Å. CCSD(T) yields about the same results with 1.149 and 1.202 Å, respectively, as reported already by Glendening and Halpern [START_REF] Glendening | Ab initio calculations of nitrogen oxide reactions: Formation of N 2 O 2 , N 2 O 3 , N 2 O 4 , N 2 O 5 , and N 4 O 2 from NO, NO 2 , NO 3[END_REF]. If we optimize the geometry within a MCSCF approach, we find a better agreement, as of 1.133 and 1.190 Å for the two bonds in an aug-cc-pvtz basis set. This deviation for MP2 or CCSD(T) results may point already to an inherent difficulty to treat the N 2 O molecule correctly with a single-reference method.

The MCSCF calculation (8 electrons in 6 active π orbitals) gives indeed a main configuration with a coefficient of 0.94, two closed-shell configurations with a coefficient of -0.13, and four configurations with four open shells with coefficients of about 0.08. Table 1 shows the respective total energies for the different calculations, once for a near-minimal 6-31G* basis set and once for the basis set employed throughout. If we extend the MCSCF wavefunction to the full valence space, i.e. 16 electrons in 12 orbitals, the expansion runs over 61677 determinants in C 2v symmetry and about 250 000 independent variables, and the subsequent ACPF (Averaged Coupled Pair Functional) calculation involves, in the aug-cc-pvdz basis set, a total of about 167 millions of uncontracted intermediate configurations.

To these we compare the total energies of the VB-SCF and BOVB calculations, which we will detail in the next section.

VB calculation on N 2 O

Harcourt and Schulz [START_REF] Harcourt | Valence Bond Structures for N 2 O and HCNO with Apparently Pentavalent Nitrogen Atoms[END_REF], and later Braïda et al. [START_REF] Braïda | A Clear Correlation between the Diradical Character of 1,3-dipoles and their Reactivity toward Ethylene or Acetylene[END_REF] published VB calculations on several 1,3-dipoles which as well comprise the NNO molecule. In both publications the VB structures of the π systems were elucidated, putting forward that open-shell VB structures are the most important ones in the expansion of the wavefunction. All calculations for these more pedagogical articles were carried out in minimal (STO-6G) or near-minimal (6-31G*) basis sets. As we would like to make interact two molecules, we push the VB calculations farther by using the aug-cc-pvdz basis. To this end we define as contributing orbitals the respective p orbitals of the three atoms, which makes thus again 8 electrons in 6 orbitals to distribute as in the previous MCSCF calculations.

In total the valence-bond program used, xmvb, creates 105 possible structures or 225 determinants. Open shells are spin-adapted in a Rumer scheme, coupling pairs of orbitals in singlets. In our case we may have at most four open shells, making it thus 4 determinants for produc-

ing 1 2 (α 1 β 2 -β 1 α 2 ) (α 3 β 4 -β 3 α 4 ) = 1 2 (α 1 β 2 α 3 β 4 -β 1 α 2 α 3 β 4 -α 1 β 2 β 3 α 4 + β 1 α 2 β 3 α 4 ) or 2 for having 1 √ 2 (α 1 β 2 -β 1 α 2 )
With the trivial expansion of a doubly occupied molecular orbital on two sites

σ σ = (s 1 + s 2 )( s1 + s2 ) = s 1 s1 + s 1 s2 + s 2 s1 + s 2 s2
we see that the open-shell VB contributions in atomic p orbitals may be important whereas the closed-shell ionic components may seem negligible, in apparent, but obvious contrast to the MCSCF results.

From the 105 generated structures we have in a simple VB calculation (all orbitals are common to all structures) about half of them with coefficient > 10 -4 (precisely 57), and all remaining structures are found with weights < 10 -15 due to numerical precision and symmetry. The numbers of 105 and 57 correspond as well to the number of configurations generated when restricting the MCSCF space to 8 electrons in 6 orbitals, without and with symmetry respected in the calculation. Of course, 57 structures do not allow simple interpretations -however, we can see some trends. Of the 57 structures we find 45 with evenly distributed electrons in x and y, the other 12 structures having 6 electrons in one direction and 2 in the other.

The most important VB configuration with representing 4 determinants, is that of doubly occupied p x and p y on oxygen, and otherwise singly occupied p orbitals of the nitrogen atoms, in perfect agreement with chemical intuition, see Figure 4. The next step is to set up a BOVB calculation, in order to add some (more) dynamical correlation. Looking at the relative importance of the coefficients, we may use the 45 VB structures having 4 electrons in the π x and 4 electrons in the π y system, as the remaining 12 ones drop to zero in their relative importance in the BOVB expansion. Among the 45 configurations are three configurations with closed shells and completely ionic atoms, which seem non-physical, but are necessary for generating the MO picture with completely delocalized closed-shell orbitals. Therefore the coefficients of these configurations are not at all vanishing.

Overlap of structures is maximally 67 %, showing convergence to a reasonable minimum. We see some small and negative Coulson-Chirgwin weights, astonishing at first glance, but numerically possible due to the Mulliken-like symmetric distribution of cross terms.

The BOVB calculation in the aug-cc-pvdz basis set includes about 13500 independent variables to optimize recovering 0.21 a.u. of correlation energy, which is about 40 % with respect to the MP2 result (see Table 1). We see as well that the converged BOVB energy is about the same as that from the MCSCF calculation in the full valence space, even though the BOVB function includes only determinants with active p orbitals. Relaxing orbitals for different occupations included thus some electron correlation, and generates a description of the full valence space by significant configurations only.

We still have a singlet wavefunction which allows to construct an average Fock matrix as in Dewar's half-electron method [START_REF] Dewar | Ground state of conjugated molecules. IX. Hydrocarbon Radicals and Radical Ions[END_REF]. However, as the active orbitals are restricted on the respective atoms, they are non-orthogonal with respect to each other and with respect to the doubly occupied orbitals. An open-shell determinant resembles thus more an UHF wavefunction, as the α and the β part span different orbital spaces. Therefore we have to generate separately for each determinant of each configuration or Rumer diagram an α and a β orbital space and corresponding Fock matrices in the atomic-orbital basis (see appendix).

Excitations

As we are interested in intermolecular interactions, we will ignore terms of spin polarization between the two moieties. In our terms this means that we look only at mono-excitations conserving the spin on each monomer. Since we are in an UHF context, it is quite improbable that a zero-gap situation arrives. Occupied and virtual orbitals should be separated energetically for each spin type.

Each excitation creates in our picture a transition dipole. As the calculation runs only over one single monomer, the ensemble of dipoles for the interaction geometry may be easily obtained not by rotating the whole molecule and the VB calculation (which will make the definition of the active space difficult), but only the transition dipoles, and their respective localization in space with respect to the molecular barycentre. The energy denominators remain, of course, invariant.

If we used eq. 7, we would have to calculate all bielectronic integrals of the dimer system, first in atomic orbitals and the four-index-transformed over molecular orbitals, for each pair of determinants, and for each relative orientation and distance.

The N 2 O dimer

Having a density matrix and excitations (dipoles) at hand for one molecule, we may assemble these data for the dimer by rotating the local frames, recalculating mono-and bi-electronic integrals to be contracted with the density matrices, and set up the SAPT calculations for the different geometries, as these were carried out in a dimer-centred basis set. pol interaction (left), and the difference of VB density interactions to the SAPT E (1) pol term (circles), the "response" term E (1) pol,r (triangles), and the bare E100 interaction term (squares), for the T geometry variation. Full symbols are for an intermolecular distance of 7.0 bohr, whereas hollow symbols represent the corresponding calculations for 10 bohr.

Electrostatic interactions

For having an estimation of the effect of the different basis sets we compare first the Hartree-Fock electrostatic interaction of the two unrelaxed monomers to the SAPT E (100) pol term, for the "T" set-up with one molecule rotated with respect to the other (Figure 5, left panel). The effect of the slight differences in basis sets is indeed negligible with about 0.1 kcal/mol at the intermolecular distance of 7 bohr (3.69 Å, close to the experimentally observed equilibrium distance). On the other panel we trace the difference of the corresponding VB density interaction to the HF density interactions and the E (1) pol energy, i.e. including intra-monomer correlation, in the present case up to second order.

The intra-molecular correlation can be added to the SAPT polarization term in two different ways as the MP2 density is not a uniquely defined object. Either the expectation value may be calculated as a sum over the contributions of all determinants of the first-order wavefunction, or the response of the density to an infinitely small external charge may be considered [START_REF] Jeziorski | Perturbation Theory Approach to Intermolecular Potential Energy Surfaces of van der Waals Complexes[END_REF]. The former will be denoted with E (1) pol , and the latter with E (1) pol,r , both including second-and third-order intra-molecular corrections.

Augmenting slightly the intermolecular distance from 7 to 10 bohr reduces consistently the differences between the different methods to evaluate the electrostatic interaction of the monomers.

In Figure 6 we give then the absolute energies, as well for corresponding MCSCF densities, for all three geometry variations. As found for the comparison of total energies, the MCSCF and the VB electrostatic interactions are very close, pointing again toward the idea that with the restricted expansion into 45 configurations with non-orthogonal and individual orbitals we arrive at a wavefunction of the same quality as for the MCSCF expansion with its about 60 000 determinants.

The dispersion terms

Having generated UHF-like orbitals for each VB determinant and spin, we obtain excitations and and their location for the N 2 O monomer as depicted in Figure 7 for the single-determinant Hartree-Fock case. We see several orbital barycenters of virtual orbitals, and dipole centres, far beyond the nitrogen end of the molecule, whereas at the opposite end orbitals and transition dipoles lie close to the oxygen nucleus in the HF case, but not for the VB configurations.

This indicates that our dipole model will behave differently for a parallel rotation and an anti-parallel rotation for fixed intermolecular distance as some dipoles will come closer when the two N ends point to each other than when the oxygen ends meet, even if numerators may remain small.

For the inter-molecular separation of 7 bohr we trace in Figure 8 the different approaches to the dispersion interaction, for the three geometry variations. For following the distance dependence of the calculated dispersion terms, we collect in Figure 9 for selected angles and distances in a double-logarithmic plot, for highlighting the power law of our model, and the dispersion term calculated with the other methods. Indeed, we observe in the limit a close coincidence of all approaches, but only our model follows the ideally expected R -6 power law. In Figure 10 we give at last the difference of the dispersion models with respect to our VB-dipole-approximation. We do not distinguish explicitly the three geometry variations, and of course, for 90 degrees the differences become much larger than for 30 degrees, where the ends of the molecules remain well separated even at moderate intermolecular distances.

Discussion

The electrostatic term

From the figures in section 3.3 we recognize that the SAPT E 100 pol term and our HF model are, by construction the same, as expected, modulo the slight basis set difference due to the use of 5 or 6 d functions, and the monomer-or dimer-centred basis sets. Including intra-molecular correlation in the SAPT approach via the two possible density models (expectation value or response) results in two opposite trends. Including intra-molecular correlation via the valencebond approach or a multi-configurational calculation gives about the same effect, and lies in the present case on average closer to the response SAPT calculation than the evaluation of the density expectation value with the first-and 2nd-order wavefunctions. The VB approach gives thus a unique answer with a shorter expansion in determinants as the comparable MCSCF wavefunction. However, a MCSCF calculation in orthogonal orbitals still remains computationally advantageous with respect to the BOVB approach if only a first-order density matrix is required.

The dispersion term

Concerning our model for the dispersion energy, we see from Figure 8 that globally some trends are reproduced, in particular the asymmetry of the angular dependence for the "T" geometry variation. If we keep in mind that we model a bi-electronic integral (ia| jb) by two interacting dipoles φ i | r | φ a and φ j | r | φ b , we see that the bielectronic integral has a different behaviour with respect to orbital symmetry than the dipoles. For instance for a pair ia of orbitals of the same symmetry, the corresponding transition dipole may become exactly zero, whereas the complete integral does not at all vanish if the jb part does not disappear. Transition quadrupoles and higher moments may as well be non-zero, and contribute thus to the expansion of the bielectronic integral if the components lie in the irreducible representation of totally symmetric objects. Thus the reproduction of the trends of the angular dependence of the calculated dispersion interaction, despite this inherent systematic difference of approaches, is already a result to be noted. Using dipoles from a Hartree-Fock calculation or from a valence-bond calculation seems, alas, of secondary importance, or result in a systematic deviation from the SAPT values.

The dipole model follows strictly, by construction, a 1/R 6 power law, as seen from Figure 9, whereas SAPT and the MP2 dispersion in localized orbitals deviate from this power law for smaller distances. In particular the significant inflection of the MP2 dispersion is due to the use of orthogonal dimer orbitals, difficult to remain sufficiently localized for representing monomer wavefunctions. We may, on the other hand, extract reasonable asymptotic informations like C 6 coefficients, even for intermolecular distances where the other methods do not yet give a simple power law.

Figure 10 shows again that the model in the present form is not, alas, reliable for interesting situations.

Conclusions

In the present contribution we explored the use of transition dipoles, generated from the different configurations in a Breathing-Orbital Valence-Bond context for describing some elements of intermolecular interactions. The non-orthogonality of orbitals belonging to different configurations is circumvented by using the underlying determinants one by one. In that sense the VB approach may present an alternative to recent developments in symmetry-adapted perturbation theory, even including non-dynamic correlation of the monomer wavefunctions. Indeed, for the electrostatic interaction energy results lie close to those employing MCSCF densities of the monomers, with a much more compact form of the wavefunction. Were it not that the optimization of the BOVB wavefunction, even in the moderate basis set used here for intermolecular interactions, demands already a computational effort beyond routine applications.

The dipole model to the dispersion energy reproduces qualitatively the trends of other ways to estimate this interaction component, and coincides asymptotically with SAPT and the MP2 estimation for large intermolecular distances, despite its inherently wrong angular dependence for reproducing a bielectronic integral by the interaction of transition dipoles. For a quantitative reproduction in the region of interest for studying potential surfaces close to an equilibrium geometry, our model reveals to be too crude an approximation, mostly due to the eccentric positions of the transition dipoles, exaggerating thus dipole-dipole interactions through very short distances.

For the interaction of two N 2 O molecules the dispersion term is not significantly modified by the multi-reference character of the monomer wavefunction. Including thus intra-molecular correlation in the SAPT second-order dispersion term makes not a qualitative difference. Nevertheless, our dipole model comes closer to the SAPT results when employing the VB-generated orbitals in the multi-determinantal ansatz. The electrostatic term confirms that the mostly closedshell MCSCF wavefunction in delocalized molecular orbitals is of the same quality as the VB wavefunction expanded in mainly open-shell configurations in localized orbitals.

For further investigations the more complicated exchange contributions to the SAPT expansion may be studied towards a more complete description of intermolecular interaction between multi-configurational monomer wavefunctions.

Transformation of VB information to excitations

Several interface programs have been written, in order to transfer information between the different codes, using a Fortran preprocessor [START_REF] Scemama | IRPF90: a programming environment for high performance computing[END_REF].

From the BOVB run of xmvb [START_REF] Song | XMVB 2.1: An ab initio Non-Orthogonal Valence Bond Program[END_REF] we take from the output the section containing the structures, and a second section containing the orbitals. From this information we create for each structure an orbital file with occupation numbers and the expansion coefficients of the corresponding orbitals. The density is written already to file by xmvb.

The geometry, orbitals and density matrix are translated and rotated for the positions of the two monomers, in order to calculate the SAPT interaction energy, and to recalculate the electrostatic interaction directly from the density matrix.

From the positions and the basis sets we generate input for the ORTHO suite of experimental ab-initio programs [START_REF] Reinhardt | Ortho, series of ab-initio programs in localized orbitals[END_REF] -here in order to calculate mono-electronic and bi-electronic integrals as well as dipole-and overlap integrals. The actual calculation of integrals is again left to Dalton [START_REF] Helgaker | Dalton, an ab initio electronic structure program[END_REF], but may be carried out as well with other integral generators.

Once the integrals and orbitals at hand, the unrestricted Fock matrices are constructed and diagonalized, leading to canonical orbitals, orbital energies and transition dipole moments, which are collected for each determinant separately.

Generating UHF molecular orbitals for the monomers

In UHF a separate Fock matrix is constructed for α and β spins [START_REF] Berthier | Extension de la méthode du champ moléculaire self-consistent à l'étude des états à couches incomplètes[END_REF]. From the two orbital sets φ α and φ β we construct two density matrices

F φ i = ǫ i φ i
We may then diagonalize each Fock matrix separately, for obtaining orbitals and corresponding orbital energies, even if the VB orbitals are not self-consistent solutions of the UHF equations (or Pople-Nesbet-Berthier equations). In order to keep the orbital spaces invariant, despite the diagonalization of the respective Fock matrices, elements coupling occupied and virtual orbitals are set to zero, i.e. the Fock matrices are diagonalized only within the occupied or virtual orbitals without any mixing. Of course, orbital energies are then not strictly ordered: energies belonging to occupied orbitals may appear higher in energy than those corresponding to unoccupied ones. 19 4) (hollow circles). In the upper panels only the difference of the latter with the former is traced.

Monomer-localized dimer orbitals via CI of single excitations

In order to evaluate the MP2 energy in localized orbitals, we first perform a Hartree-Fock calculation on each monomer separately, generating canonical orbitals. These serve as starting orbitals for the dimer step, where we use the iteration on the coefficients of the configuration interaction of singly excited determinants to correct occupied and virtual orbitals [START_REF] Daudey | Direct determination of localized SCF orbitals[END_REF][START_REF] Reinhardt | Direct Determination of Localized Hartree Fock Orbitals as a step towards N scaling procedures[END_REF], maintaining the localization of the orbitals spaces (occupied and virtual) on the respective monomers as much as possible. The MP2 energy may be evaluated using equation (4), neglecting thus the effect of non-zero off-diagonal Fock matrix elements. For including these we minimize the Hylleraas functional [START_REF] Hylleraas | Über den Grundterm der Zweielektronenprobleme von H -, He, Li + , Be ++ usw[END_REF], leading to CEPA-like equations with the Fock matrix instead of the full Hamilton matrix. The correlation energy is evaluated at the end as a sum over amplitudes t ab i j multiplied with appropriate (anti-symmetrized) bi-electronic integrals E Corr = i jab t ab i j (ia| jb) [START_REF] Harcourt | Valence Bond Structures for N 2 O and HCNO with Apparently Pentavalent Nitrogen Atoms[END_REF] which again may be decomposed into the six excitation classes mentioned above (Figure 1). As we started with canonical monomer orbitals which are only slightly perturbed during the dimer SCF procedure, the impact of the remaining Fock-matrix elements is expected to be small. Computer time, however, is significantly increased when solving the linear system of equations instead of the straight-forward evaluation of the perturbation terms. For two situations (the T geometry variation at R = 7.0 a.u. and different angles, and for 90 degrees and variable R) we show in Figure 11 the difference of the two variants.

-exploring non-dynamical correlation for intermolecular interactions -using non-orthogonal molecular orbitals -simple dipole-dipole interaction for dispersion -development of SAPT for MCSCF functions Highlights

Figure 1 :

 1 Figure 1: The 6 classes of di-excitations in the MP2 correlation energy of the dimer system.

Figure 2 :

 2 Figure2: The result for the H 2 dimer in a parallel setup, in the minimal STO-9G basis set. The dipole approximation reproduces exactly the SAPT or MP2 results. For short distances orthogonality and penetration give the deviation between dipole/SAPT and the MP2 energy. Note that the energy scale is logarithmic.

Figure 3 :

 3 Figure3: Definition of the local frames and angles. Note that the angle is measured differently in the different publications. Our angle of 30 degrees corresponds to 60 degrees in Havenith et al[START_REF] Hecker | IR-High resolution spectroscopy of the N 2 O dimer: The torsional mode[END_REF].

Figure 4 :

 4 Figure 4: The distribution of 8 electrons among 6 atom-centred p orbitals results in a leading contribution of 4 singly occupied orbitals on the nitrogen p orbitals.

Figure 5 :

 5 Figure5: Difference of the interactions of the Hartree-Fock densities and the SAPT E(100) pol interaction (left), and the difference of VB density interactions to the SAPT E(1) pol term (circles), the "response" term E(1) pol,r (triangles), and the bare E100 interaction term (squares), for the T geometry variation. Full symbols are for an intermolecular distance of 7.0 bohr, whereas hollow symbols represent the corresponding calculations for 10 bohr.

Figure 6 :Figure 7 :

 67 Figure6: The electrostatic interaction (in kcal/mol) for the three series (rotating molecules in parallel, P, anti-parallel, AP, and only one of the two, T), at a fixed intermolecular distance of R = 7.0 bohr. Displayed are the interaction of the unrelaxed Hartree-Fock densities or E100 SAPT (full circles), SAPT term including intra-molecular correlation (full and hollow triangle point upwards) and the interaction of correlated monomer densities (full and hollow triangles pointing down, for the VB wavefunction and for the MCSCF wavefunction)

Figure 8 :

 8 Figure8: Angular dependence of the calculated dispersion energy, for the two SAPT approaches (with and without intra-molecular correlation terms, full and hollow circles), for our model with HF orbitals (upper triangles) and with the valence-bond-generated orbitals (diamonds), and the MP2 dispersion term in localized orbitals (lower triangles). All calculations are done for an intermolecular interaction of 7 bohr.

Figure 9 :

 9 Figure 9: The distance dependence of the different models, for selected angles and the three geometry variations, in a double-logarithmic plot. The legend is the same as for the previous figure.

Figure 10 :

 10 Figure 10: Energy difference with respect to our VB model, as a function of the intermolecular distance, for the three geometry variations and different angles. Symbols correspond to the legends for the previous two figures.

  The two Fock matrices are coupled through the electrostatic interaction, leading toF α ρσ = h ρσ + k∈α ((ρσ|kk) -(ρk|kσ)) + k∈β (ρk|kσ) = h ρσ + γδ P α γδ ((ρσ|γδ) -(ργ|δσ)) + γδ P β γδ (ρσ|γδ) =Ξ α ρσ F β ρσ = h ρσ + k∈β ((ρσ|kk) -(ρk|kσ)) + k∈α

Figure 11 :

 11 Figure 11: Comparison of the dispersion terms calculated via the Hylleraas functional (full circles) and the straightforward application of equation (4) (hollow circles). In the upper panels only the difference of the latter with the former is traced.

Table 1 :

 1 Total energies for the different calculations in a nearly-minimal 6-31G* and in the aug-cc-pvdz basis set. MCSCF (8/6) and (16/12) stand for multi-reference calculations with, respectively, 8 electrons in 6 orbitals or 16 electrons in 12 orbitals. All calculations are carried out at the experimental geometry.

	function.

3.1. The N 2 O monomer 3.1.1. "Standard" calculations for the monomer If we optimize the linear geometry of a single N 2 O molecule (NNO

Acknowledgments

The present work has been carried out in the Laboratoire de Chimie Théorique, University Pierre and Marie Curie in Paris, and UMR 7616 of the CNRS. No other grants were involved. Discussions with B. Braïda (Paris) and his help on the VB program are gratefully acknowledged, as well as the exchange with E. Alikhani (Paris) and A.L. Tchougréeff (Aachen, Moscow) on experimental and theoretical data.

Technical details

The VB calculations are carried out with the program xmvb [START_REF] Song | XMVB 2.1: An ab initio Non-Orthogonal Valence Bond Program[END_REF], version 2.1, in a 6-31G* and an aug-cc-pvdz basis set. The starting orbitals are generated with the preint tool of xmvb.

The version of 1996 of SAPT [START_REF] Jeziorski | SAPT: A Program for Many-Body Symmetry-Adapted Perturbation Theory Calculations of Intermolecular Interaction Energies[END_REF] is coupled to the Dalton program package [START_REF] Helgaker | Dalton, an ab initio electronic structure program[END_REF] which serves as integral generator. For the MCSCF calculations we employed Molpro, as of version 2008.2 [START_REF] Werner | MOLPRO is a package of ab-initio programs written[END_REF].

The p orbitals on each atom are taken as active orbitals, defining the σ skeleton as closed shells. We have thus a space of 6 atom-centered localized orbitals, carrying 8 electrons.