# On the characterization of flowering curves using Gaussian mixture models

Frédéric Proïa, Alix Pernet, Tatiana Thouroude, Gilles Michel, Jérémy Clotault

► To cite this version:

Frédéric Proïa, Alix Pernet, Tatiana Thouroude, Gilles Michel, Jérémy Clotault. On the characterization of flowering curves using Gaussian mixture models. Journal of Theoretical Biology, 2016, 402, pp.75-88. 10.1016/j.jtbi.2016.04.022. hal-01388746

HAL Id: hal-01388746

https://hal.science/hal-01388746v1

Submitted on 12 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
ON THE CHARACTERIZATION OF FLOWERING CURVES
USING GAUSSIAN MIXTURE MODELS

FRÉDÉRIC PROÏA, ALIX PERNET, TATIANA THOUROUDE, GILLES MICHEL, AND JÉRÉMY CLOTAULT

ABSTRACT. In this paper, we develop a statistical methodology applied to the characterization of flowering curves using Gaussian mixture models. Our study relies on a set of rosebushes flowering data, and Gaussian mixture models are mainly used to quantify the reblooming properties of each one. In this regard, we also suggest our own selection criterion to take into account the lack of symmetry of most of the flowering curves. Three classes are created on the basis of a principal component analysis conducted on a set of reblooming indicators, and a subclassification is made using a longitudinal *k*-means algorithm which also highlights the role played by the precocity of the flowering. In this way, we obtain an overview of the correlations between the features we decided to retain on each curve. In particular, results suggest the lack of correlation between reblooming and flowering precocity. The pertinent indicators obtained in this study will be a first step towards the comprehension of the environmental and genetic control of these biological processes.

# 1. INTRODUCTION AND MOTIVATIONS

As it is explained by Putterill et al. [25], matching the flowering period with the best climatic conditions is a crucial step for wild plants to obtain a high fertility rate. In agriculture, the amount of seeds and fruits produced by plants is directly related to their ability to produce a great number of flowers, hence flowering is extremely important for high-yields crops. For ornamental plants, obtaining a large number of flowers over the longest period of the year is an important breeding objective.

Plants present a large diversity of flowering patterns between taxa and suitable parameters are necessary to summarize these flowering profiles. Flowering curves, counting the number of flowers observed for a plant at regular time intervals, can be obtained from field scorings. Statistical methods have been rarely used to efficiently describe and compare flowering curves. As an example, regression curves have been used to fit flowering curves, but only for *once-flowering* plants whose curve shape skewed from normality (see Clark and Thompson [6]). Especially in horticulture, annual flowering curves are sometimes much more complex. *Reblooming* – or *recurrent flowering* – plants are able to flower and fructify several times over the year. Such

<footnote>*Key words and phrases.* Flowering curves, Reblooming behavior, Recurrent flowering, Gaussian mixture models, Longitudinal *k*-means algorithm, Principal component analysis, Characterization of curves, Classification of curves.</footnote>
plants are found among several ornamental species, like irises, hydrangeas, daylilies or roses, but also in fruit-producing species like strawberry or raspberry plants.

For roses (*Rosa* sp. or genus *Rosa*), flowering traits are particularly important, either for cut or garden roses. In Occident, the nineteenth century represents a golden age for rosebush breeding. It involved the creation of many cultivars, with the introduction of new traits in created hybrids, as explained by Oghina-Pavie [23]. Very early in this century, breeding activities have aimed at obtaining earlier – or later – flowering cultivars to increase the range of flowering periods (Oghina-Pavie, pers. comm.). Later in this century, the reblooming trait became the most important trait in rosebush breeding. Current modern roses result from crosses between reblooming Chinese roses and once-flowering European roses, obtained during the nineteenth century, according to Wylie [31]. By the number of created cultivars and by the diversification of flowering profiles, rosebush genetic resources of the nineteenth century are probably the base of the most interesting methodologies to be developed in characterizing flowering curves.

<figure><img src="image_2.png" /><figcaption>FIGURE 1.1. Flowering stages for a rose (photo credit: Ballerie, 2012). On the left, [a] shows the detachment of one petal from the bud and [b] shows the detachment of one sepal from the bud. On the middle, a standard state of open flower with no withered petal is depicted. On the right, [c] points out the last non-withered petal. Were considered all flowers whose developmental stage lied between the left and right pictures.</figcaption></figure>

The biological sample analyzed in this article is composed of 329 exploitable flowering curves obtained in 2012 in the rose garden “Loubert” (Les Rosiers-sur-Loire, France). The studied genotypes were predominantly bred during the nineteenth century. For each genotype, the number of open flowers was counted almost each week between May, 10th and November, 15th. For the most widespread case, were considered as open flowers developmental stages between these two following stages: (1) flower bud with at least one sepal detached from petals and at least one petal detached from the others (except for simple flowers, having five petals) and (2) flower with at least one petal which remains with original aspect and colour (see Figure 1.1 above). The plant shape (sphere, cylinder or cone), circumference and height were measured for calculation of flower density (number of flowers per m²). The
mean number of flowers within an inflorescence was also counted. The dataset originally contained some irregularities and missing values, and different temporal lags between consecutive measures. All these issues have been carefully dealt with by the authors, but occasional presence of residual artificial values cannot be excluded.

In rosebush, the contemporary works of Iwata et al. [17] and Kawamura et al. [19] have highlighted the fact that flowering process is tightly linked to the branching process of the plant. In once-flowering cultivars, inflorescences are produced in the spring by the development of shoots from axillary buds of shoots from the previous year. Later in the year, new indeterminate shoots are produced and remain vegetative (having no flower). Inflorescence will develop the year after from axillary buds, borne by these vegetative shoots. In reblooming cultivars, either axillary buds will give inflorescence, or new determinate shoots terminated by an inflorescence will emerge successively from older shoots. Therefore, the best way to characterize rose flowering profile would be to differentiate the number of flowers produced by each shoot developed along the year. As an illustration of decomposition of the flowering shoot by shoot, Durand et al. [11] tried to model biennial bearing in apple trees. For a large sample of elderly rosebushes with many shoots like the one studied in this article, this represents a huge and laborious work. Statistical methods are therefore needed to characterize flowering profiles (flowering date, flowering intensity, reblooming magnitude) from flowering curves obtained by counting flowers along the year in the whole plant. As for the characterization of reblooming, it is especially challenging to distinguish a long unique flowering period from several partially overlapping ones, corresponding to successive floral initiations.

