Mario Rienzo

Alain Biz

Dionyssis Pongas

St Phanie Michineau

Jonathan Melka

Hon Lai Chan

Lucien Sambin

Jin Bo Su

Jean-Luc Dubois-Rand

Luc Hittinger

Alain Berdeaux

Bijan Ghaleh

Impaired left ventricular function in the presence of preserved ejection in chronic hypertensive conscious pigs

Keywords: Angiotensin II, Animals, Diastole, Female, Hemodynamics, Hypertension, physiopathology, Hypertrophy, Left Ventricular, chemically induced, Myocardial Contraction, Swine, Ventricular Function, Left isovolumic contraction, isovolumic relaxation, hypertension, left ventricular function

Systolic function is often evaluated by measuring ejection fraction and its preservation is often assimilated with the lack of impairment of systolic left ventricular (LV) function. Considering the left ventricle as a muscular pump, we explored LV function during chronic hypertension independently from increased afterload conditions. Fourteen conscious and chronically instrumented pigs received continuous infusion of either angiotensin II (n 8) or saline (n 6) during 28 days. Hemodynamic recordings were = = regularly performed in the presence and 1h after stopping angiotensin II infusion to evaluate intrinsic LV function. Throughout the protocol, mean arterial pressure steadily increased by 55 4 mmHg in angiotensin II-treated animals. There were no significant ± changes in stroke volume, LV fractional shortening or LV wall thickening, indicating the lack of alterations in LV ejection. In contrast, we observed maladaptive changes with 1) the lack of reduction in isovolumic contraction and relaxation durations with heart rate increases, 2) abnormally blunted isovolumic contraction and relaxation responses to dobutamine and 3) a linear correlation between isovolumic contraction and relaxation durations. None of these changes were observed in saline-infused animals. In conclusion, we provide evidence of impaired LV function with concomitant isovolumic contraction and relaxation abnormalities during chronic hypertension while ejection remains preserved and no sign of heart failure is present. The evaluation under unloaded conditions shows intrinsic LV abnormalities.

Introduction

Systolic function is often evaluated by measuring ejection fraction and its preservation is often interpreted as an absence of impairment of systolic left ventricular (LV) function 15 . Usual description of LV function distinguishes systolic from diastolic alterations based on a [ ] classic hemodynamic model of the heart 36 . However it has to be emphasized that the classical so-called diastole and systole are tightly LV function as a hydraulic input-output system as opposed to a muscular pump perspective 7, 23, 38 .

[ ]

In this setting and considering the heart as a muscular pump, we aimed at investigating LV function during chronic hypertension in a large mammal. Previous large animal hypertensive models used renal wrapping 13, 27 or angiotensin II delivered through implanted

[ ]
osmotic pumps 29, 30 . This hypertension was irreversible and the influence of increased afterload could not be easily distinguished

[ ]
per se from intrinsic LV alterations. Moreover there is generally missing information about early phases of LV dysfunction development.

Alternatively, we focused on the development of LV dysfunction during chronic hypertension in the absence of any sign of heart failure and investigated LV contraction, relaxation and filling independently of increased afterload condition. In the present study, hypertension was induced by continuously infusing angiotensin II during 4 weeks with external peristaltic pumps in chronically instrumented conscious pigs. During the whole protocol, hemodynamic status was regularly investigated and, at each scheduled recording, infusion of angiotensin II was stopped for an hour in order to allow arterial blood pressure to return to pre-treatment values. This novel approach allowed both the characterization of LV function during hypertension and the evaluation of the intrinsic LV properties under unloaded conditions, , in the i.e.

absence of angiotensin II infusion.

Methods

The experiments were conducted in accordance with the regulations concerning the use of animals in research.

Surgical instrumentation

Under anesthesia, fourteen female pigs (20 30 kg) underwent left thoracotomy. Fluid-filled Tygon catheters were implanted in the proximal descending thoracic aorta, in the left atrium and in the pulmonary artery for measurement of arterial blood pressure and angiotensin II infusion. A solid-state pressure transducer (P5A, Konigsberg Instruments, Pasadena, CA) was introduced into LV. A flow probe (Transonic Systems Inc, Ithaca, NY) was implanted around the aortic root. Piezoelectric crystals were implanted 1) on opposing LV anterior and posterior endocardial surfaces to measure LV internal diameter, and 2) on opposing LV endocardial and epicardial anterior free wall surfaces to measure wall thickness. All catheters and wires were exteriorized between the scapulae. Animals were sedated with diazepam (0.2 0.4 mg/kg iv) for postoperative care. Enrofloxacine (2 mg/kg/day, im) and long acting amoxicillin (15 mg/kg every two days, im) were administered during the ten days after surgery. Buprenorphine (0.3 mg sc) was administered for five days after surgery. The position of all catheters and crystals was confirmed at autopsy.

Hemodynamic measurements

All hemodynamic data were recorded (1 kHz), digitized and analyzed using HEM v4.2 software (Notocord Systems, Croissy sur Seine, France). Aortic and left atrial pressures were measured with P23XL pressure transducers (Becton-Dickinson, Franklin Lakes, NJ, USA). Cardiac output was measured using a T206 blood flow-meter (Transonic Systems Inc., Ithaca, NY, USA). Left ventricular pressure was cross-calibrated with the left atrial and aortic pressures. The change in LV pressure over time (dP/dt) was computed from the LV pressure signal. Percentage of wall thickening was defined as end-systolic minus end-diastolic thicknesses divided by end-diastolic thickness 100. The signal of LV wall thickness and LV internal diameter were suitable in 11 out of the 14 instrumented pigs.

×

LV end-diastole was defined as the initiation of the upstroke of LV pressure tracing after atrial contraction and indicated by the initial increase in LV dP/dt. Aortic valve opening and closure were identified by crossing LV and aortic pressure waveforms. Furthermore these determinations were confirmed by matching the waveforms with the aortic flow signal. The isovolumic contraction time was defined as the time interval between end-diastole and aortic valve opening. The ejection time was defined as the time interval between aortic valve opening and closure. Mitral valve opening was identified by crossing LV and atrial pressure waveforms. The isovolumic relaxation period was defined as the period elapsed from aortic valve closure and mitral valve opening. The filling period was computed from the end of isovolumic relaxation to end-diastole.

Mean ejection LV wall stress

Cylindrical wall stress was calculated as: stress 1.36 (LVP x ID/2h), where LVP is LV pressure, ID is internal diameter (short axis)

= ×
and is wall thickness. The integral of the systolic wall stress over time, so-called mean ejection wall stress, was calculated during the h ejection period.

Isovolumic relaxation time constant

The time constant of isovolumic LV pressure (LVP) decay was calculated using pressure data points during the isovolumic relaxation τ period and was computed using the logistic method 25, 31 :

[ ]
where P is an amplitude constant and P is the pressure asymptote. The non-linear least squares problem was solved utilizing the A B

Levenberg-Marquardt algorithm.

Protocol

The experiments were conducted 2 3 weeks after surgery when pigs were healthy and apyretic. All animals were studied in the conscious state staying in a sling and care was taken that at least 2 respiratory cycles were included in each recording. A first set of recordings was performed to obtain measurements at Day 0. Then, responses to dobutamine (10 g/kg/min, 5 min) were evaluated.

μ Thereafter, 8 pigs received a continuous angiotensin II infusion (30 ng/kg/min) for 4 weeks using external peristaltic portable pumps. This dose was chosen on the basis of dose-response curves established in preliminary experiments (angiotensin II, 10 to 100 ng/kg/min iv).

Hemodynamic data was regularly recorded at days 3, 7, 14, 21 and 28. Each recording was performed in the presence and 1h after stopping the angiotensin II infusion in order to minimize the impact of changes in loading conditions, , to evaluate the intrinsic LV i.e. contractile properties. Dobutamine infusions were repeated at Day 28 in the absence of angiotensin II.

Each animal served as its own control. In addition, six other pigs received saline infusion and served as matched controls.

Histology 3 13

Slices from the left ventricle were fixed in 4 formalin and embedded in paraffin for histology. Five m thick paraffin-embedded % μ sections were stained with Sirius red for visualization of collagen fibers. Sections were mounted in Eukitt and examined with a bright field microscope. Interstitial fibrosis was quantified and expressed as a percentage of the field area (10 fields analyzed at X10 magnification).

Statistical Analysis

All results are means SEM. Statistical analysis was performed using one-way analyses of variance for repeated measures. When ± overall differences were detected, individual comparisons were performed by Student s -test for paired observations with Bonferroni s ' t ' correction. A value of < 0.05 was considered statistically significant. P

Results

Hemodynamic

Hemodynamic parameters in angiotensin II-and saline-treated pigs are shown in Tables 1 and2, respectively. The angiotensin II-treated animals did not develop any clinical signs of heart failure such as dyspnea, oedema or ascitis.

In angiotensin II-treated animals, arterial blood pressure significantly increased during the 4 weeks of angiotensin II infusion and this effect remained stable (Figure 1A), confirming the constant effectiveness of angiotensin II infusion. The rises in systolic, diastolic and mean blood pressures averaged 61 3, 50 3 and 54 3 mmHg over the 4 weeks, respectively. At Day 28, systolic, diastolic and mean

± ± ±
arterial blood pressures reached 184 5, 134 3 and 158 5 mmHg, respectively (Figure 1B).

± ± ±

In order to evaluate the intrinsic LV function independently from an increase in arterial blood pressure, angiotensin II infusion was stopped for 1h during each recording. Values of systolic, diastolic and mean arterial pressure returned to Day 0 levels (122 3 120 3 ± vs ± mmHg, 85 4 82 3 mmHg and 105 3 101 3 mmHg, respectively) (Figure 1C). Calculation of mean ejection wall stress confirmed ± vs ± ± vs ± that the loading conditions were matched between Day 0 and Day 28 when the infusion of angiotensin II was stopped (247 23 270 39 ± vs ± g/cm , respectively, p NS).

2

=

Heart rate values measured at Day 28 either with or without angiotensin II were significantly increased Day 0 (93 4 and 107 8 vs ± ± vs 80 2 beats/min, respectively).

±

None of these effects were observed in saline treated pigs.

Development of LV hypertrophy

Chronic infusion of angiotensin II induced a progressive and significant increase in LV end diastolic wall thickness (Figure 2A). As compared to age and weight-matched control pigs (n 6), LV to body weight ratio was significantly increased in angiotensin II-treated = animals (Figure 2B), indicating LV hypertrophy. Histological analysis also showed a clear increase in interstitial fibrosis (Figure 2C).

Collagen content evaluated with Sirius red was increased from 6.2 0.7 to 10.2 1.2 in saline-angiotensin II-treated animals

± % ± % vs (p<0.05).

Isovolumic relaxation alterations

We observed abnormal evolution of isovolumic relaxation over the 4-weeks of angiotensin II infusion (Table 3). Along with the increase in heart rate at Day 28, total cycle was significantly reduced by 22 in the absence of angiotensin II. However, isovolumic % relaxation time remained paradoxically unchanged between Day 0 and Day 28, , the ratio of isovolumic relaxation to total systolic time i.e. was significantly increased by 21 . That is, the time devoted to isovolumic relaxation during each minute increased by 36 (8.0 0.3 % % ± vs 5.9 0.1 s per min at Day 28 Day 0, respectively). Similar alterations were observed in the presence of angiotensin II at Day 28 as ± vs compared to Day 0. Thus, regardless of the presence or absence of angiotensin II, the relative time of isovolumic relaxation was abnormally prolonged at Day 28. These changes were not observed in control pigs (Table 4).

We further analyzed the isovolumic relaxation period by calculating the isovolumic time constant . Paradoxically, remained τ τ

unchanged at Day 28 both in the presence and absence of angiotensin II as compared to Day 0 (17 1 and 18 1 17 1 ms), ,

± ± vs ± i.e.
isovolumic relaxation failed to accelerate with the increase in heart rate.

Importantly as increasing heart rate differently influences the various phases of the cardiac cycle, we further analyzed the hemodynamic recordings at different levels of heart rate measured at Day 0 before any angiotensin II infusion. As shown in Figure 3A, the physiological increase in heart rate in these pigs (before starting angiotensin II infusion, n 8) induced significant decreases in isovolumic = relaxation time, as opposed to the lack of decrease after 4 weeks of angiotensin II infusion.

Concerning dobutamine infusion, at Day 0, isovolumic relaxation time decreased by 18 with the normal increase in triple product of % infusion ( 11 ) (Figure 3B). Saline-treated animals exhibited similar response between Day 0 and Day 28 ( 16 and 17 ,

-% -% -% respectively).

Preserved ejection

During the 4 weeks of chronic angiotensin II infusion, stroke volume, LV fractional shortening and LV wall thickening were similar both in the presence and absence of angiotensin II as compared to Day 0 (Figure 4). Concomitantly and due to the rise in heart rate, cardiac output increased at Day 28 as compared to Day 0 (3.5 0.4 L/min and 3.3 0.5 L/min with or without angiotensin II at Day 28, respectively, ± ± 2.7 0.3 L/min at Day 0, p<0.05). vs

±

Altered isovolumic contraction

We next analyzed isovolumic contraction to determine whether it was altered along with isovolumic relaxation. We first focused on time intervals (Table 3). Despite significant decrease in cycle time, absolute values of isovolumic contraction time remained paradoxically unchanged at Day 28 in the absence of angiotensin II as compared to Day 0, , it was significantly increased by 30 when related cycle i.e. % time, respectively. Thus, the time devoted to isovolumic contraction during each minute was significantly increased by 38 (5.9 0.3 4.3 % ± vs 0.2 s per min at Day 28 and Day 0, respectively). These changes were not observed in control pigs (Table 4).

±

At Day 0 (before starting angiotensin II infusion), these parameters were also analyzed at different levels of heart rate. Isovolumic contraction time decreased when heart rate was rising (Figure 5A). This parameter (expressed as of cycle) remained unchanged (8 1, 8

% ± ±
1 and 8 1 of cycle for heart rates at 80 2, 3 1 and 102 1 beats/min, respectively). This contrasts with the changes of these parameters

± % ± ± ±
observed after 4 weeks of angiotensin II infusion.

Finally, we investigated isovolumic contraction using dobutamine infusion. At Day 0 (Figure 5B), isovolumic contraction time decreased by 45 while the triple product of heart rate, LV dP/dt and maximal LV pressure was increased. After 4 weeks of angiotensin % max II infusion (in the absence of angiotensin II), this response was significantly blunted ( 27). This shows a maladaptive response to

-%
inotropic stimulation with a lack of reduction in isovolumic contraction time. Saline-treated animals exhibited similar response between Day 0 and Day 28 ( 42 and 44 , respectively).

-% -%

Diastolic filling alterations

In parallel to the maladaptive response of isovolumic contraction and relaxation, the absolute diastolic filling time was significantly reduced (Table 3). Saline-treated pigs showed similar values between Day 0 and Day 28 (Table 4).

Relationship between isovolumic contraction and relaxation

As illustrated in Figure 6A, there was a linear relationship between the isovolumic relaxation time values ( of cycle time) and % corresponding times of isovolumic contraction ( of cycle time) either in the presence of the absence of angiotensin II. Values were % measured at Days 0, 3, 7, 14, 21 and 28. This indicates that the development of alterations in isovolumic relaxation mirrored those of isovolumic contraction (Figure 6B). The timing of alterations did not differ between contraction and relaxation. For example, isovolumic contraction time and isovolumic relaxation time (expressed as of cycle) were both increased up to 8.9 0.6 and 12.7 0.4 at Day 14, % ± % ± % from 7.2 0.3 and 9.8 0.1 at Day 0, respectively.

± % ± %

Discussion

The present study demonstrates impaired left ventricular function in the presence of a preserved ejection in chronic hypertensive pigs as demonstrated by maladaptative changes in isovolumic contraction and relaxation, regardless of changes in loading conditions. These two periods were tightly coupled and simultaneously altered. Indeed, we observed 1) a paradoxical increase in relative durations of isovolumic contraction and relaxation concomitant with the increase in heart rate over the 4 weeks of angiotensin II infusion, 2) both abnormal isovolumic contraction and relaxation responses to dobutamine and 3) a linear correlation between isovolumic contraction and relaxation durations. The link between contraction and relaxation has been previously described in intact animal and human studies investigating the effects of acute load changes, ischemic heart disease, pulmonary hypertension or systemic chronic hypertension 5, 6, 8, In this study, we investigated our animals during chronic hypertension and did not observe any of the usual signs of heart failure. As previously described 9, 11, 17 , LV hypertrophy correlates with impaired LV isovolumic relaxation and delayed mitral valve opening. In 28, despite the increase in heart rate, was abnormally unchanged as compared to Day 0. This contrasts with the normal heart which τ accelerates LV relaxation at higher rates as shown by our results obtained at Day 0 at various levels of heart rate and by previous reports 4

[ . We further investigated isovolumic relaxation during dobutamine infusion and demonstrated again the inability of the heart to adequately ] adapt to stress.

We next focused on stroke volume, fractional shortening and LV wall systolic thickening. All these parameters were unchanged and cardiac output was not reduced, indicating preserved ejection. However, none of these parameters specifically investigate the isovolumic phase. We therefore turned our attention towards the isovolumic contraction phase searching for alterations independent from ejection. In contrast to the normal heart, isovolumic contraction time remained unchanged. The ratios of isovolumic contraction to total systolic time was paradoxically and significantly increased after 4 weeks of angiotensin II infusion, , the systolic isovolumic contraction failed to i.e. accelerate with the increase in heart rate. One could argue that increasing heart rate would differently influence the various phases of the cardiac cycle and therefore explain our results 16 . This is unlikely as the analysis of these time intervals with different levels of heart rate [ ] in basal conditions (at Day 0 before starting angiotensin II infusion) clearly shows opposite pattern. Moreover in the presence of dobutamine, the isovolumic contraction time failed to adequately decrease. Again, this contrasts with the normal heart in which there is a linear relationship between this pre-ejection period and R-R interval 34 .

[ ]

In the present study, changes in the two LV isovolumic phases were simultaneous and tightly coupled. There was a linear correlation between isovolumic contraction and relaxation under normal and stress conditions with dobutamine. In agreement with the concept that the heart is a muscular pump rather than a hydraulic pump 3 , there was symmetry between isovolumic contraction and relaxation as [ ] illustrated in Figure 6B. One could speculate that increased fibrosis and remodeling 1, 24 could participate to these alterations but further

[ ]
investigations are needed.

Several limitations of the present study should be stressed. First, we did not specifically investigate the cellular mechanisms explaining the alteration in isovolumic contraction. Abnormalities in calcium handling and cell shortening may occur 10, 12, 28, 32 . Second, we did

[ ]
not perform echocardiographic studies to measure myocardial deformation and to evaluate left ventricular asynchrony. One could speculate that inappropriate prolonged isovolumic contraction is linked to increased systolic asynchrony with nonuniform activation. This would induce exacerbated ventricular nonuniformity and diastolic asynchrony with altered isovolumic relaxation. Furthermore, isovolumic contraction and cross-bridge detachment during relaxation are energy-consuming processes. In this regard, increased isovolumic time and slow LV relaxation may further jeopardize LV function during the development of the disease and the reduction of myocardial energy reserve. This pig model of hypertension does not induce heart failure but we believe that it constitutes the basis for future studies per se and the development of new pharmacological strategies. Finally, it remains unknown whether our results are limited to our experimental setting or rather could be generalized. The evaluation of isovolumic contraction should be re-visited in other animal models and/or humans. The use of pressure volume catheters may have missed subtle alterations in isovolumic contraction as the parameters usually investigated with PV loop do not specifically explore this phase and studies in patients with conductance catheter report the use of low frequency of acquisition (250 Hz) which may reduce the sensitivity of the detection 35 .