Mixture models have for a long time been popular in life sciences, especially in biology and genetics, in fact since the seminal works of Pearson in the late nineteenth century. We guide the reader to the far from exhaustive mixtures applications to biology by Hale and Knott [16] in 1992, Lynch and Walsh [20] in 1998, Detilleux and Leroy [10] in 2000, Boettcher et al. [3] in 2005, Choi et al. [5] in 2010, Shekofteh et al. [29] in 2015, and references therein. The ability of a Gaussian mixture to split an apparently chaotic whole phenomenon into simple components and to highlight hidden structures was our main motivation to apply such models to flowering curves. Indeed, a rosebush is made of stems among which one may start to bloom while another starts to lose its flowers. On a whole set of branches, the resulting phenomenon is not suitably explained from a deterministic approach consisting in counting all variations from one week to the other. On the contrary, the *waves* mechanism of Gaussian mixture models seems to form a relevant alternative, as we will see in the Appendix. The paper is organized as follows. Sections 2 and 3 are devoted to the statistical tools that we intend to customize and to their applications on our dataset, respectively for the characterization and the classification of the flowering curves. In particular, a theoretical background is supplied, when necessary. Some concluding remarks are given in Section 4 and a schematic example is provided in the Appendix, to justify our choice of Gaussian mixture models.
## 2. AN APPLICATION OF GMM TO FLOWERING CURVES

This section is devoted to the application of the Gaussian Mixture Model – short- ened from now on GMM – on a set of flowering curves according to a statistical methodology that we will detail. Firstly, we need to supply a short theoretical background about GMM (see [21]–[22] for more details).

### 2.1. The Gaussian Mixture Model.
Consider a set $(X_1, \dots, X_n)$ of $n$ real-valued random variables that we want to divide in $k$ classes. For all $1 \le i \le n$, we denote by $Z_i$ the latent random variable in $\{1, \dots, k\}$ corresponding to the class of $X_i$. We suppose that $(Z_1, \dots, Z_n)$ are independent and have the same distribution as a random variable $Z$ such that, for all $1 \le j \le k$,

$$ \mathbb{P}(Z = j) = \pi_j. $$

In addition, we suppose that for all $1 \le i \le n$, the random variable $X_i | \{Z_i = j\}$ has a $N(\mu_j, \sigma_j^2)$ distribution and accordingly, $\pi_j$ stands for the proportion of the class $j$ in the whole population. For all $x \in \mathbb{R}$, the distribution of the mixture is

$$ f_{\text{GM}}(x) = \sum_{j=1}^{k} \pi_j f(x | \mu_j, \sigma_j^2) $$

where $f(\cdot | \mu_j, \sigma_j^2)$ is the Gaussian distribution function with parameters $\mu_j$ and $\sigma_j^2$. If we consider that, given a subdivision in $k$ classes (that is, conditionally on the latent variables), the sequence $(X_1, \dots, X_n)$ is made of independent variables, then the (incomplete) log-likelihood for a set $\beta_k = (\pi_1, \dots, \pi_k, \mu_1, \dots, \mu_k, \sigma_1^2, \dots, \sigma_k^2)$ of $3k$ parameters is given for any observation $x = (x_1, \dots, x_n) \in \mathbb{R}^n$ by

$$ (2.1) \qquad \ln \ell_{\text{GM}}(x | \beta_k) = \sum_{i=1}^{n} \ln \left( \sum_{j=1}^{k} \pi_j f(x_i | \mu_j, \sigma_j^2) \right) = \sum_{i=1}^{n} \ln f_{\text{GM}}(x_i). $$

The classic approach (see e.g. [7]–[32]) to estimate the $3k-1$ parameters, considering the relation $\pi_1 + \dots + \pi_k = 1$, is to run the so-called Expectation-Maximisation algorithm [9] to maximize the above log-likelihood. The resulting estimator $\tilde{\beta}_k$ is finally used to classify the observations via the Bayes’ theorem. Namely,

$$ \tilde{\mathbb{P}}(Z_i = j | X_i = x_i) = \frac{\tilde{\pi}_j f(x_i | \tilde{\mu}_j, \tilde{\sigma}_j^2)}{\sum_{\ell=1}^{k} \tilde{\pi}_\ell f(x_i | \tilde{\mu}_\ell, \tilde{\sigma}_\ell^2)} $$

usually leading to the posterior classification rule given, for all $1 \le i \le n$, by

$$ \tilde{Z}_i = \arg \max_{1 \le j \le k} \tilde{\mathbb{P}}(Z_i = j | X_i = x_i). $$

For a given number of classes $k$, the Bayesian information criterion [27] is

$$ (2.2) \qquad \text{BIC}(k) = -2 \ln \ell_{\text{GM}}(x | \tilde{\beta}_k) + (3k - 1) \ln n $$

where the log-likelihood is given in (2.1). It is then a natural solution to select

$$ \tilde{k} = \arg \min_{1 \le k \le k_m} \text{BIC}(k) $$

for an arbitrary upper bound $k_m$.
<figure><img src="image_3.png" /><figcaption>FIGURE 2.1. Example of GMM on simulated data according to $k = 2$, $\pi = (0.3, 0.7)$, $\mu = (2, 5)$, $\sigma^2 = (2, 0.5)$. The coloured curves are the Gaussian components whereas the black one is the resulting mixture. The evolution of BIC is shown alongside for $k_m = 8$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 2.2. Example of GMM on simulated data according to $k = 3$, $\pi = (0.5, 0.2, 0.3)$, $\mu = (1, 3, 8)$, $\sigma^2 = (0.3, 1, 0.8)$. The coloured curves are the Gaussian components whereas the black one is the resulting mixture. The evolution of BIC is shown alongside for $k_m = 8$.</figcaption></figure>
In the following, the R library mclust [13]–[12] is used to run the GMM procedures. On Figures 2.1–2.2 above, two examples of GMM are represented on simulated data with n = 1000. The coloured curves stand for the Gaussian components whereas the black one is the resulting mixture. The evolution of BIC is shown alongside for k_m = 8, leading to $\tilde{k} = 2$ and $\tilde{k} = 3$, respectively. A learning set is available on which we assume that the phenotypic behavior in terms of blooming is well-known. Working on these learning curves will help us to evaluate some parameters, such as α in the penalization that we will introduce. To apply a GMM on a flowering curve, we suggest to handle the curve as a probability distribution function and to simulate a sample in accordance with. In this context, the fitted mixture model characterizes the temporal probability distribution of the amount of flowers on the rosebush, along the year.

## 2.2. The flowering curve as a distribution function.
First, we need to precise that what we call a flowering curve is a discontinuous set of measures. The continuous line between all points that is represented on the figures throughout the study is only a visual tool, it is never technically used. In particular, the lack of information between each measure was precisely our main motivation to use the step interpolation function that we are going to describe. Let us consider the observed path $(y_{i,1}, \dots, y_{i,d}) \in \mathbb{R}^d$ associated with the i-th curve of the dataset, and $(t_{i,1}, \dots, t_{i,d}) \in \mathbb{N}^d$ corresponding to the instants of measure. For all $1 \le \ell \le d$, we build the step function $g_i$ according to