[ ]

In conclusion, we provide evidence of impaired LV function with isovolumic contraction and relaxation abnormalities during chronic hypertension while ejection remains preserved and no sign of heart failure is present. The evaluation under unloaded conditions shows intrinsic LV abnormalities. No statistical differences were detected.

Acknowledgements:

  21 and it is well known that concomitant alterations of systolic and diastolic ventricular function occur during ischemic heart [ ] disease 22, 26 . Moreover several human studies have reported some alterations in parameters describing systolic function, 2 or long axis systolic DTI velocity 33, 37 while ejection fraction was preserved. Most of these studies analyzed [ ] [ ]

[

  18, 19, 20, 21, 37 . However, they did not distinguish between load influences and intrinsic LV changes nor analyze the consequences of ] stress conditions such as dobutamine.

  model, isovolumic relaxation time remained unchanged despite the increase in heart rate observed over the 4 weeks of angiotensin II infusion. This was accompanied by a reduced diastolic filling time. We also evaluated LV pressure decay during isovolumic relaxation by a mathematical model based on a logistic equation, yielding a constant which shows a stable estimate of relaxation 25, 31 . At Day τ [ ]

Figure 1 Figure 4 Figure 6

 146 Figure 1 Arterial blood pressure in angiotensin II-treated animals A) mean arterial pressure measured at Days 0, 3, 7, 14, 21 and 28 of the protocol, in the presence (open circle) and absence (closed circle) of angiotensin II; B) systolic (SAP), diastolic (DAP) and mean (MAP) arterial blood pressures measured at Day 0 (open bars) and Day 28 (grey bars) in the presence of angiotensin II; C) systolic (SAP), diastolic (DAP) and mean (MAP) arterial blood pressures measured in the absence of angiotensin II, at Day 0 (open bars) and Day 28 (grey bars). , p<0.05, time effect over the 4 weeks of the protocol (ANOVA). , p<0.05 † *

  This work was supported by grants from INSERM (Mario Rienzo, poste d accueil 2008), the Soci t Fran aise d Hypertension Art rielle Elliott EB , Hasumi H , Otani N , Matsuda T , Matsuda R , Kaneko N , Smith GL , Loughrey CM . 2011 ; K201 (JTV-519) alters the spatiotemporal properties of diastolic Ca2 release and the associated diastolic contraction during -adrenergic stimulation in rat ventricular cardiomyocytes

	+	β		. Basic Res Cardiol . 106 : 1009 -22
	10.1007/s00395-011-0218-4
	11 . Falc o-Pires
		'	é é	ç	'	é
	(2010), the Fondation de l Avenir (ET9-529) and the R gion Ile de France (CODDIM). We thank the Laboratoire Roche for providing ' é
	diazepam (Valium ) and Sanofi-Aventis for providing heparin.
	1 . Booysen HL , Norton GR , Opie LH , Woodiwiss AJ . 2012 ; Reverse chamber remodelling following adrenergic-induced advanced cardiac dilatation and pump
	dysfunction . Basic Res Cardiol . 107 : 238 -10.1007/s00395-011-0238-0
	2 . Bruch C , Schmermund A , Marin D , Katz M , Bartel T , Schaar J , Erbel R . 2000 ; Tei-index in patients with mild-to-moderate congestive heart failure . Eur Heart J . 21 :
	1888 -95 10.1053/euhj.2000.2246
	3 . Brutsaert DL . 2006 ; Cardiac dysfunction in heart failure: the cardiologist s love affair with time '	. Prog Cardiovasc Dis . 49 : 157 -81 10.1016/j.pcad.2006.08.010
	4 . Colin P , Ghaleh B , Hittinger L , Monnet X , Slama M , Giudicelli JF , Berdeaux A . 2002 ; Differential effects of heart rate reduction and -blockade on left ventricular β
	relaxation during exercise . Am J Physiol . 282 : H672 -H679 10.1152/ajpheart. 00547. 2001
	5 . 535 -46
	10.1007/s00395-009-0017-3
	6 . D Angelo

® References: Correia-Pinto J , Henriques-Coelho T , Roncon-Albuquerque R , Louren o ç AP , Melo-Rocha G , Vasques-N voa ó F , Gillebert

TC , Leite-Moreira AF . 2009 ; Time course and mechanisms of left ventricular systolic and diastolic dysfunction in monocrotaline-induced pulmonary hypertension . Basic Res Cardiol . 104 : ' R , Shah N , Rubler S . 1975 ; Diastolic time intervals in ischemic and hypertensive heart disease: A comparison of isovolumic relaxation time and rapid filling time with systolic time intervals . Chest . 68 : 56 -61 7 . De Keulenaer GW , Brutsaert DL . 2007 ; Diastolic heart failure: a separate disease or selection bias? . Prog Cardiovasc Dis . 49 : 275 -83 10.1016/j.pcad.2006.08.002 8 . Dodek A , Burg JR , Kloster FE . 1975 ; Systolic time intervals in chronic hypertension: Alterations and response to treatment . Chest . 68 : 51 -55 10.1378/chest.68.1.51 9 . Douglas PS , Berko B , Lesh M , Reichek N . 1989 ; Alterations in diastolic function in response to progressive left ventricular hypertrophy . J Am Coll Cardiol . 13 : 461 -7 10.1016/0735-1097(89)90527-5 10 .

Table 1

 1 Hemodynamic parameters in angiotensin II-treated pigs

	n	Day 0	Day 28
		Ang -	Ang +

± * Values are mean s.e.m. ± Values are mean s.e.m. ± * p<0.05 Day 0 vs

Table 3

 3 Diastolic time intervals in angiotensin II-treated pigs

	Day 0	Day 28
	Ang -	Ang +

±

Diastolic filling time to total cycle time ( )

±

Isovolumic contraction time to total systolic time ( )

Table 4

 4 Time intervals in saline-treated pigs

	Day 0	Day 28
	Saline -	Saline +

± Ejection time to total systolic time ( ) % 68 1 ± 67 1 ± 67 1 ± Isovolumic relaxation time to total systolic time ( ) % Values are mean s.e.m., n 6 ± =

13

13

13

13

13

13

13

13

13


Jeanne Pellerin
email: jeanne.pellerin@gocad.org

Bruno Lévy
email: bruno.levy@inria.fr

Guillaume Caumon
email: guillaume.caumon@gocad.org

A Voronoi-Based Hybrid Meshing Method

In this research note we introduce a method to build volumetric hybrid meshes suitable for finite element computations. The generated mesh is conformal to a set of input surfaces; thin sections are meshed with prisms, pyramids, and tetrahedra, while the remaining of the domain is meshed with Delaunay tetrahedra.

Introduction

In geological modeling, building volumetric finite elements meshes suitable for finite element computations is of crucial importance for many applications such as assessing natural resources or evaluating environmental risks and natural hazards. Geological models are multi-material models that can have a high number of interfaces, these interfaces can come close to each other in certain zones and intersect at low-angles defining thin features that may have an impact on large scale numerical simulations. Meshing such models with Delaunay refinement based methods typically generates a very large number of elements, to avoid it we use prismatic, pyramidal, and tetrahedral elements in these zones.

We propose to adapt the elements of a tetrahedral dominant mesh to the local model configuration by generalizing the Voronoi-Delaunay dual relationship. The volume, defined by a set of triangulated surfaces, is sampled with a given number of points that optimize a function based on Centroidal Voronoi Tessellation (Sect. 3). From the Voronoi diagram of these points we identify the zones that can be meshed with their Delaunay triangulation, and the thin sections that cannot. The latter are filled with triangular prisms, square pyramids, and tetrahedra (Sect. 4). Some examples of hybrid meshes generated by the method are given in Sect. 5.

Goals

The input of our method is a set of triangulated surfaces defining a valid boundary representation model. Among the surfaces, we distinguish those constituting the external boundary, delimiting the zone of interest, from the internal boundaries (optional), that are interfaces between the different materials. The external boundary Fig. 1. Top: 200 points randomly distributed inside 3 nested spheres, their Voronoi diagram, its intersection with the external sphere and the Delaunay triangulation of the points. Bottom: After the optimization of the point positions, the sampling of the volume and the sampling of the spheres are isotropic, and the Delaunay tetrahedra are well-shaped. Note that the density of points is not sufficient to recover the thin layers sandwiched between the spheres. Our method uses prisms and pyramids to recover them (see Fig. 2, 3, and further). must define the interior and the exterior of the model. In the current implementation, we assume that all interfaces cut the model into two parts.

The output of our method is a tet/prism/pyramid mesh that conforms to the geometry of the input boundary surfaces and to their topological relationships. In other words, a high-quality triangular remesh of the boundary surfaces appears as a subset of the facets of the volumetric elements. No preliminary identification of thin sections is needed; the size and resolution of the final mesh are controlled by the number of points (vertices of the final mesh) allocated to sample the model.

Sampling the Domain and the Boundary Surfaces

The first main step of our approach is to sample adequately the model and its boundary surfaces with a given number of points. A sampling, set of points S, can be characterized by its Voronoi diagram that associates to each point p ∈ S a cell containing the points of the space closer to this point p than to any other point of the sampling S. When each point of the sampling is at the centroid of its Voronoi cell, the Voronoi diagram is a Centroidal Voronoi Tessellation (CVT); and the dual Delaunay triangulation has well-shaped isotropic elements [START_REF] Du | Centroidal Voronoi Tesselations: Applications and Algorithms[END_REF] (Fig. 2). When each point of the sampling which Voronoi cell intersects the surfaces is at the centroid of the intersection of its Voronoi cell with the surfaces, we have a Restricted Centroidal Voronoi Tessellation.

Our current strategy to sample the model and its boundaries is to randomly distribute the points in the model, optimize their positions to obtain a CVT using Lloyd relaxation and finally move the points near the boundary to obtain a restricted CVT to the surfaces as in [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted voronoi diagram[END_REF]. Further work is needed to define a sampling quality relevant to the hybrid meshing method we propose. It is a crucial point, since the quality and validity of the final mesh elements depend upon it.

Hybrid Mesh Building

Once the points that define the vertices of the final hybrid mesh have been adequately placed, we build the edges, facets and cells. We first generate the triangular facets that remesh the model boundaries as in [START_REF] Pellerin | Topological Control for Isotropic Remeshing of Non-Manifold Surfaces with Varying Resolution: Application to 3D Structural Models[END_REF]. For each connected component of the intersection between a Voronoi cell and the model boundaries, we add one point; for each intersection between an edge of the Voronoi diagram and the model boundaries we build a triangle linking the points associated to the three neighboring restricted Voronoi cells (Fig. 2).

Then, the prisms filling the thin sections are built between the triangles corresponding to two successive Voronoi edge-model boundary intersections (Fig. 2). Pyramids and tetrahedra are also obtained in the thin sections if the triangles defining the element share one or two vertices. After meshing the thin sections, we build the Delaunay triangulation of the points whose Voronoi cells do not intersect the model boundaries. The last elements to build are those linking the thin sections and the zones filled by tetrahedra. We add one hybrid cell for each vertex of the Voronoi diagram shared by 4 cells among which at least one intersects the model boundaries. We determine successively the vertices, edges, and facets of the cell to build from the intersections between the edges of the Voronoi diagram and the model boundaries. In a simple case the cell is a tetrahedron, but in a more complex configuration a post-processing is necessary to divide that cell into prisms, pyramids and tetrahedra. Fig. 2. Prism building in a thin section between two close spheres. The Voronoi edges (black) and the intersections of the Voronoi cells with the spheres (grey) define the triangles remeshing the spheres and the prisms to build in the thin section.

Results

The first results of our meshing method applied to non-geological models are promising particularly for smooth surfaces; meshes of the three nested spheres (Fig. 3) reproduce boundary layer meshes at different resolutions, thin models like the heptagonal toroid (Fig. 4) can be meshed with relatively few elements. For surfaces with sharp angles, the placement of the points has a crucial importance (Fig. 5), and the results could be improved by changing the point sampling strategy. We need to adapt the evaluation of the point sampling quality to our mesh extraction method and integrate it in the optimization of the point positions. In the current implementation, we check the errors of the meshing algorithm to assess the validity of the final mesh, this is indeed not sufficient and we need to work on a thorough check of the mesh validity and on an evaluation of the quality of its elements.

Fig. 3 .

 3 Fig. 3. Hybrid meshing of the three nested spheres (Fig.2) at different resolutions. The number of points used to sample the model (50; 200; 1000) determines the thin features. At all resolutions, we recover the two thin sections.

Fig. 4 .

 4 Fig. 4. Hybrid meshing of a smooth manifold surface of genus 22 with tetrahedra (grey), prisms (white) and pyramids (black).

Fig. 5 .

 5 Fig. 5. Hybrid meshing of models with sharp features. Thin sections are filled with prisms, pyramids and tetrahedra.

Acknowledgments

This work is supported by the Gocad Consortium (www.gocad.org) and the European Research Council (project GOODSHAPE, ERC-StG-205693).


Julien Jacques
email: julien.jacques@polytech-lille.fr

Cristian Preda
email: cristian.preda@polytech-lille.fr

Functional data clustering: a survey

Keywords: Functional data, Nonparametric clustering, Model-based clustering, Functional principal component analysis Données fonctionnelles, Classification automatique, Méthodes non paramétriques, Modèles génératifs, Analyse en composantes principales fonctionnelles

The main contributions to functional data clustering are reviewed. Most approaches used for clustering functional data are based on the following three methodologies: dimension reduction before clustering, nonparametric methods using specific distances or dissimilarities between curves and model-based clustering methods. These latter assume a probabilistic distribution on either the principal components or coefficients of functional data expansion into a finite dimensional basis of functions. Numerical illustrations as well as a software review are presented.

Introduction

The aim of the cluster analysis is to build homogeneous groups (clusters) of observations representing realisations of some random variable X. Clustering is often used as a preliminary step for data exploration, the goal being to identify particular patterns in data that have some convenient interpretation for the user. In the finite dimensional setting, X is a random vector with values in R p , X = (X 1 , . . . , X p ), p ≥ 1. Earliest methods, such as hierarchical clustering [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] or k-means algorithm [START_REF] Hartigan | Algorithm as 1326 : A k-means clustering algorithm[END_REF] are based on heuristic and geometric procedures. More recently, probabilistic approaches have been introduced to characterize the notion of cluster through their probability density [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF][START_REF] Celeux | Gaussian parsimonious clustering models[END_REF][START_REF] Mclachlan | Finite mixture models[END_REF].

In recent years, researchers concentrated their efforts to solve problems (regression, clustering) when p is large, in his absolute value or with respect to the size of some sample drawn from the distribution of X. The curse of dimensionality was and is still a very active topic. A particular case is that of random variables taking values into an infinite dimensional space, typically a space of functions defined on some continous set T . Thus, data is represented by curves (functional data) and the random variable underlying data is a stochastic process X = {X(t), t ∈ T }. If this type of data was for longtime inaccessible for statistics (because of technological limitations), today it becomes more and more easy to observe, to store and to process large amounts of such data in medicine, economics, chemometrics and many others domains (see [START_REF] Ramsay | Functional data analysis[END_REF] for an overview).

Clustering functional data is generally a difficult task because of the infinite dimensional space that data belong to. The lack of a definition for the probability density of a functional random variable, the defintion of distances or estimation from noisy data are some examples of such difficulties. Different approaches have been proposed along the years. The most popular approach consists of reducing the infinite dimensional problem to a finite one by approximating data with elements from some finite dimensional space. Then, clustering algorithms for finite dimensional data can be performed. On the other hand, nonparametric methods for clustering consist generally in defining specific distances or dissimilarities for functional data and then apply clustering algorithms as hierarchical clustering or k-means. Recently, model-based algorithms for functional data have been developed.

The aim of this paper is to propose a review of these clustering approaches for functional data. It is organized as follows. Section 2 introduces functional data and functional principal component analysis as the main tool for analysing and clustering functional data. Section 3 reviews the different clustering methods for functional data : two-stage methods which reduce the dimension before clustering, nonparametric methods and model-based methods. Section 4 discusses the common problems of selecting the number of clusters and of choosing appropriate representation for functional data approximation. Section 5 presents some software for clustering functional data. The numerical results of the application of some reviewed methods on real data are presented in Section 6. Some open problems related to functional data clustering end the paper.

Functional Data Analysis

Functional data analysis (FDA) extends the classical multivariate methods when data are functions or curves. Some examples of such data are presented in Figure 1 The first contributions to functional data analysis concern the factorial analysis and are mainly based on the Karhunen-Loeve expansion of a second order L 2 -continuous stochastic process [START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Fonctions aléatoires de second ordre[END_REF]. A pioneer work paper on the subject is due to Deville [START_REF] Deville | Méthodes statistiques et numériques de l'analyse harmonique[END_REF] -a one hunderd pages paper in the Annales de l'INSEE -with some applications in economy. In [START_REF] Dauxois | Asymptotic theory for the principal component analysis of a vector random function: some applications to statistical inference[END_REF] and [START_REF] Antoniadis | Joint estimation of the mean and the covariance of a Banach valued Gaussian vector[END_REF] the authors obtained asymptotic results for the elements derived from factorial analysis. The contributions of Besse [START_REF] Besse | Etude descriptive d'un processus[END_REF] and of Saporta [START_REF] Saporta | Méthodes exploratoires d'analyse de données temporelles[END_REF] extends to functional data the principal component analysis, the canonical analysis of two functional variables, the multiple correspondence analysis for functional categorical data and the linear regression on functional data. An important contribution to functional categorical data is due to [START_REF] Boumaza | Contribution a l'étude descriptive d'une fonction aléatoire qualitative[END_REF].

More recently, important contributions to regression models for functional data are due to the research group working on functional statistics in Toulouse (STAPH1 ). Let us remind also the monographs on functional data by Ramsay and Silverman [START_REF] Ramsay | Applied functional data analysis[END_REF][START_REF] Ramsay | Functional data analysis[END_REF] developing theory and applications on functional data, the book of Bosq [START_REF] Bosq | Linear processes in function spaces[END_REF] for modeling dependent functional random variables and the recent book of Ferraty and Vieu [START_REF] Ferraty | Nonparametric functional data analysis[END_REF] on nonparametric models for functional data containing a review of the most recent contributions on this topic.

Functional Data

According to [START_REF] Ferraty | Nonparametric functional data analysis[END_REF], a functional random variable X is a random variable with values in an infinite dimensional space. Then, functional data represents a set of observations {X 1 , . . . , X n } of X.

The underlying model for X i 's is generally an i.i.d. sample of random variables drawn from the same distribution as X.

A well accepted model for this type of data is to consider it as paths of a stochastic process X = {X t } t∈T taking values in a Hilbert space H of functions defined on some set T . Generally, T represents an interval of time, of wavelengths or any other continuous subset of R. We restrict our presentation to the case where H is a space of real-valued functions. For multivariate functional data (elements of H are R p -valued functions, p ≥ 2) the reader can refer for instance to [START_REF] Jacques | Model-based clustering for multivariate functional data[END_REF] for a recent work on multivariate functional data clustering.

The main source of difficulty when dealing with functional data consists in the fact that the observations are supposed to belong to an infinite dimensional space, whereas in practice one only has sampled curves observed into a finite set of time-points. Indeed, it is usual that we only have discrete observations X ij of each sample path X i (t) at a finite set of knots {t ij : j = 1, . . . , m i }.

Because of this, the first step in FDA is often the reconstruction of the functional form of data from discrete observations. The most common solution to this problem is to consider that sample paths belong to a finite dimensional space spanned by some basis of functions (see, for example, [START_REF] Ramsay | Functional data analysis[END_REF]). An alternative way of solving this problem is based on nonparametric smoothing of functions [START_REF] Ferraty | Nonparametric functional data analysis[END_REF].

Let us consider a basis Φ = {φ 1 , . . . , φ L } generating some space of functions in H and assume that X admits the basis expansion

X i (t) = L ℓ=1 α iℓ φ ℓ (t) (2.1)
for some L ∈ N, with α iℓ ∈ R.

The sample paths basis coefficients are estimated from discrete-time observations by using an appropriate numerical method:

If the sample curves are observed without error

X ij = X i (t ij ) j = 1, . . . , m i ,
an interpolation procedure can be used. For example, [START_REF] Escabias | Modeling environmental data by functional principal component logistic regression[END_REF] propose quasi-natural cubic spline interpolation to reconstruct annual temperatures curves from monthly values.

On the other hand, if the functional predictor is observed with error

X ij = X i (t ij ) + ε ij j = 1, . . . , m i ,
least squares smoothing is used after choosing a suitable basis as, for example, trigonometric functions, B-splines or wavelets (see [START_REF] Ramsay | Functional data analysis[END_REF] for a detailed study). In this case, the basis coefficients of each sample path X i (t) are approximated by

α i = (Θ ′ i Θ i ) -1 Θ ′ i Xi , with α i = ( α i1 , . . . , α iL ) ′ , Θ i = (φ ℓ (t ij )) 1≤j≤mi,1≤ℓ≤L and Xi = (X i1 , . . . , X imi ) ′ .

Functional Principal Component Analysis

From the set of functional data {X 1 , . . . , X n }, one can be interested in optimal representation of curves into a function space of reduced dimension. The main tool to answer this request, Functional Principal Component Analysis (FPCA), is presented in this section. Althought the practical interest of FPCA for interpretation and data presentation (graphics), it is one of the main tools considered when clustering functional data.

In order to address this question in a formal way, we need the hypothesis that considers X as a L 2 -continuous stochastic process:

∀t ∈ T , lim h→0 E |X(t + h) -X(t)| 2 = 0.
The L 2 -continuity is a quite general hypothesis, as most of the real data applications satisfy this one.

Let µ = {µ(t) = E[X(t)]} t∈T denotes the mean function X.

The covariance operator V of X:

V : L 2 (T ) → L 2 (T ) f V -→ Vf = T 0 V (•, t)f (t)dt,
is an integral operator with kernel V defined by:

V (s, t) = E [(X(s) -µ(s))(X(t) -µ(t))] , s, t ∈ T .
Under the L 2 -continuity hypothesis, the mean and the covariance function are continuous and the covariance operator V is a Hilbert-Schmidt one (compact, positive and of finite trace).

The spectral analysis of V provides a countable set of positive eigenvalues {λ j } j≥1 associated to an orthonormal basis of eigenfunctions {f j } j≥1 :

Vf j = λ j f j , (2.2)
with λ 1 ≥ λ 2 ≥ . . . and T 0 f j (t)f j ′ (t)dt = 1 if j = j ′ and 0 otherwise.

The principal components {C j } j≥1 of X are random variables defined as the projection of X on the eigenfunctions of V:

C j = T 0 (X(t) -µ(t))f j (t)dt.
The principal components {C j } j≥1 are zero-mean uncorrelated random variables with variance λ j , j ≥ 1. With these definitions, the Karhunen-Loeve expansion [START_REF] Karhunen | Über lineare Methoden in der Wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Fonctions aléatoires de second ordre[END_REF] holds:

X(t) = µ(t) + j≥1 C j f j (t), t ∈ T . (2.3)
Truncating (2.3) at the first q terms one obtains the best approximation in norm L 2 of X(t) by a sum of quasi-deterministic processes [START_REF] Saporta | Méthodes exploratoires d'analyse de données temporelles[END_REF] ,

X (q) (t) = µ(t) + q j=1 C j f j (t), t ∈ T . ( 2
μ(t) = 1 n n i=1 x i (t) and V (s, t) = 1 n -1 n i=1 (x i (s) -μ(s))(x i (t) -μ(t)).
In [START_REF] Deville | Méthodes statistiques et numériques de l'analyse harmonique[END_REF] it has been shown that μ and V converges to µ and V in L 2 -norm with convergences rate of O(n -1/2 ).

As previously discussed, the functional data are generally observed at discrete time points and a common solution to reconstruct the functional form of data is to assume that functional data belong to a finite dimensional space spanned by some basis of functions. Let α i = (α i1 , . . . , α iL ) ′ be the expansion coefficient of the observed curve x i in the basis Φ = {φ 1 , . . . , φ L }, such that:

x i (t) = Φ(t) ′ α i with Φ(t) = (φ 1 (t), . . . , φ L (t)) ′ .
Let à be the n × L-matrix, whose rows are the vectors α ′ i , and M (t) = (x 1 (t), . . . , x n (t)) ′ the vector of the values x i (t) of functions x i at times t ∈ T (1 ≤ i ≤ n). With these notations, we have

M (t) = ÃΦ(t).
(2.5)

Under the basis expansion assumption (2.1), the estimator V of V , for all s, t ∈ T , is given by:

V (s, t) = 1 n -1 (M (s) -μ(s)) ′ (M (t) -μ(t)) = 1 n -1 Φ(s) ′ A ′ AΦ(t), (2.6)
where M (s)μ(s) means that the scalar μ(s) is subtracted to each elements of M (s), and A = (I n -1I n (1/n, . . . , 1/n)) Ã where I n and 1I n are respectively the identity n × n-matrix and the unit column vector of size n. From (2.2) and (2.6), each eigen-function f j belongs to the linear space spanned by the basis Φ:

f j (t) = Φ(t) ′ b j (2.7) with b j = (b j1 , . . . , b jL ) ′ . Using the estimation V of V , the eigen problem (2.2) becomes T 0 V (s, t)f j (t)dt = λ j f j (s),
which, by replacing V (s, t) and f j (s) by their expressions given in (2.6) and (2.7), is equivalent to

1 n -1 Φ(s) ′ A ′ A T 0 Φ(t)Φ(t) ′ dt W b j = λ j Φ(s) ′ b j , (2.8)
where W = T 0 Φ(t)Φ(t) ′ dt is the symmetric L × L matrix of the inner products between the basis functions. Since (2.8) is true for all s, we have:

1 n -1 A ′ AW b j = λ j b j .
By defining u j = W 1/2 b j , the multivariate functional principal component analysis is reduced to the usual PCA of the matrix

1 √ n-1 AW 1/2 : 1 n -1 W 1/2 ′ A ′ AW 1/2 u j = λ j u j .
The coefficient b j , j ≥ 1, of the eigen-function f j are obtained by b j = (W 1/2 ) -1 u j , and the principal component scores, are given by

C j = AW b j j ≥ 1.
Note that the principal components scores C j are also the solutions of the eigenvalues problem:

1 n -1 AW A ′ C j = λ j C j .

Preprocessing functional data

Curves are generally observed at discrete instants of time. For this reason a first step when working with functional data is to reconstruct the functional form of data.

A second important step in functional data analysis is, generally, data registration [42, chap. 7].

It consists in centring and scaling the curves in order to eliminate both phase and amplitude variations into the curve's dataset. But, in our opinion, for clustering purpose registration is not necessarily. Indeed, the amplitude and phase variability of curves can be interesting elements to define clusters. For instance, in the well-known Canadian weather dataset (temperature and precipitation curves for Canadian weather stations [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF][START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF][START_REF] Ramsay | Functional data analysis[END_REF]), the geographical interpretation of the clusters of weather stations is mainly due to amplitude variability. Nevertheless, several works perform curves registration before or simultaneously with clustering [START_REF] Liu | Simultaneous curve registration and clustering for functional data[END_REF][START_REF] Sangalli | k-mean alignment for curve clustering[END_REF] aiming to obtain new clusters which are not related to phase and amplitude variations. But, in this tentative, the conclusion is often the absence of cluster after registration. For instance, the Growth dataset [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF][START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF][START_REF] Ramsay | Functional data analysis[END_REF][START_REF] Tuddenham | Physical growth of california boys and girls from birth to eighteen years[END_REF], which consists of growth curves for girls and boys, is considered by [START_REF] Liu | Simultaneous curve registration and clustering for functional data[END_REF] for a clustering study simultaneously with data registration. The result being the absence of cluster, they failed in retrieving the gender of subjects, contrary to other methods [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF][START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF] which does not perform data registration.

Major functional data clustering approaches

Clustering functional data received particular attention from statisticians in the last decade. We present in this section a classification of the different approaches to functional data clustering into three groups. This classification, illustrated in Figure 2, is described below.

A first approach, quoted as raw-data clustering in Figure 2, consists in using directly discretization of the functions at some time points. This approach is the most simple one, since the functions are generally already observed at discrete instants of time. In this situation, there is no need to reconstruct the functional form of the data. Because of the large size of the discretization, clustering techniques for high-dimensional vectorial data must be used. These techniques are not discussed in this paper, and we refer to [START_REF] Bouveyron | Model-based clustering of high-dimensional data : A review[END_REF] for a complete review on the subject. Thus, the first category of methods discussed in the sequel (Section 3.1) is two-stage methods, which first reduce the dimension of the data and second perform clustering.

The second category concerns nonparametric clustering methods and will be reviewed in Section 3.2. These methods consist generally in using specific distances or dissimilarities between curves Inria combined to classical non-probabilistic clustering algorithms designed for finite-dimensional data.

The third category is model-based clustering techniques which assume a probability distribution underlying data. For functional data the notion of probability density generally does not exist [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF]. Therefore, one can consider models involving probability density for some finite dimensional coefficients describing data. These coefficients can be either coefficients of curves into a basis approximation (splines, wavelets...) or principal components scores resulting from functional principal component analysis of the curves. These methods will be presented in Section 3.3.

Functional Data Discrete observations

x x p p p p p p p p p p p p p

t t i i i i i i i i i i i i i i i ( ( Q Q Q Q Q Q Q Q Q Q Q Q Q Q Q Raw-
( ( Q Q Q Q Q Q Q Q Q Q Q Q Q feature extraction (FPCA) Clustering Clustering
Figure 2: Classification of the different clustering methods for functional data.

Two-stage approaches

The two-stage approaches for functional data clustering consist of a first step, quoted as filtering step in [START_REF] James | Clustering for sparsely sampled functional data[END_REF], in which the dimension of data is reduced, and of a second step in which classical clustering tools for finite dimensional data are used. The reducing dimension step consists generally in approximating the curves into a finite basis of functions. Spline basis [START_REF] Wahba | Spline models for observational data[END_REF] is one of the most common choice because of their optimal properties. For instance, B-splines are considered in [START_REF] Abraham | Unsupervised curve clustering using B-splines[END_REF][START_REF] Rossi | Clustering functional data with the som algorithm[END_REF]. Another dimension reduction technique is the functional principal component analysis (see Section 2.2), for which, from a computational point of view, one needs generally to use also a basis approximation of the curves (see Section 2.3). Functional data being summarized either by their coefficients in a basis of functions or by their first principal component scores, usual clustering algorithms can be used to estimate clusters of functional data. In [START_REF] Abraham | Unsupervised curve clustering using B-splines[END_REF] and [START_REF] Peng | Distance-based clustering of sparsely observed stochastic processes, with applications to online auctions[END_REF] the k-means algorithm is used on B-splines coefficients [START_REF] Abraham | Unsupervised curve clustering using B-splines[END_REF] and on a given number of principal component scores [START_REF] Peng | Distance-based clustering of sparsely observed stochastic processes, with applications to online auctions[END_REF]. In [START_REF] Peng | Distance-based clustering of sparsely observed stochastic processes, with applications to online auctions[END_REF] the number of principal component scores is selected according to the percentage of explained variance, which is an usual criterion in principal component analysis. Let remark also that in [START_REF] Peng | Distance-based clustering of sparsely observed stochastic processes, with applications to online auctions[END_REF], the principal component scores are not directly used but transformed in a low-dimensional space thanks to a multi-dimensional scaling [START_REF] Cox | Multidimensional Scaling[END_REF]. In [START_REF] Rossi | Clustering functional data with the som algorithm[END_REF] and [START_REF] Kayano | Functional Cluster Analysis via Orthonormalized Gaussian Basis Expansions and Its Application[END_REF] an unsupervised neural network, Self-Organised Map [START_REF] Kohonen | Self-Organizing Maps[END_REF], is applied respectively on B-spline and Gaussian coefficient's basis.

Nonparametric approaches

Nonparametric approaches for functional data clustering are divided into two categories: methods who apply usual nonparametric clustering techniques (k-means or hierarchical clustering) with specific distances or dissimilarities, and methods which propose new heuristics or geometric criteria to cluster functional data.

In the first category of methods, several works consider the following measures of proximity between two curves x i and x i ′ :

d ℓ (x i , x i ′ ) = T (x (ℓ) i (t) -x (ℓ) i ′ (t)) 2 dt 1/2 .
where x (ℓ) is the ℓ-th derivative of x. In [START_REF] Ferraty | Nonparametric functional data analysis[END_REF] the authors propose to use hierarchical clustering combined with the distance d 0 -the L 2 -metric -or with the semi-metric d 2 . In [START_REF] Ieva | Multivariate functional clustering for the analysis of ecg curves morphology[END_REF] the k-means algorithm is used with d 0 , d 1 and with (d 2 0 + d 2 1 ) 1/2 . In [START_REF] Tarpey | Clustering functional data[END_REF] the authors investigate the use of d 0 with k-means for Gaussian processes. In particular, they prove that the cluster centres are linear combinations of FPCA eigenfunctions. The same distance d 0 with k-means is considered in [START_REF] Tokushige | Crisp and fuzzy k-means clustering algorithms for multivariate functional data[END_REF] defining time-dependent clustering. These methods are summarized in Table 2.

proximity measure clustering d 0 d 1 d 0 + d 1 d 2 k-means [27, 51, 53] [27] [27] hierarchical clustering [20] [20]
Table 2: Classical nonparametric clustering methods with proximity measures specific to functional data.

Remark: Following the method used to estimate the distance d 0 , nonparametric methods can be assimilated to raw-data clustering or to a two-stage methods. Indeed, if d 0 is approximated Inria using directly the discrete observations of curves -using for instance the function metric.lp() of the fda.usc package for the R software -, nonparametric methods are equivalent to raw-data clustering methods. Similarly, if an approximation of the curves into a finite basis is used to approximate d 0 -with function semimetric.basis() of fda.usc -, nonparametric methods are equivalent to two-stage methods with the same basis approximation.

The second category of nonparametric methods proposes new heuristics to cluster functional. In [START_REF] Hébrail | Exploratory Analysis of Functional Data via Clustering and Optimal Segmentation[END_REF] two dynamic programming algorithms which simultaneously perform clustering and piecewise estimation of the cluster centres are proposed. Recently, [START_REF] Yamamoto | Clustering of functional data in a low-dimensional subspace[END_REF] develops a new procedure to identify optimal clusters of functions and optimal subspaces for clustering, simultaneously. For this purpose, an objective function is defined as the sum of the distances between the observations and their projections plus the distances between the projections and the cluster means (in the projection space). An alternate algorithm is used to optimize the objective function.

Model-based approaches

Model on

Type of model Reference FPCA scores Gaussian (parsimonious sub-models) [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF] Gaussian [START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF][START_REF] Jacques | Model-based clustering for multivariate functional data[END_REF] Gaussian spherical (k-means) [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF] basis expansion coefficients Gaussian (parsimonious) [START_REF] James | Clustering for sparsely sampled functional data[END_REF] Gaussian with regime changes [START_REF] Samé | Model-based clustering and segmentation of times series with changes in regime[END_REF] Bayesian [START_REF] Giacofci | Wavelet-based clustering for mixed-effects functional models in high dimension[END_REF][START_REF] Heard | A quantitative study of gene regulation involved in the immune response of anopheline mosquitoes: an application of Bayesian hierarchical clustering of curves[END_REF][START_REF] Ray | Functional clustering by Bayesian wavelet methods[END_REF] Table 3: Model-based clustering approaches for functional data.

Model-based clustering techniques for functional data are not so straightforward as in the finite-dimensional setting, since the notion of density probability is generally not defined for functional random variable [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF]. Thus, such techniques consists in assuming a density probability on a finite number of parameters describing the curves. But contrary to two-stage methods, in which the estimation of these coefficients is done previously to clustering, these two tasks are performed simultaneously with model-based techniques.

We divide model-based clustering techniques for functional data into two sets of methods, summarized in Table 3: those modelling the FPCA scores and those modelling directly the expansion coefficients in a finite basis of functions.

Model-based functional clustering techniques using principal components modelling

In [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF], an approximation of the notion of probability density for functional random variables is proposed. This approximation is based on the truncation (2.4) of the Karhunen-Loeve expansion, and uses the density of principal components resulting from a FPCA of the curves. After an independence assumption on the principal components (which are uncorrelated), they consider a non-parametric kernel-based density estimation and use it to estimate the mean and the mode of some functional dataset. Using a similar approximation of the notion of density for functional random variables, [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF] and [START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF] assume a Gaussian distribution of the principal components, and define model-based clustering techniques by the mean of the following mixture model [START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF]:

f (q) X (x; θ) = K k=1 π k q k j=1 f Cj |Z k =1 (c jk (x); λ jk ),

RR n 8198

where θ = (π k , λ 1k , . . . , λ q k k ) 1≤k≤K are the model parameters and q k is the order of truncation of the Karhunen-Loeve expansion (2.3), specific to cluster k. The main interest of this model, called funclust by the authors, is in the fact that principal component scores c jk (x) of x are computed per cluster, thanks to an EM-like algorithm, which iteratively computes the conditional probabilities of the curves to belong to each cluster, performs FPCA per cluster by weighting the curves according to these conditional probabilities, and computes the truncation orders q k thanks to the scree-test of Cattell [START_REF] Cattell | The scree test for the number of factors[END_REF]. In [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF], the q k 's are fixed to the maximum number of positive eigenvalues, L, which corresponds to the number of basis functions used in FPCA approximation (see Section 2.3), and some parsimony assumptions on the variance λ jk are considered to define a family of parsimonious sub-models, quoted as funHDDC as an extension of the HDDC method for finite dimensional data [START_REF] Bouveyron | High Dimensional Data Clustering[END_REF]. The choice between these different sub-models is performed thanks to the BIC criterion [START_REF] Schwarz | Estimating the dimension of a model[END_REF].

Previously to this work, [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF] have considered a k-means algorithm based on a distance defined as the L 2 distance between truncations of the Karhunen-Loeve expansion at a given order q k . This model, named k-centres, is a particular case of [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF][START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF] assuming narrowly that the variance λ jk are all equals within and between clusters.

Model-based functional clustering techniques using basis expansion coefficients modelling

To our knowledge, the first model-based clustering algorithm has been proposed in [START_REF] James | Clustering for sparsely sampled functional data[END_REF], under the name fclust. The authors consider that the expansion coefficients of the curves into a spline basis of functions are distributed according to a mixture Gaussian distributions with means µ k , specific to each cluster, and common variance Σ:

α i ∼ N (µ k , Σ).
Contrary to the two-stage approaches, in which the basis expansion coefficients are considered fixed, they are considered as random variable, what allows inter alia to proceed efficiently with sparsely sampled curves. Parsimony assumptions on the cluster means µ k allow to define parsimonious clustering models and low-dimensional graphical representation of the curves. The use of spline basis is convenient when the curves are regular, but are not appropriate for peak-like data as encountered in mass spectrometry for instance. For this reason, [START_REF] Giacofci | Wavelet-based clustering for mixed-effects functional models in high dimension[END_REF] recently proposes a Gaussian model on a wavelet decomposition of the curves, which allows to deal with a wider range of functional shapes than splines.

An interesting approach has also been considered in [START_REF] Samé | Model-based clustering and segmentation of times series with changes in regime[END_REF], by assuming that the curves arise from a mixture of regressions on a basis of polynomial functions, with possible changes in regime at each instant of time. Thus, at each time point t ij , the observation X i (t ij ) is assumed to arise from one of the polynomial regression models specific to the cluster X i belongs to. Some Bayesian models have also been proposed. On the one hand, [START_REF] Heard | A quantitative study of gene regulation involved in the immune response of anopheline mosquitoes: an application of Bayesian hierarchical clustering of curves[END_REF] consider that the expansion coefficients are distributed as follows:

α i |σ k ∼ N (µ, σ k Σ) and σ k ∼ IG(u, v),
where IG is the Inverse-Gamma distribution. On the other hand, [START_REF] Ray | Functional clustering by Bayesian wavelet methods[END_REF] propose a hierarchical Bayesian model assuming further that Σ is modelled by two sets of random variables controlling the sparsity of the wavelets decomposition and a scale effect.

Synthesis

We now present a short synthesis underlying the advantage and disadvantage of each categories of methods.

Inria

The use of raw-data clustering is probably the worst choice since it does not take into account the "time-dependent" structure of data, which is inherent to functional data. The two-stage methods consider the functional nature of the data since the first stage consists of approximating the curves into a finite basis of function. The main weakness of these methods is that the filtering step is done previously to clustering, and then independently of the goal of clustering.

Nonparametric methods have the advantage of their simplicity: these methods are easy to understand and to implement. But their strength is also their weakness, since complex cluster structures can not be efficiently modelled by such approach. For instance, using k-means assumes in particular a common variance structure for each cluster, which is not always a realistic assumption.

In our opinion, the best methods are model-based clustering ones, because they take into account the functional nature of data, they perform simultaneously dimensionality reduction and clustering, and they allow to model complex covariance structure by modelling more or less free covariance operator with more or less parsimony assumptions. Both sub-categories of methods discussed in Section 3.3.2 and 3.3.1 are very efficient, and they both have their own advantages.

On the one hand, model-based approaches built on the modelling of the basis expansion coefficients allow to model the uncertainty due to the approximation of the curve into a finite basis of function, what can be important especially for sparsely sample curve. On the other hand, the model-based approaches built on the principal component modelling define a general framework which can for instance be efficiently extended to the clustering of multivariate curves [START_REF] Jacques | Model-based clustering for multivariate functional data[END_REF] or categorical curves.

Model selection

A common problem to clustering studies is the choice of the number of clusters. We present in Section 4.1 several criterion used for model selection in the functional data framework. A second model selection problem occurs for methods using an approximation of the curves into a finite basis of function, i.e. two-stage methods and model-based ones: the choice of an appropriate basis. Section 4.2 discusses this issue.

Choosing the number of clusters

If classical model selection tools, as BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF], AIC [START_REF] Akaike | A new look at the statistical model identification[END_REF] or ICL [START_REF] Biernacki | Assessing a mixture model for clustering with the inegrated completed likelihood[END_REF] are frequently used in the context of model-based clustering to select the number of clusters (see for instance [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF][START_REF] Giacofci | Wavelet-based clustering for mixed-effects functional models in high dimension[END_REF][START_REF] Samé | Model-based clustering and segmentation of times series with changes in regime[END_REF][START_REF] Serban | Multilevel functional clustering analysis[END_REF]), more specific criteria have also been introduced. First of all, Bayesian model for functional data clustering [START_REF] Heard | A quantitative study of gene regulation involved in the immune response of anopheline mosquitoes: an application of Bayesian hierarchical clustering of curves[END_REF][START_REF] Ray | Functional clustering by Bayesian wavelet methods[END_REF] defines a framework in which the number of clusters can be directly estimated. For instance, [START_REF] Heard | A quantitative study of gene regulation involved in the immune response of anopheline mosquitoes: an application of Bayesian hierarchical clustering of curves[END_REF] considered a uniform prior over the range {1, . . . , n} for the number of clusters, which is then estimated when maximizing the posterior distribution. More empirical criteria have also been used for functional data clustering. In the two-stage clustering method presented in [START_REF] Kayano | Functional Cluster Analysis via Orthonormalized Gaussian Basis Expansions and Its Application[END_REF], the clustering is repeated several times for each number of clusters and that leading to the highest stability of the partition is retained. Even more empirical and very sensitive, [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF][START_REF] Ieva | Multivariate functional clustering for the analysis of ecg curves morphology[END_REF] retain the number of clusters leading to a partition having the best physical interpretation. In [START_REF] James | Clustering for sparsely sampled functional data[END_REF], an original model selection criterion is considered. Proposed initially in [START_REF] Sugar | Finding the number of clusters in a dataset: an informationtheoretic approach[END_REF], this criterion is defined as the averaged Mahalanobis distance between the basis expansion coefficients α i and their closest cluster centre. In [START_REF] Sugar | Finding the number of clusters in a dataset: an informationtheoretic approach[END_REF], it is shown for a large class of mixture distributions that this criterion choose the right number of clusters asymptotically with the dimension (here the number L of basis functions).

Choosing the approximation basis

Almost all clustering algorithms for functional data needs the approximation of the curves into a finite dimensional basis of functions. Therefore, there is a need to choose an appropriate basis and thus, the number of basis functions. In [START_REF] Ramsay | Functional data analysis[END_REF], the authors advise to choose the basis according to the nature of the functional data: for instance, Fourier basis can be suitable for periodic data, whereas spline basis is the most common choice for non-periodic functional data. The other solution is to use less subjective criteria such as penalized likelihood criteria BIC, AIC or ICL.

The reader can for instance to refer to [START_REF] James | Clustering for sparsely sampled functional data[END_REF][START_REF] Samé | Model-based clustering and segmentation of times series with changes in regime[END_REF][START_REF] Serban | Multilevel functional clustering analysis[END_REF] for such use.

Software

Whereas there exist several software solutions for finite dimensional data clustering, the software devoted to functional data clustering is less developed. Under the R software environment, two-stage methods can be performed using for instance the functions kmeans or hclust of the stats package, combined with the distances available from the fda or fda.usc packages.

Alternatively, several recent model-based clustering algorithms have been implemented by their authors and are available under different forms:

R functions for funHDDC [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF] and funclust [START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF] are available from request from their authors. An R package is currently under construction and will be available in 2013 on the CRAN2 website, an R function for fclust [START_REF] James | Clustering for sparsely sampled functional data[END_REF] is available directly from James's webpage, the package curvclust for R [START_REF] Giacofci | Wavelet-based clustering for mixed-effects functional models in high dimension[END_REF] is probably the most finalized tool for curves clustering in R, and implements the wavelets-based methods [START_REF] Giacofci | Wavelet-based clustering for mixed-effects functional models in high dimension[END_REF] described in Section 3.3.2.

A MATLAB toolbox, Curve Clustering Toolbox [START_REF] Gaffney | Probabilistic Curve-Aligned Clustering and Prediction with Mixture Models[END_REF], implements a family of two-stage clustering algorithms combining mixture of Gaussian models with spline or polynomial basis approximation.

Numerical illustration

The evaluation of clustering algorithms is always a difficult task [START_REF] Guyon | Clustering: Science or art[END_REF]. In this review, we only illustrate the ability of the clustering algorithms previously discussed to retrieve the class labels of classification benchmark datasets.

The data

Three real datasets are considered: the Kneading, Growth, and ECG datasets. These three datasets are plotted in Figure 3. The Kneading dataset comes from Danone Vitapole Paris Research Center and concerns the quality of cookies and the relationship with the flour kneading process. The kneading dataset is described in detail in [START_REF] Lévéder | Discrimination de courbes de prétrissage[END_REF]. There are 115 different flours for which the dough resistance is measured during the kneading process for 480 seconds. One obtains [START_REF] Lévéder | Discrimination de courbes de prétrissage[END_REF][START_REF] Preda | PLS classification of functional data[END_REF]. They are known to be hard to discriminate, even for supervised classifiers, partly because of the medium quality class. Taking into account that the resistance of dough is a smooth curve but the observed one is measured with error, and following previous works on this data [START_REF] Lévéder | Discrimination de courbes de prétrissage[END_REF][START_REF] Preda | PLS classification of functional data[END_REF], least squares approximation on a basis of cubic B-spline functions (with 18 knots) is used to reconstruct the true functional form of each sample curve. The Growth dataset comes from the Berkeley growth study [START_REF] Tuddenham | Physical growth of california boys and girls from birth to eighteen years[END_REF] and is available in the fda package of R. In this dataset, the heights of 54 girls and 39 boys were measured at 31 stages, from 1 to 18 years. The goal is to cluster the growth curves and to determine whether the resulting clusters reflect gender differences. The ECG dataset is taken from the UCR Time Series Classification and Clustering website3 . This dataset consists of 200 electrocardiogram from 2 groups of patients sampled at 96 time instants, and has already been studied in [START_REF] Olszewski | Generalized Feature Extraction for Structural Pattern Recognition in Time-Series Data[END_REF]. For these two datasets, the same basis of functions as for the Kneading dataset has been arbitrarily chosen (20 cubic B-splines).

Experimental setup

All the clustering algorithm presented in Section 3 can not tested, in particular because they are not all implemented in a software. The following clustering algorithms for functional data are considered:

two-stage methods:

the classical clustering methods for finite dimensional data considered are k-means, hierarchical clustering and Gaussian Mixture Models (package mclust [START_REF] Banfield | Model-based Gaussian and non-Gaussian clustering[END_REF]) and two methods dedicated to the clustering of high-dimensional data: HDDC [START_REF] Bouveyron | High Dimensional Data Clustering[END_REF] and MixtP-PCA [START_REF] Tipping | Mixtures of principal component analyzers[END_REF],

these methods are applied on the FPCA scores with choice of the number of components with the Cattell scree test, directly on discretizations of the curves at the observation times points, and on the coefficients in the cubic B-spline basis approximation.

non-parametric method: k-means with distance d 0 and d 1 [START_REF] Ieva | Multivariate functional clustering for the analysis of ecg curves morphology[END_REF],

model-based clustering methods: Funclust [START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF], FunHDDC [START_REF] Bouveyron | Model-based clustering of time series in group-specific functional subspaces[END_REF], fclust [START_REF] James | Clustering for sparsely sampled functional data[END_REF], k-centres [START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF] (results for the Growth dataset are available in their paper, but not software allow to proceed with the two other datasets), and curvclust [START_REF] Giacofci | Wavelet-based clustering for mixed-effects functional models in high dimension[END_REF].

The corresponding R codes are given in Appendix A.

Results

The correct classification rates (CCR) according to the known partitions are given in Table 4.

Even if no numerical study can conclude to which method is the best, the present results suggest several comments:

A first comment concerns the use of the different types of clustering methods: two-stage, nonparametric or model-based approaches. Two-stage methods can sometimes perform very well (to estimate the class label), but the main problem is that, in the present unsupervised context, we have no possibility to choose between working with the discrete data, with the spline coefficients or with the FPCA scores. For instance, HDDC and MixtPPCA are very well performing on the Growth dataset using the FPCA scores, but they are very poor using the discrete data or the spline coefficients. If nonparametric methods suffer from a similar limitation, due to the choice of the distance or dissimilarity to use, model-based clustering methods, which also require the choice of an appropriate basis, allow generally to use penalized likelihood criteria such as BIC to evaluate the different basis choices. In that sense, model-based approaches provide more flexible tools for functional data clustering.

Concerning the model-based clustering methods, FunHDDC and Funclust are among the best methods on these datasets. On the contrary, fclust and curvclust lead to relatively poor clustering results. This is probably due to the nature of the data, which are regularly sampled and without peak whereas fclust and curvclust are especially designed for respectively irregularly sampled curves and peak-like data.

Conclusion and future challenge

This paper has presented a review of the main existing algorithms for functional data clustering. A classification of these methods has been proposed, into three main groups: 1. two-stage methods which perform dimension reduction before clustering, 2. nonparametric methods using specific distances or dissimilarities between curves and 3. model-based clustering methods which assume a probabilistic distribution on either the FPCA scores or the coefficients of curves into a basis approximation. A critical analysis has been proposed, which highlights the advantages of model-based clustering methods. Some numerical illustrations and a short software review have also been presented, and the corresponding R codes given in the appendix may help the reader in applying these clustering algorithms to his own data. Literature on functional data clustering generally consider the case of functional data as realizations of a stochastic process X = {X t } t∈T , with X t ∈ R, which is the subject of the present paper. Recently, some authors are interested in the case of multivariate functional data, i.e. X t ∈ R p , in which a path of X is a set of p curves. An example of bivariate functional data is given in [START_REF] Ramsay | Functional data analysis[END_REF] with temperature and precipitation curves of Canadian weather stations. Few works have defined clustering algorithms for such multivariate functional data: [START_REF] Ieva | Multivariate functional clustering for the analysis of ecg curves morphology[END_REF][START_REF] Jacques | Model-based clustering for multivariate functional data[END_REF][START_REF] Yamamoto | Clustering of functional data in a low-dimensional subspace[END_REF]. Another case of interest is qualitative functional data [START_REF] Boumaza | Contribution a l'étude descriptive d'une fonction aléatoire qualitative[END_REF], in which X t lives in a categorical space. The marital status of individuals, the status of some patients with respect to some diseases are some examples of such data. To the best of our knowledge, there are no works considering this type of data in the functional data context and in particular, in the clustering topic.

A R codes for curve clustering

In this appendix are given the R codes used to perform functional data clustering on growth dataset. A.2 Clustering with Funclust and FunHDDC

A.1 Data loading

The corresponding computer code are available from request to their authors. Funclust and FunHDDC can be applied directly on the fd object fdata: >res=funclust(fd,K=K) and >res=fun hddc(fd,K=K,model='AkjBkQkDk') FunHDDC proposing several sub-models, each of one have to be tested -'AkjBkQkDk', 'AkjBQkDk', 'AkBkQkDk','AkBQkDk', 'ABkQkDk', 'ABQkDk'-, and the one leading to the highest BIC criterion is retained (available from res°bic).

For both methods, the clusters are stored in res°cls.

A.3 Clustering with fclust

The corresponding computer code is available from James's webpage.

  : the top Figure (a) plots the evolution of some stock-exchange index is observed during one hour; the bottom Figure (b) presents the knee flexion angle observed over a gait cycle.

  Knee flexion angle (degree) over a complete gait cycle.

Figure 1 :

 1 Figure 1: Some examples of functional data.

Figure 3 :

 3 Figure 3: Kneading, Growth and ECG datasets.

2 . 3 2 . 4 8 3 3 . 2 1 16 7

 2324832116 First, the data have to be stored in a list as follows: > nr=nrow(data) > N = ncol(data) > fdat = list() > fdat°x = as.vector(data) > fdat°curve = rep(1:N,rep(nr,N)) > fdat°timeindex = rep(as.matrix(seq(1,nr,1)),N) > grid = seq(1, nr, length = nr) And then, the clustering can be estimated by: > testfit=fitfclust(data=fdat,grid=grid,K=K) the cluster being available from fclust.pred(testfit)°classA.4 Clustering with curvclustFirst, the values of functional data discretization are registered in a list Y, then transformed in an object of the class CClustData: > library('curvclust') > fdat= list() > for (j in 1:ncol(data)) fdat[[j]] =data[,j] > CCD = new("CClustData",Y=fdat,filter.number=1) Dimension reduction is then performed:> CCDred = getUnionCoef(CCD)The number of clusters is specified in the class CClustO:> CCO = new("CClustO") > CCO["nbclust"] = K > CCO["Gamma2.structure"] = "none"and clustering is performed thanks to the function getFCM: > CCR = getFCM(CCDred,CCO) Inria > summary(CCR) The cluster are finally estimated by maximum a posteriori: > cluster = apply(CCR["Tau"],1,which.max) Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.2 Functional Principal Component Analysis . . . . . . . . . . . . . . . . . . . . . . 5 Computational methods for FPCA . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Preprocessing functional data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Major functional data clustering approaches 8 3.1 Two-stage approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Nonparametric approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.3 Model-based approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3.1 Model-based functional clustering techniques using principal components modelling 11 3.3.2 Model-based functional clustering techniques using basis expansion coefficients modelling 12 3.4 Synthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 4 Model selection 13 4.1 Choosing the number of clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.2 Choosing the approximation basis . . . . . . . . . . . . . . . . . . . . . . . . . . 14 The data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 6.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion and future challenge 16 A R codes for curve clustering 17 A.1 Data loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 A.2 Clustering with Funclust and FunHDDC . . . . . . . . . . . . . . . . . . . . . . . 18 A.3 Clustering with fclust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 A.4 Clustering with curvclust . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 Inria RESEARCH CENTRE LILLE -NORD EUROPE Parc scientifique de la Haute-Borne 40 avenue Halley -Bât A -Park Plaza 59650 Villeneuve d'Ascq Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

  .4) Let {x 1 , . . . , x n } be the observation of the sample {X 1 , . . . , X n }. The estimators for µ(t) and V (s, t), for s, t ∈ T , are:

	2.3 Computational methods for FPCA
	Inria

Table 1

 1

	summarizes these two-stage approaches.
	Let remark that there exist several other approaches developed in specific context. For instance,
	[49] decomposes a dataset of curves using a functional analysis of variance (ANOVA) model: tak-
	ing into account repeated random functions the authors propose a clustering algorithm assuming
	a mixture of Gaussian distributions on the coefficients of the ANOVA model.
		type of basis functions
	clustering	B-spline Gaussian eigenfunctions
	k-means	[1]	[39]
	Self-Organised Map	[44]	[32]

Table 1 :

 1 Summary of two-stage clustering approaches for functional data.

Table 4 :

 4 Correct classification rates (CCR) in percentage for Funclust, FunHDDC (best model according BIC), fclust, kCFC, curvclust and usual non-functional methods on the Kneading, Growth and ECG datasets.

  First, the values of the functional data at the observations time points are loaded in the matrix data, and the true label in the vector cls:

	> splines <-create.bspline.basis(rangeval=c(1, max(t)), nbasis = 20,norder=4)
	> fdata <-Data2fd(data, argvals=t, basisobj=splines)
	The number of clusters is 2 for this dataset:
	> K=2
	> library(fda)
	> data=cbind(matrix(growth°hgtm,31,39),matrix(growth°hgtf,31,54))
	> cls=c(rep(1,39),rep(2,54))
	The functional form is reconstructed using spline basis (for FPCA-based methods), and stored
	in an object of the class fd of the fda package:
	> t=growth°age

RR n 8198

http://univ-tlse3.fr/STAPH Inria

http://cran.r-project.org/ Inria

http://www.cs.ucr.edu/∼eamonn/time series data/ RR n 8198