$$ (2.3) \qquad g_i(x) = y_{i,\ell} \quad \text{for } x \in \left[ t_{i,\ell} - \frac{t_{i,\ell} - t_{i,\ell-1}}{2} ; t_{i,\ell} + \frac{t_{i,\ell+1} - t_{i,\ell}}{2} \right] $$

with the convention that $t_{i,0} = t_{i,1} - 1$ and that $t_{i,d+1} = t_{i,d} + 1$. On Figures 2.3 and 2.4 below, some examples are provided from the dataset. The blue lines are the flowering curves and the magenta rectangulations stand for the associated step functions $g_i$. Hence, the step function $\tilde{g}_i$ defined for all $t_{i,0} \le x \le t_{i,d+1}$ by

$$ \tilde{g}_i(x) = \frac{g_i(x)}{A_i} \quad \text{where} \quad A_i = \int_{t_{i,0}}^{t_{i,d+1}} g_i(x) \, dx = \frac{1}{2} \sum_{\ell=1}^{d} y_{i,\ell} (t_{i,\ell+1} - t_{i,\ell-1}), $$

can be seen as a probability distribution function. A GMM applied on a sample $(\tilde{Y}_{i,1}, \dots, \tilde{Y}_{i,n_s})$ of $n_s$ independent random variables distributed according to $\tilde{g}_i$ provides a tool to characterize the recurrent blooming of the i-th rosebush of the dataset. The smoothing effect on the chaotic behavior of most of the flowering curves is a substantial improvement compared to the deterministic strategy consisting in counting all changes of variations: the induced waves mechanism of GMM seems somehow more adapted to the environmental and genetic reality of the plant, as it is shown in the schematic example provided in the Appendix. Through this illustration, we intend to highlight that the biological process has to be explained by the superimposition of hidden waves of flowering, and not simply by the evaluation of changes in the number of flowers over time. A weather effect can also involve a locally chaotic behavior, for example rain and wind are likely to produce a sudden fall of petals leading to biased measurements. Indeed, the counting process does not
distinguish the ability of the rosebush to produce flowers from their resistance to time and weather.

<figure><img src="image_5.png" /><figcaption>FIGURE 2.3. Examples of rectangulations (in magenta) of two flowering curves (in blue) using step functions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 2.4. Examples of rectangulations (in magenta) of two flowering curves (in blue) using step functions.</figcaption></figure>

**2.3. An adapted criterion.** The main issue arising in the GMM procedure is that, due to environment and biological effects, most often the generated sample fails in the usual statistical testing procedures for Gaussianity (such as Shapiro's) on the learning curves having a unique peak of blooming (clearly perceptible on the curves represented throughout the paper). This phenomenon is also observed in [6] and references within, though to a far lesser extent. The resulting asymmetry may lead to inappropriate conclusions from the GMM algorithm, in particular when some clusters become very close to each other to take into account all irregularities
of the curve. To reduce this phenomenon, we suggest to use a selection criterion penalizing the $3k - 1$ parameters to estimate (via the ordinary BIC), but also the smallest difference between the estimated means of the $k$ Gaussian components. Using the same notations as in the definition of BIC in (2.2) given an estimation on $k$ clusters, let us define

$$ (2.4) \qquad \text{BIC}^*(k) = (c + \text{BIC}(k)) (1 + e^{-\alpha d_k}) $$

where $\alpha, c \in \mathbb{R}^+$, $d_1 = +\infty$ and for $k \ge 2$,

$$ (2.5) \qquad d_k = \min_{\substack{1 \le j_1, j_2 \le k \\ (j_1 \neq j_2)}} |\tilde{\mu}_{j_1} - \tilde{\mu}_{j_2}|. $$

Figure 2.5 shows the evolution of the penalization coefficient according to $d_k$, for $\alpha \in \{0.5, 1, 1.5, 2, 2.5\}$. As we can see, this coefficient aims to sharply penalize any model where the smallest difference between the estimated means becomes less than 1.5, on the whole. According to us, the probability is higher that such a situation corresponds to a lack of symmetry of the current flowering.

<figure><img src="image_7.png" /><figcaption>FIGURE 2.5. Evolution of the penalization coefficient $(1 + e^{-\alpha d_k})$ according to $d_k$, for different values of $\alpha$.</figcaption></figure>

As for the usual BIC, our choice relies on

$$ \tilde{k}^* = \arg \min_{1 \le k \le k_m} \text{BIC}^*(k) $$

for an arbitrary upper bound $k_m$. The estimation of $c$ is not difficult, it only has to ensure that $\tilde{c}^* + \text{BIC}(k) > 0$ and one can choose for example

$$ \tilde{c}^* = \max_{1 \le k \le k_m} |\text{BIC}(k)| - \min_{1 \le k \le k_m} \text{BIC}(k). $$

To evaluate $\tilde{\alpha}^*$, we make $\alpha$ vary on a grid and experiments are conducted on the learning set. Note that $c$ and $\alpha$ depends on $i$, meaning that each curve has its own
parameters to be estimated in the criterion. From a practical point of view as we will see in the next section, $\tilde{\alpha}^*$ only changes from a class of curves to another. On Figures 2.6–2.9 are represented four examples of flowering curves having different blooming behaviors supplied with their associated $g_i$ step functions. On the right, the evolution of BIC and $BIC^*$ for $k_m = 8$ is also given. As one can see on the graphs, $BIC^*$ plays a moderation role and suggests to select $\tilde{k}^* = 1, 2, 3$ and $3$ respectively, whereas BIC suggests to choose the quite unrealistic values $\tilde{k} = 3, 8, 8$ and $8$ respectively, for the aforementioned reasons. In fact as the reader will observe, for some curves (as the one of Figure 2.7), BIC suggests to select the maximal number of components whereas the common sense would have been to choose the same value of $k$ using BIC or $BIC^*$. However as we deal with hundreds of curves in each dataset, it is essential for us that an algorithmic procedure enables to select $k$, with no *human* intervention.

<figure><img src="image_8.png" /><figcaption>FIGURE 2.6. Example of flowering curve with its rectangulation and evolution of the associated BIC (in black) and $BIC^*$ (in blue) from GMM estimation.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIGURE 2.7. Example of flowering curve with its rectangulation and evolution of the associated BIC (in black) and BIC* (in blue) from GMM estimation.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 2.8. Example of flowering curve with its rectangulation and evolution of the associated BIC (in black) and BIC* (in blue) from GMM estimation.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 2.9. Example of flowering curve with its rectangulation and evolution of the associated BIC (in black) and BIC* (in blue) from GMM estimation.</figcaption></figure>

The estimation of $c$ was automatically conducted and values $\tilde{\alpha}^* = 2, 2.5, 2.5$ and $2.5$ were chosen for $\alpha$. The results of the GMM algorithm on these curves for $\tilde{k}^*$ clusters are given on Figure 2.10–2.13 together with the histogram of the generated samples $(\tilde{y}_{i,1}, \dots, \tilde{y}_{i,n_s})$ with $n_s = 1000$.

<figure><img src="image_12.png" /><figcaption>FIGURE 2.10. Example of the GMM algorithm running on a flowering curve, for $\tilde{k}^* = 1$ selected cluster. On the right, the coloured curve is the estimated Gaussian component and the histogram of the generated sample is superimposed.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 2.11. Example of the GMM algorithm running on a flowering curve, for $\tilde{k}^* = 2$ selected clusters. On the right, the coloured curves are the estimated Gaussian components and the histogram of the generated sample is superimposed.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 2.12. Example of the GMM algorithm running on a flowering curve, for $\tilde{k}^* = 3$ selected clusters. On the right, the coloured curves are the estimated Gaussian components and the histogram of the generated sample is superimposed.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIGURE 2.13. Example of the GMM algorithm running on a flowering curve, for $\tilde{k}^* = 3$ selected clusters. On the right, the coloured curves are the estimated Gaussian components and the histogram of the generated sample is superimposed.</figcaption></figure>

### 3. INDICATORS FOR PCA AND CLASSIFICATION

The dataset of flowering curves is very heterogeneous, this is the main motivation to find a set of indicators allowing to classify each rosebush. As for GMM, we also need to shortly summarize the main principles of the k-means for longitudinal data algorithm [15] – shortened from now on KML. The objective is first to subdivide the dataset into classes according to the reblooming behavior of each curve, and then to provide a representative mean curve for each cluster of a subclassification.

3.1. The k-means for longitudinal data. We consider a set of $n$ random vectors of $\mathbb{R}^d$ defined as $Y_i = (Y_{i,1}, \dots, Y_{i,d})$, for all $1 \le i \le n$. The KML algorithm merely works like an usual k-means algorithm on a set of paths $(y_1, \dots, y_n)$ that may be thought as $n$ time-related trajectories of length $d$. After convergence of the classification algorithm in $k$ clusters, denote by

$$B_k = \frac{1}{n} \sum_{j=1}^{k} \tilde{n}_j (\bar{y}_j - \bar{y})(\bar{y}_j - \bar{y})'$$

where $\tilde{n}_j$ is the size of cluster $j$, $\bar{y}_j$ the corresponding mean trajectory and $\bar{y}$ the mean trajectory of the whole set. Define also

$$W_k = \frac{1}{n} \sum_{j=1}^{k} \sum_{i=1}^{\tilde{n}_j} (y_i^j - \bar{y}_j)(y_i^j - \bar{y}_j)$$

where $y_i^j$ is the $i$-th path in cluster $j$. The between-variance matrix $B_k$ can be seen as an estimator of $\mathbb{V}(\mathbb{E}[Y|Z])$, where $Y$ stands for a representative random vector of an independent population and $Z$ is the latent classification random variable
associated with Y, having k modalities. Similarly, the within-variance matrix $W_k$ is an estimator of $E[V(Y|Z)]$. The number of clusters is selected by maximization of the Calinski-Harabasz criterion [4] given by

$$CH(k) = \frac{(n-k) \text{tr}(B_k)}{(k-1) \text{tr}(W_k)}$$

where $B_k$ and $W_k$ are estimated for k clusters. The CH selection is then

$$\tilde{k} = \arg \max_{2 \le k \le k_m} CH(k)$$

for an arbitrary upper bound $k_m$. The centroids of the $\tilde{k}$ clusters will be seen as the representative curves of each class. To run KML algorithm, we will use the R library kml [14]. The opportunity to deal with temporal missing values was our main motivation to make use of KML instead of standard k-means for random vectors. Indeed as mentioned in the introduction, a non-negligible amount of data is missing in our curves, and a consistent completion was essential. In Figures 3.1-3.2 below, an example of KML is represented on simulated data with $n = 35$ ($n_1 = 20$, $n_2 = 15$), $d = 20$ and two different patterns with an additive noise, constrained to stay nonnegative. On the first figure, the patterns (chosen to look like flowering curves) are given, the generated curves together with their centroids found by the algorithm appear alongside. On the second figure, the evolution of the size of each cluster and the evolution of the CH criterion associated are given, for $k_m = 6$. According to the CH criterion, $\tilde{k} = 2$ is obviously selected.

<figure><img src="image_16.png" /><figcaption>FIGURE 3.1. Example of KML running on a set of $n_1 = 20$ and $n_2 = 15$ simulated trajectories according to the patterns in the left, for $k = 2$ and an additive noise. The centroids are highlighted.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>FIGURE 3.2. Evolution of the clusters size on the barplots and evolution of the CH criterion related to the above example, for $k_m = 6$.</figcaption></figure>

3.2. A set of indicators and a reblooming classification. Let us denote by $t_{i,R}$ the first instant of significative reblooming of the i-th rosebush, that is the first time where the density exceeds a given threshold after the end of its first flowering. This threshold is calculated as a percentage (25% in our examples) of the peak value during the first significative flowering, itself detected using a similar algorithm. Consider the areas

$$ (3.1) \qquad P_i = \int_{t_{i,0}}^{t_{i,R}} g_i(x) \, dx \quad \text{and} \quad R_i = \int_{t_{i,R}}^{t_{i,d+1}} g_i(x) \, dx $$

with the convention that $t_{i,R} = t_{i,d+1}$ if the i-th rosebush does not have a significative reblooming (we recall that d is the observation vector size), and where $g_i$ is the step function given in (2.3). $P_i$ stands for the area under the first significative flowering and $R_i$ for the area under the whole reblooming period, possibly zero. On Figures 3.3 and 3.4, we present four examples of such detection, leading to $t_{i,R} = 13$ (top left), 15 (top right), 33 (bottom left) and 14 (bottom right). The areas $P_i$ and $R_i$ are tinted, respectively in blue and red (note that, strictly speaking, it is only an approximation of the areas which is coloured, we have not represented the rectangulations on these graphs, which are the effective tools used to compute the areas).
<figure><img src="image_18.png" /><figcaption>FIGURE 3.3. Examples of separation of the main peak area and of the reblooming area for two flowering curves. The first significative flowering is tinted in blue and the reblooming area detected by our algorithm is tinted in red.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIGURE 3.4. Examples of separation of the main peak area and of the reblooming area for two flowering curves. The first significative flowering is tinted in blue and the reblooming area detected by our algorithm is tinted in red.</figcaption></figure>

Now we consider the reblooming magnitude that we define as

$$ (3.2) \qquad M_i = \frac{R_i}{P_i + R_i}. $$

This indicator is going to play a primordial role in our classification. Literally, we have built the OF/R₁/R₂ classes by using a k-means algorithm on the first
PCA plane (see Figure 3.6 later) according to the indicators that we are going to
describe. OF stands for *once-flowering*, R₁ and R₂ for *weak* and *strong reblooming*.
The numeric indicators that we have decided to retain to form a basis of comparison
between the whole curves are given below. They come either from the plant features
or from the statistical analysis. We use the notations $\tilde{\pi}_{max} = \max(\tilde{\pi}_1, \dots, \tilde{\pi}_{\tilde{k}^*})$,
$\tilde{\mu}_{min} = \min(\tilde{\mu}_1, \dots, \tilde{\mu}_{\tilde{k}^*})$ and $\tilde{\sigma}_{min}^2 = \min(\tilde{\sigma}_1^2, \dots, \tilde{\sigma}_{\tilde{k}^*}^2)$, corresponding to the GMM
estimations. Among all sets of indicators on which a PCA has been conducted, the
following one has given the most meaningful results.

* **Nb.Clust**: the number $\tilde{k}^*$ of clusters selected by our BIC* criterion in the GMM applied on the curve.

* **LMax.P**: $\ln(\tilde{\pi}_{max})$, the log-weight affected to the main flowering.

* **Rat.P**: $(1 - \tilde{\pi}_{max})/\tilde{\pi}_{max}$, the ratio between the weight affected to the reblooming area and the weight of the main flowering.

* **Min.M**: $\tilde{\mu}_{min}$, the lowest estimated mean.

* **DifMax.M**: $d_{\tilde{k}^*}$ as it is defined in (2.5) for $\tilde{k}^* \ge 2$, the highest difference between the estimated means. By convention, difference between the last and first instants of measure for $\tilde{k}^* = 1$.

* **Min.V**: $\tilde{\sigma}_{min}^2$, the lowest estimated variance.

* **Area.Peak.Inflo**: the area under the main flowering, defined as $P_i$ in (3.1), divided by the mean number of flowers within an inflorescence.

* **Reb.Mag**: the reblooming magnitude defined as $M_i$ in (3.2).

* **Max.Peak**: the highest value of the first significative flowering.

* **Max.Reb**: the highest value of the reblooming period.

* **First.Reb**: the first instant of the reblooming period called $t_{i,R}$ in (3.1), last instant of measure if there is no significative reblooming.

* **Rat.Peak**: the ratio between **Max.Reb** and **Max.Peak**.

* **First.Flo**: the first instant of flowering.

* **Cont.Flo**: an indicator in {0, 1} characteristic of plants having a continuous significative flowering. In the framework of this study, a continuous flowering is a feature of a curve taking most of the time (> 80%) positive values.
<figure><img src="image_20.png" /><figcaption>FIGURE 3.5. Projection of the indicators on the main factorial plane of the PCA on the left, and barplot of the percentage of inertia associated with each axis on the right.</figcaption></figure>

The PCA analysis reveals a dominant eigenvalue expressed on the first axis and another eigenvalue far less significative, expressed on the second axis. In the context of the study, we will content ourselves with this couple of eigenvalues. On Figure 3.5, all our indicators have been projected on the main factorial plane and we make the following observations: the reblooming and the precocity seem to form the main orthogonal features of a rosebush, in order of importance. To be complete, let us precise that more indicators have been initially tested (LMin.P, Max.M, Max.V, Moy.M and Moy.V (mean values), DifMax.M, Inflo (intensity of inflorescence), Area.Peak, etc.) With equivalent results (more than 0.95 of empirical correlation between Area.Peak and Max.Peak, for example), the smallest set has been preferred, for purposes of parsimony. One can already notice that a few curves are considered as OF whereas their reblooming magnitudes are non-zero: according to us, this is due either to artefacts for biological or environmental reasons, or imprecisions in gathering data.

3.2.1. *The reblooming indicators.* Six indicators are strongly correlated to describe the reblooming behavior of a rosebush: Nb.Clust, Reb.Mag, Rat.P, Rat.Peak, Cont.Flo and LMax.P, the reasons being obvious for most of them. According to our model, the reblooming deepens when $\tilde{\pi}_{max}$ decreases, this explains the fact that Rat.P is positively and LMax.P negatively correlated with the reblooming behavior. First.Reb and DifMax.M play a negative role: DifMax.M is affected by the large values of once-flowering curves whereas First.Reb increases when reblooming starts later. The positive correlation of Max.Reb with the reblooming phenomenon seems also quite obvious.
3.2.2. *The precocity indicators.* Four indicators are related to the precocity of the rosebush: First.Flo, Min.M, Area.Peak.Inflo and Max.Peak. Some of them have a trivial explanation, whereas it is quite unexpected for us to observe that the main flowering seems positively related to the precocity: a rosebush having a later first flowering produces on average a more abundant first flowering. A more favourable weather might be a trail to explain it.

Finally, the last indicator Min.V does not appear to play any role in the reblooming or the precocity of the rosebush. To sum up, the reblooming features are clearly described by the classification OF/R₁/R₂. In the following section, we aim to look for a subclassification of the curves based on the precocity indicators. As we will see, this is unfortunately far less convincing. There is an effect of the second axis on the subclasses but it is not as perceptible as the effect of the first axis on the classes.

<figure><img src="image_21.png" /><figcaption>FIGURE 3.6. Projection of the individuals on the main factorial plane of the PCA. From a k-means classification, OF are coloured in blue, R₁ in red and R₂ in green. The checked points correspond to the flowering curves of Figure 3.7 and the centroids are indicated.</figcaption></figure>

The fuzzy boundary between green (R₂) and red (R₁) individuals is not problem-
atic in our framework since R₁ and R₂ are only separated on a *quantitative* basis.
The presence of blue (OF) individuals in the red (R₁) area, however separated on a
*qualitative* basis, seems somehow more annoying. They correspond to the aforemen-
tioned curves having a non-zero reblooming but weak enough to be considered as
artefact. For R₂, the lack of concentration of the individuals in the factorial plane is
justified by the strong heterogeneity of the reblooming curves, despite the numerous descriptive indicators. By way of example, we have represented on Figure 3.7 the flowering curves checked on Figure 3.6.

<figure><img src="image_22.png" /><figcaption>FIGURE 3.7. Examples of the peak/reblooming areas associated with the flowering curves checked in Figure 3.6. The first one has a massive flowering despite a low reblooming magnitude whereas the second one has the highest reblooming magnitude of the dataset.</figcaption></figure>

Both of them are $R_2$ curves. The first one has a quite low reblooming magnitude (≈ 0.14) but the flowering is so abundant that the peak indicators, though weakly correlated with the first axis, project the rosebush to the bottom-right corner of the factorial plane. The second one has the highest reblooming magnitude of the dataset (≈ 0.86) and is accordingly also located at the right of the plane. These examples show the huge heterogeneity of the reblooming classes. This prior classification leads to consider a substantial part of the dataset as once-flowering rosebushes whereas reblooming classes are dominated by weakly reblooming plants. The sizes are the following: on $n = 329$ exploitable flowering curves, $n_0 = 152$ are OF, $n_1 = 127$ are $R_1$ and $n_2 = 50$ are $R_2$ (see Figure 3.11 below for a recap chart of the subclasses).

3.3. **A subclassification of the curves using KML.** The next step is to run the KML algorithm in each class of curves (OF/$R_1/R_2$), to highlight similar behaviors. We start by standardizing the dataset so as to restrict any scale effect. For all $1 \le l \le d$, let

$$z_{i,l} = \frac{y_{i,l}}{\|y_i\|_{\infty}}$$

where $y_i = (y_{i,1}, \dots, y_{i,d})$ is the observed path associated with the $i$-th curve of the dataset, and $\|\cdot\|_{\infty}$ is the usual infinity norm of $\mathbb{R}^d$. On the basis of the CH criterion (see Section 3.1) and looking for reasonable sizes in each subclass, $\tilde{k} = 4, 3$ and $2$ clusters are selected for OF, $R_1$ and $R_2$ respectively. The representative curves of each subclass are given on Figure 3.8–3.10.
<figure><img src="image_23.png" /><figcaption>FIGURE 3.8. Representative curves of the $\tilde{k} = 4$ subclasses of the OF standardized dataset, from the KML algorithm.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 3.9. Representative curves of the $\tilde{k} = 3$ subclasses of the $R_1$ standardized dataset, from the KML algorithm.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 3.10. Representative curves of the $\tilde{k} = 2$ subclasses of the $R_2$ standardized dataset, from the KML algorithm.</figcaption></figure>
The proportions of the whole classification are represented thereafter in Figure 3.11. Among the $n_0 = 152$ curves in the OF class, $n_{01} = 65$ are $C_1$ (blue centroid), $n_{02} = 49$ are $C_2$ (red centroid), $n_{03} = 24$ are $C_3$ (green centroid) and $n_{04} = 14$ are $C_4$ (magenta centroid). Similarly, among the $n_1 = 127$ curves in the $R_1$ class, $n_{11} = 71$ are $C_1$ (blue centroid), $n_{12} = 29$ are $C_2$ (red centroid) and $n_{13} = 27$ are $C_3$ (green centroid). Finally, among the $n_2 = 50$ curves in the $R_2$ class, $n_{21} = 34$ are $C_1$ (blue centroid) and $n_{22} = 16$ are $C_2$ (red centroid).

<figure><img src="image_26.png" /><figcaption>FIGURE 3.11. Descriptive statistics: pie chart associated with the proportions of the reblooming classes OF/$R_1/R_2$ in the dataset and with their KML subclasses.</figcaption></figure>

The main indicator of selection is clearly the reblooming magnitude, but the representative curves highlight the orthogonal indicator lying in the precocity of the first flowering peak. It is especially perceptible on the once-flowering (OF) and weakly reblooming ($R_1$) curves, whereas it seems to vanish on the highly reblooming ($R_2$) curves where the chaotic behavior yields difficulties to clearly identify the beginning of the flowering.

4. CONCLUDING REMARKS

According to the authors, this work is a further evidence of the well-established consistence of Gaussian mixture models with biological data. The ability of a Gaussian mixture to smooth the observations and to highlight hidden structures form a convincing answer to the main features of flowering curves. In addition, the *waves* mechanism is not only a suitable statistical modelling, it also has a biological and genetic interpretation. The modified BIC* criterion we have suggested to use may be seen as an artificial correction of the potentially asymmetry of the flowering phenomenon, and somehow interpreted as a simplistic solution. Hence, non-Gaussian mixture models could be used for further investigations. The authors do not exclude the possibility that a peak of flowering may be better explained by an almost Gaussian but slightly asymmetric distribution. A deep statistical investigation of the symmetry in the OF dataset should lead to sufficient evidence, nevertheless the parameters estimation will pose some significant issues.

By comparison, the subclassification using the *k*-means longitudinal algorithm is far less meaningful. Not surprisingly, the precocity of the rosebushes is the main
highlighted feature and, as we have seen in this work, a reasonable interpretation is only possible in the OF class. The authors have also considered the eventuality of a deformation model in which, for all $1 \le \ell \le d$,

$$Y_{i,\ell} = a_i R_{i,\ell-t_i} + \varepsilon_{i,\ell}$$

where $(Y_{i,1}, \dots, Y_{i,\ell})$ is the random flowering curve of the $i$-th rosebush, $a_i$ and $t_i$ are real parameters, $(R_{i,1}, \dots, R_{i,\ell})$ is the representative curve of the associated KML subclass and $(\varepsilon_{i,1}, \dots, \varepsilon_{i,\ell})$ is a white noise sequence. Thus, each flowering curve was seen as a linear deformation by centering/distension of its representative KML curve. Owing to the heterogeneity of the $R_1/R_2$ datasets, results obtained from a standard least squares approach were not satisfactory for our purposes. The estimations of $a_i$ and $t_i$ could have been relevant alternatives to Max.Peak and First.Flo in our PCA, but probably they would have led to the same kind of conclusions. Nevertheless, it seems that the two main orthogonal features characterizing a flowering curve are the reblooming phenomenon (number of clusters, magnitude, etc.) and the precocity (first instant of flowering, lowest estimated mean, etc.).

The fact that reblooming and precocity are found independent in this sample may have a genetical reason. Recurrent blooming is mainly controlled by a recessive locus that was recently identified as RoKSN by Iwata et al. [17]. The date of the first flowering has a more complex genetic determinism with different loci controlling this trait, as shown by Kawamura et al. [18] and Roman et al. [26]. The QTL having the main effect was proposed to correspond to the gene RoFT (see Otagaki et al. [24]). The work of Spiller et al. [30] has established that RoKSN and RoFT map respectively to linkage groups 3 and 4. Due to the genetic independence between these genes, specific alleles for these genes have not been conjointly inherited, and therefore these phenotypic traits are not correlated. Previously, Kawamura et al. [18] and Roman et al. [26] have shown a weak correlation between recurrent blooming and the precocity by studying $F_1$ progenies. In these progenies, reblooming rosebushes flower earlier than once-flowering ones. This difference can be explained by a genetic linkage between the recurrent blooming locus and a precocity locus (genes involved in gibberellic acid are in the vicinity and are good candidates). In the current study, such association is not found because the genetic basis is larger (more than 300 individuals) than in $F_1$ progeny (two parents) and linkage disequilibrium between two linked genes is more likely to be reduced by the number of meioses during rose history.

Environmental and genetic effects are not entirely explicable and separable due to the lack of repetitions of the measures. Indeed, a weather effect can not be neglected as flowering is highly controlled by environmental factors (see Andrès and Coupland [1] for a review). Likewise, in our probabilistic model all rosebushes are considered as independent. This hypothesis is questionable: we have every biological reasons to think that two plants growing in the same environmental conditions show some similitudes. Some rosebushes are also genetically close, which can lead to artificial associations between flowering parameters. Another direction for a future study lies in the consistence and the persistence of our conclusions on a temporal and
spatial prolongation, and in the comparative flowering behavior of genetically close rosebushes. Such an extended study could also provide answers to some questions of great biological interest, like the relation between the intensity of the first flowering of a rosebush and its reblooming nature.

Since the work of Semeniuk [28] in 1971 in *Rosa wichurana* and recent molecular characterization by Iwata *et al.* [17] in 2012, reblooming is considered to be controlled by a recessive allele of a single gene (*RoKSN*). Therefore, reblooming has until now been treated as a qualitative trait with as modalities ‘recurrent’ and ‘non-recurrent’ blooming roses (see de Vries and Dubois [8]) or sometimes ‘once-flowering’, ‘occasionally reblooming’ and ‘continuous flowering’ roses (see Iwata *et al.* [17]). The complexity of flowering curves and the analysis displayed in Figure 3.6 show that the reblooming trait should be treated as a semi-quantitative trait (once-flowering roses and roses with various reblooming intensities). The large sample size and its large diversity by comparison with previous studies probably explain the detection of the heterogeneity in flowering curves. Data produced and statistical tools developed in this work pave the way for the investigation of the environmental and genetic causes of this heterogeneity. This work is a preliminary step for finding new genes or new alleles controlling lower reblooming intensity than the one conferred by the Chinese copia allele of *RoKSN*. As an example, *Rosa damascena* ‘Four Seasons’, with occasional reblooming, was cultivated in the Middle East then in Europe long before the introduction of the first Chinese roses at the end of eighteenth century. Genetic methods like association mapping (see Balding [2] for a review), aiming at finding correlations between genetic markers and phenotypic traits in germplasm, may allow to complete the knowledge of both reblooming genetic control and rose breeding history.

**Acknowledgments.** The acquisition of flowering curves was supported by the "Région Pays de la Loire" in the framework of the FLORHIGE project, and by the department BAP of INRA in the framework of the SIFLOR project. We thank Thérèse and Raymond Loubert for providing access to their rose garden. We also thank Rachid Boumaza for his advices, Laure Ballerie for her work, Fabrice Foucher for his critical and pertinent view along the project, and the courageous people spending days and days making flowering measurements. Finally, we thank the Associate Editor and the two anonymous Reviewers for their suggestions and constructive comments which helped to improve substantially the paper.
APPENDIX. THE WAVES MECHANISM OF FLOWERING

<figure><img src="image_27.png" /><figcaption>FIGURE A.1. Example of six flowering stages (S1-S6) of a reblooming rosebush along the year (credit: free clipart) in which the evolution of four branches are represented (B₁, B₂, B₃, B₄).</figcaption></figure>

On Figure A.1 above, a simulation of the flowering stages of a reblooming rosebush along the year is given. Through this example, we aim to highlight the fact that merely counting flowers over time on a plant is not always relevant to study its flowering process. The six stages are described thereafter.
* S1 : before its first instant of flowering, the rosebush has developed two branches ($B_1$ and $B_2$).

* S2 : branch $B_2$ gives an abundant flowering.

* S3 : while $B_2$ is still blooming, branch $B_1$ produces a scattered flowering and a small branch $B_3$ is appearing at the root of $B_1$.

* S4 : branch $B_3$ immediately comes into flower whereas $B_2$ is withering, and a sizable branch $B_4$ is growing from the base of the rosebush.

* S5 : branches $B_1$ and $B_2$ are pruned after their whole decline while $B_4$ displays a sudden flowering.

* S6 : at the end of the flowering process, the rosebush is pruned.

The gap between the flowering of $B_1$ and $B_2$ explains the spread of the first flowering event, common to both once-flowering and reblooming roses. Flowering of $B_3$ and $B_4$ corresponds to reblooming events, only observable in reblooming roses. On Figure A.2, we have represented the shape of this simulated flowering curve, to show that a naive evaluation of the number of flowers on the plant along the year clearly leads to conclude that this rosebush is once-flowering. The superimposition of the reblooming waves detected by GMM allows the opposite conclusion which is biologically more relevant.

<figure><img src="image_28.png" /><figcaption>FIGURE A.2. Shape of the simulated flowering curve of the example above (black dotted line), the waves detected by the GMM algorithm are coloured.</figcaption></figure>
REFERENCES

[1] ANDRÈS, F., AND COUPLAND, G. The genetic basis of flowering responses to seasonal cues. *Nat. Rev. Genet.* **13-9** (2012), 627–639.

[2] BALDING, D. J. A tutorial on statistical methods for population association studies. *Nat. Rev. Genet.* **7-10** (2006), 781–791.

[3] BOETTCHER, P. J., MORONI, P., PISONI, G., AND GIANOLA, D. Application of a finite mixture model to somatic cell scores of italian goats. *J. Dairy. Sci.* **88-6** (2005), 2209–2216.

[4] CALINSKI, T., AND HARABASZ, J. A dendrite method for cluster analysis. *Commun. Stat. 3-1* (1974), 1–27.

[5] CHOI, H., QIN, Z. S., AND GHOSH, D. A double-layered mixture model for the joint analysis of DNA copy number and gene expression data. *J. Comput. Biol.* **17-2** (2010), 121–137.

[6] CLARK, R. M., AND THOMPSON, R. Estimation and comparison of flowering curves. *Plant. Ecol. Divers.* **4-2-3** (2011), 189–200.

[7] DAY, N. E. Estimating the components of a mixture of normal distributions. *Biometrika.* **56-3** (1969), 463–474.

[8] DE VRIES, D. P., AND DUBOIS, L. A. M. Inheritance of the Recurrent Flowering and Moss Characters in F₁ and F₂ Hybrid Tea × *R. centifolia muscosa* (Aiton) Seringe Populations. *Gartenbauwissenschaft.* **49-3** (1984), 97–100.

[9] DEMPSTER, A. P., LAIRD, N. M., AND RUBIN, D. B. Maximum likelihood from incomplete data via the EM algorithm. *J. Roy. Stat. Soc. B.* **39-1** (1977), 1–38.

[10] DETILLEUX, J., AND LEROY, P. L. Application of a mixed normal mixture model to the estimation of mastitis-related parameters. *J. Dairy. Sci.* **83** (2000), 2341–2349.

[11] DURAND, J., GUITTON, B., PEYHARDI, J., HOLTZ, Y., GUÉDON, Y., TROTTIER, C., AND COSTES, E. New insights for estimating the genetic value of segregating apple progenies for irregular bearing during the first years of tree production. *J. Exp. Bot.* **64-16** (2013), 5099–5113.

[12] FRALEY, C., AND RAFTERY, A. E. Model-based clustering, discriminant analysis and density estimation. *J. Am. Stat. Assoc.* **97** (2002), 611–631.

[13] FRALEY, C., RAFTERY, A. E., MURPHY, T. B., AND SCRUCCA, L. *mclust Version 4 for R: Normal Mixture Modeling for Model-Based Clustering, Classification, and Density Estimation*, 2012.

[14] GENOLINI, C., ALACOQUE, X., SENTENAC, M., AND ARNAUD, C. kml and kml3d: R packages to cluster longitudinal data. *J. Stat. Softw.* **65-4** (2015), 1–34.

[15] GENOLINI, C., AND FALISSARD, B. Kml: k-means for longitudinal data. *Comput. Stat.* **25** (2010), 317–328.

[16] HALEY, C. S., AND KNOTT, S. A. A simple regression method for mapping quantitative trait loci in line crosses using flanking markers. *Heredity.* **69-4** (1992), 315–324.

[17] IWATA, H., GASTON, A., REMAY, A., THOUROUDE, T., JEAUFFRE, J., KAWAMURA, K., HIBRAND-SAINT OYANT, L., ARAKI, T., DENOYES, B., AND FOUCHER, F. The *TFL1* homologue *KSN* is a regulator of continuous flowering in rose and strawberry. *Plant. J.* **69-1** (2012), 116–125.

[18] KAWAMURA, K., HIBRAND-SAINT OYANT, L., CRESPEL, L., THOUROUDE, T., LALANNE, D., AND FOUCHER, F. Quantitative trait loci for flowering time and inflorescence architecture in rose. *Theor. Appl. Genet.* **122-4** (2011), 661–675.

[19] KAWAMURA, K., HIBRAND-SAINT OYANT, L., THOUROUDE, T., JEAUFFRE, J., AND FOUCHER, F. Inheritance of garden rose architecture and its association with flowering behaviour. *Tree. Genet. Genomes.* **11-2** (2015), 1–12.

[20] LYNCH, M., AND WALSH, B. *Genetics and analysis of quantitative traits*. Sinauer, Sunderland, 1998.

[21] McLACHLAN, G., AND BASFORD, K. *Mixture models: inference and applications to clustering*. Dekker, 1988.
[22] McLACHLAN, G., AND PEEL, D. *Finite Mixture Models*. Wiley, 2000.

[23] OGHINA-PAVIE, C. Rose and pear breeding in nineteenth-century France: The practice and science of diversity, in new perspectives on the history of life sciences and agriculture. D. Phillips, and S. Kingsland, Eds. (Springer International Publishing). 53–72.

[24] OTAGAKI, S., OGAWA, Y., HIBRAND-SAINT OYANT, L., FOUCHER, F., KAWAMURA, K., HORIBE, T., AND MATSUMOTO, S. Genotype of *FLOWERING LOCUS T* homologue contributes to flowering time differences in wild and cultivated roses. *Plant. Biology.* 17-4 (2015), 808–815.

[25] PUTTERILL, J., LAURIE, R., AND MACKNIGHT, R. It's time to flower: the genetic control of flowering time. *BioEssays.* 26-4 (2004), 363–373.

[26] ROMAN, H., RAPICAULT, M., MICLOT, A. S., LARENAUDIE, M., KAWAMURA, K., THOUROUDE, T., CHASTELLIER, A., LEMARQUAND, A., DUPUIS, F., FOUCHER, F., LOUSTAU, S., AND HIBRAND-SAINT OYANT, L. Genetic analysis of the flowering date and number of petals in rose. *Tree. Genet. Genomes.* 11-4 (2015), 1–13.

[27] SCHWARZ, G. Estimating the dimension of a model. *Ann. Statist.* 6-2 (1978), 461–464.

[28] SEMENIUK, P. Inheritance of recurrent blooming in *Rosa wichuraiana*. *J. Hered.* 62-3 (1971), 203–204.

[29] SHEKOFTEH, Y., JAFARI, S., CLINTON SPROTT, J., REZA HASHEMI GOLPAYEGANI, M., AND ALMASGANJ, F. A Gaussian mixture model based cost function for parameter estimation of chaotic biological systems. *Commun. Nonlinear. Sci. Numer. Simulat.* 20 (2015), 469–481.

[30] SPILLER, M., LINDE, M., HIBRAND-SAINT OYANT, L., TSAI, C.-J., BYRNE, D. H., SMULDERS, M. J. M., FOUCHER, F., AND DEBENER, T. Towards a unified genetic map for diploid roses. *Theor. Appl. Genet.* 122-3 (2011), 489–500.

[31] WYLIE, A. The history of garden roses. *J. R. Hortic. Soc.* 79 (1954), 555–571.

[32] XU, L., AND JORDAN, M. I. On convergence properties of the EM algorithm for Gaussian mixtures. *Neural. Comput.* 8-1 (1996), 129–151.

E-mail address: Frederic.Proia@univ-angers.fr

E-mail address: Alix.Pernet@angers.inra.fr

E-mail address: Tatiana.Thouroude@angers.inra.fr

E-mail address: Gilles.Michel@angers.inra.fr

E-mail address: Jeremy.Clotault@univ-angers.fr

LABORATOIRE ANGEVIN DE RECHERCHE EN MATHÉMATIQUES – UMR 6093, UNIVERSITÉ D'ANGERS, DÉPARTEMENT DE MATHÉMATIQUES, FACULTÉ DES SCIENCES, 2 BOULEVARD LAVOISIER, 49045 ANGERS CEDEX, FRANCE.

INSTITUT DE RECHERCHE EN HORTICULTURE ET SEMENCES – UMR 1345, INRA, SFR 4207 QUASAV, 42 RUE GEORGES MOREL, 49071 BEAUCOUZÉ CEDEX, FRANCE.

INSTITUT DE RECHERCHE EN HORTICULTURE ET SEMENCES – UMR 1345, UNIVERSITÉ D'ANGERS, SFR 4207 QUASAV, 42 RUE GEORGES MOREL, 49071 BEAUCOUZÉ CEDEX, FRANCE.