A M Wubs

E Heuvelink

J A Dieleman

J J Magan

Alain Palloix

F A Van Eeuwijk

QTL analyses on

Keywords: genotype-to-phenotype model, light use efficiency, leaf area index, total shoot biomass

published or not. The documents may come

INTRODUCTION

Main focus in plant breeding programs is development of high yielding cultivars. This includes experiments in different contrasting environments, to evaluate performance of genotypes under different circumstances [START_REF] Van Eeuwijk | Statistical models for genotype by environment data: from conventional ANOVA models to ecophysiological QTL models[END_REF]. Under these conditions, yield often shows large genotype by environment interactions, a consequence of QTL by environment interaction at the level of underlying QTL. This makes direct improvement of yield by selection on QTL for yield difficult.

It has been proposed that the combination of crop growth modelling and QTL analysis might help to overcome this problem [START_REF] Hammer | Models for navigating biological complexity in breeding improved crop plants[END_REF][START_REF] Yin | Role of crop physiology in predicting gene-to-phenotype relationships[END_REF]. In this approach, the complex target trait, for example yield, is dissected into its underlying component traits. These component traits are assumed to have less genotype-by-environment interaction, and therefore less QTL by environment interaction, than the target trait. With the use of a crop simulation model, the target trait can be reconstructed from the component traits. [START_REF] Yin | Coupling estimated effects of QTLs for physiological traits to a crop growth model: predicting yield variation among recombinant inbred lines in barely[END_REF] simulated yield of barley using QTL-fitted values for the component traits partitioning, specific leaf area and pre-and post-flowering periods. However, there was a low correlation between simulated and observed yield and biomass. Other research has focussed on other aspects of crop growth and development, like time to flowering [START_REF] Yin | QTL analysis and QTL-based prediction of flowering phenology in recombinant inbred lines of barley[END_REF][START_REF] Uptmoor | Crop model based QTL analysis across environments and QTL based estimation of time to floral induction and flowering in Brassica oleracea[END_REF] or leaf elongation (Chenu et al., 2009). For these traits, the results of this approach were more successful.

The EU-SPICY project (www.spicyweb.eu) aimed at developing a set of tools for more efficient and faster plant breeding. Pepper was used as a model crop. One of the aims was to develop a genotype-to-phenotype model, in which yield would be reconstructed from underlying components. As a first step into modelling yield, a simple model for the production of total shoot biomass, i.e. total biomass of stem, leaves and fruit, was created, as accurate simulation of yields depends largely on accurate simulation of total biomass production. The steps involved in this process are presented and discussed in this paper.

MATERIALS AND METHODS

Experiments

In 2009, four phenotyping experiments were performed at two locations and in two seasons. The two locations were Wageningen, The Netherlands (latitude 52°N), and El Ejido, Almería, Spain (latitude 36°N). At each location two experiments were conducted, one in spring and one in autumn (Table 1). In these experiments, 149 genotypes of a Capsicum annuum intraspecific Recombinant Inbred Line (RIL) population (6 th generation) of the cross 'Yolo Wonder' x 'Criollo de Morelos 334' (Barchi et al., 2007) were used, as well as both parents and the F 1 .

In the Netherlands, experiments were performed in four Venlo-type greenhouse compartments (12×12 m) with glass cover. Per compartment, there were 44 plots. A plot consisted of eight plants of one genotype in an arrangement of 4×2 plants, from which the inner four plants were used for observations and final destructive harvest. 'Yolo Wonder', CM334 and the F 1 were placed in a single plot in all four compartments. Four sets of three genotypes were placed in two compartments (1-2, 2-3, 3-4, 4-1), where they appeared in two plots per compartment. Of the remaining 137 genotypes, 70 genotypes were placed in one plot in one compartment, while 67 genotypes were placed in one plot in two compartments. Plants had one main stem. Side shoots were stopped at the second internode, i.e. they had 3 leaves and 3 flowers. Plant density was app. 6.4 plants per m 2 . Set points for greenhouse climate were 18/21°C night/day with 16°C from sunset to midnight. CO 2 was supplied up to app. 600 ppm when windows opening was less than 5%. At higher ventilation rates, CO 2 was supplied to a level of 400 ppm. Throughout the experiment, 5 min data of inside temperature, humidity and CO 2 concentration were measured. Radiation was measured outside at 10 minute intervals. Transmissivity of the greenhouse was 62%. Additional assimilation light (SON-T) was given when outside global radiation was below 150 W m -2 ; lamps were switched off when global radiation was above 250 W m -2 with a maximum of 16 h light.

In Spain, experiments were conducted in a greenhouse of 40×60 m with a plastic cover. The 152 genotypes were placed in 2 blocks of 192 plots each (extra plots were filled with dummy plants). Plots consisted of 10 plants in a single row, of which the inner three plants were used for measurements on plant growth and development. Plants had two main stems. In the first experiment, side shoots were stopped at the second internode and flowers were removed from the upper two nodes, i.e. side shoots had 3 leaves and 1 flower. In the second experiment, side shoots were stopped at the second internode as well but no flowers were removed, i.e. the side shoots bore 3 leaves and 3 flowers. Plant density was approx. 3 plants per m 2 . In the first experiment, the heating system was used when temperature was lower than 14°C, whereas in the second experiment the heating system was not used. No CO 2 was supplied. Climate data were registered every 5 min for greenhouse air temperature, humidity, CO 2 concentration, and inside global radiation.

At the planting date of each experiment, 6 plants per genotype were harvested destructively. Per plant, dry weights of leaves and stems were recorded, as well as the leaf area. During the experiment, fruits were harvested (at least 50% red), and their fresh weight (all fruits) and dry weight (selection of fruits) were measured. At the end of the experiment three plants per genotype (NL) or per plot (SP) were harvested destructively, and their leaf area and dry weights of leaves, stems and fruits were measured. Not all genotypes were harvested at the same time, due to difference in plant height.

Description of Simulation Model

A simple LINTUL type simulation model [START_REF] Van Ittersum | On approaches and application of the Wageningen crop models[END_REF] was used to simulate growth of shoot biomass. Total shoot biomass depended on light interception and light use efficiency (Eqn. 1).

*

Derivation of Crop Parameters

For the calculation of the rate of change of leaf area index, the increase in leaf area index between the final and initial destructive harvest was calculated. The time between these two harvests was expressed in thermal time (Eqn. 2). Rate of change of leaf area index was the ratio between the increase of leaf area index and thermal time between initial and final harvest.

Light use efficiency is the biomass produced per unit of intercepted radiation. First, the leaf area index for every day was calculated using the initial leaf area index, the average daily temperature and the rate of change of leaf area index. This resulted in daily values for the leaf area index. Using this, the daily intercepted radiation was calculated. A value of 0.7 for the extinction coefficient k was used for all genotypes. Daily intercepted radiation was summed over the total growth period, resulting in the total intercepted radiation per m 2 . Total shoot biomass was expressed as the total shoot biomass per m 2 , which included the fruits which had been harvested during cultivation. Total shoot biomass production was the increase in shoot biomass between final and initial harvest. The light use efficiency is the ratio between total shoot biomass production and total intercepted radiation.

The two component traits were calculated for each of the 149 genotypes in the four experiments. Rate of change of leaf area index was calculated per plant, whereas light use efficiency was calculated per plot. Using the statistical set-up of the experiment, genotypic means were derived from the data per plant/plot using mixed models, correcting for the position of the plant/plot in the greenhouse. From this analysis, the generalized heritability [START_REF] Oakey | Joint modelling of additive and non-additive genetic line effects in single field trials[END_REF] for the component traits was calculated. Genetic correlations were calculated between the four experiments, as well as correlations between the two component traits within an experiment. Genotypic means, heritability and genetic correlations were also calculated for the target variable total shoot biomass. For this variable, harvest date was included as a covariate in the derivation of genotypic means, to correct for differences in the date of final destructive harvest between the genotypes.

QTL Analysis

The genetic map of the population consisted of 455 markers, mainly SNPs, distributed over 12 linkage groups. For the component traits and the target trait, a multienvironment QTL analysis was performed in GenStat (14 th version), following the methodology described by Boer et al. (2007). After the QTL analysis, estimated values for the component traits and the target trait were obtained for the four environments.

Simulations

Two type of simulations were done: 1) with the genotypic means for the component traits derived from the experiment; 2) with the values from the QTL analysis on the component traits. Simulated shoot biomass was compared to the measured genotypic total shoot biomass at the average harvest date. Explained variance by the simulations was calculated from the correlation coefficient between observed and simulated data.

RESULTS

In all four experiments, the heritability of rate of change of leaf area index was high (Table 2). The heritability for light use efficiency was lower, especially in the Netherlands. Genetic correlations between the four experiments ranged between 0.67 and 0.79 for rate of change of leaf area index, indicating that there was a low genotype-byenvironment interaction. The genetic correlations for light use efficiency ranged between 0.4 and 0.64. However, heritabilities and genetic correlations for the target trait of the study, total shoot biomass, were within the same range for the component traits (heritabilities between 0.56 and 0.87 and genetic correlations between 0.56 and 0.73).

For both component traits, seven QTL were found. Effect of the QTL for rate of change of leaf area index was more directed to the parent CM334, while the effect of the QTL for light use efficiency was more determined by parent 'Yolo Wonder'. Two QTL were the same (within 10 cM) for both component traits while another one was closely located (within 20 cM). In all environments the effect of the majority of the QTL was in the same direction, although the size of the effect differed. The QTL explained 40, 51, 42 and 51% of the observed variance rate of change of leaf area index in N1, N2, S1 and S2, respectively, and 37, 27, 23 and 39% of the variance for light use efficiency in N1, N2, S1 and S2, respectively. Comparing the QTL of the component traits with the seven QTL of the target trait total shoot biomass shows that two QTL of rate of change of leaf area index were exactly the same as for total shoot biomass, while three others were located within 15cM. For light use efficiency, one QTL was the same as for total shoot biomass and one other was located within 15cM. There was only one QTL of total shoot biomass not found in the component traits.

Using the observed component traits to simulate total shoot biomass observed in the four experiments resulted in accurate simulation (Fig. 1A), where 93, 100, 99 and 94% of the observed variation in total shoot biomass is explained by the simulations for N1, N2, S1 and S2, respectively. Using QTL-fitted-values for the component traits in the simulations explained 40, 43, 27 and 34% of the variation observed in the total shoot biomass in N1, N2, S1 and S2 (Fig. 1B). The simulation using component traits based on QTL analyses explained on average more of the variation in observed biomass than a direct QTL analysis on total shoot biomass did, where 30, 29, 34 and 37% was explained by the QTL fit in N1, N2, S1 and S2 respectively.

DISCUSSION

When observed component traits were replaced by their QTL-fitted values, the model explained 27 to 43% of the variation in observed total shoot biomass. In the Netherlands, this was 10 to 14% higher than for the direct QTL analysis on total shoot biomass, but in Spain it was 3 to 7% lower. This result confirms the expectation that estimating the total shoot biomass from the QTL fits of its component traits can give better predictions than directly estimating total shoot biomass from a QTL fit.

The present crop model was very simple and contained fewer component traits than for example the models of [START_REF] Yin | Coupling estimated effects of QTLs for physiological traits to a crop growth model: predicting yield variation among recombinant inbred lines in barely[END_REF] and Chenu et al. (2009). Models used in these kind of exercises don't necessary have to be complex [START_REF] Yin | Role of crop physiology in predicting gene-to-phenotype relationships[END_REF][START_REF] Hammer | Models for navigating biological complexity in breeding improved crop plants[END_REF]. The right balance between detail and parsimony should among others result in stable QTL for the component traits [START_REF] Hammer | Models for navigating biological complexity in breeding improved crop plants[END_REF]. The two component traits showed for most QTL either a consistent positive or a negative effect across environments, and could therefore be regarded as stable. However, this was also true for the target trait total shoot biomass. This is one of the first attempts where the production of assimilates (represented here by LUE) was subjected to QTL analysis. Seven QTL were found for this component trait, but they were not very large. In order to calculate the parameter light use efficiency, the target variable total shoot biomass, is used. This component trait was therefore not independently determined from the target variable. Replacing the light use efficiency by leaf photosynthesis would avoid this, however, obtaining reliable photosynthesis curves for 149 genotypes requires numerous measurements which is not feasible in practice.

Calculation of the component traits takes into account the variation in cultivation period between the genotypes and in environmental factors like temperature and radiation, as the component traits are expressed per degree-day (rate of change of leaf area index) and per unit of radiation (light use efficiency). In contrast, when deriving the genotypic means for total shoot biomass, differences in duration of the cultivation period have to be taken into account, for which a suitable correction has to be made.

The extinction coefficient k, which could be seen as the third component trait, was chosen not to be genotype specific. Determination of this parameter demands light measurements across the canopy, for which the plots were not big enough, or destructive and non-destructive measurements on widely spaced plants [START_REF] Park | Predicting the growth interactions between plants in mixed species stands using a simple mechanistic model[END_REF]. The extinction coefficient depends on the type of light, leaf angle distribution and leaf characteristics (Goudriaan and Van Laar, 1994); the present value is for an average situation.

Concluding shoot biomass can be reasonably well simulated with a crop growth model using QTL based component traits, and this gave better results than directly fitting a QTL model on total shoot biomass. The next challenging step is to extend the model with a function simulating yield, either by a simple partitioning coefficient, or a more complex procedure based on sink strengths [START_REF] Marcelis | Sink strength as a determinant of dry matter partitioning in the whole plant[END_REF], and to validate the method with independent data.

Literature Cited

Barchi, L., Lefebvre, V., Lanterri, S., Nagy, I., Grandbastien, M.A. and Palloix, A. 2007.

A high resolution intra-specific linkage map of pepper (Capsicum annuum L.) and the selection of reduced RILs subsets for fast mapping. Genome 50:51-60. Boer, M.P., Wright, D., Feng, L., Podlich, D.W., Lang, L., Cooper, M. and Van Eeuwijk, F.A. 2007. A mixed-model quantitative trait loci (QTL) analysis for multipleenvironment trial data using environmental covariables for QTL-by-environment interactions, with an example in maize. Genetics 177:1801-1813. Chenu, K., Chapman, S.C., Tardieu, F., McLean, G., Welcker, C. and Hammer, G.L. 2009. Simulating the yield impacts of organ-level quantitative trait loci associated with drought response in maize; a 'gene-to-phentype" modelling approach. Genetics 183:1507-1525. Goudriaan, J. and Van Laar, H.H. 1994. Modelling Potential Crop Growth Processes.

Tables

FiguresFig. 1 .

 1 Figures

  Shoot is the total shoot biomass per m 2 at a specified time. I t is the radiation at day t, k the extinction coefficient and LAI the leaf area index at a specified time. Light use efficiency (LUE) is expressed in g biomass produced per MJ of intercepted radiation. (1

	* 1	*	(1)

*

represents the fraction of intercepted radiation on day t (Beer's law). The model assumes that the light use efficiency and extinction coefficient are constant over time.

Leaf area growth depended only on temperature, and increased linearly with temperature sum (Eqn. 2).

* ∑ (2)

LAI t is the leaf area index on day t, LAI 0 is the initial leaf area index, LAI rate is the increase of leaf area index per degree-day. T i is the average 24 h temperature at day i, and T b is the base temperature below which no development takes place. A base temperature of 10°C was used. Time step of the model is one day.

Table 1 .

 1 Time schedule of the experiments and average temperature and radiation of the four phenotyping experiments. N1 spring experiment in the Netherlands, N2 autumn experiment in the Netherlands, S1 spring experiment in Spain, S2 autumn experiment in Spain. **Global radiation inside the greenhouse, second value in the Netherlands represent PAR from assimilation lights.

		Location Season	Date of sowing	Date of planting	Average date of final	Average air temp.	Radiation (MJ m -2 d -1 )**
						harvest	(°C)
	N1* Netherlands Spring	9 Dec 08	5 Feb 09	26 May 09	20.9	7.5+1.01
	N2 Netherlands Autumn 11 Jun 09	15 Jul 09	22 Sep 09	22.4	8.5+0.69
	S1	Spain	Spring 20 Nov 08 14 Jan 09	8 Jun 09	19.2	11.0
	S2	Spain	Autumn 9 Jun 09	16 Jul 09	20 Nov 09	22.5	7.1

*

Table 2 .

 2 Heritabilities of the total shoot biomass and component traits, and correlations between component traits with an experiment.

			Heritabilities		Correlation
		Rate of change of	Light use	Total shoot	between
		leaf area index	efficiency	biomass	component traits
	N1	0.71	0.52	0.77	0.17
	N2	0.88	0.66	0.56	0.32
	S1	0.90	0.73	0.87	-0.34
	S2	0.94	0.75	0.81	0.00

Proc. IV th IS on HortiModel 2012 Eds.: Weihong Luo et al. Acta Hort. 957, ISHS 2012


Arnaud Guyader

Rémi Guisnel

Fabienne Simonneau

Bernard Rocand

Camille Le Bras

Agnès Grapin

Philippe Chatelet

Stéphane Dussert

Florent Engelmann

Laurence Feugey

First results on cryopreservation by dormant bud technique of a set of Malus and Pyrus cultivars from the INRA Biological Resources Centre

Introduction

The Pip Fruit Biological Resources Centre of INRA (National Institute for Agricultural Research) located in Angers (France) is in charge of the preservation, management, characterization and promotion of traditional and scientific genetic resources of apple (8,493 accessions), pear (825 accessions), quince (64 accessions) and related species. The creation of a cryobank should allow us to optimize the security of long term preservation of these germplasm collections. The development of a cryopreservation program for apple and pear germplasm was initiated in 2010 in the UMR GenHort (Angers) thanks to our participation in the national project CRYOVEG, funded by the French public body IBISA (Biology Infrastructure Health and Agronomy), which aims at developing or optimizing cryopreservation techniques for different plant species maintained in French Biological Resources Centres, and at establishing at the national level a scientific and technical cryopreservation network. Two training periods in USDA NCGRP, Fort Collins USA (25 th January to 29 th January 2010) and in JKI Dresden Germany (1 st February to 5 th February 2010) funded by IBISA and COST respectively allowed us to improve our skills and knowledge about this technique. The aims of this project were two-fold:

• to validate the reference protocol in our experimental environment (equipment, plant material and climatic aspects), • to evaluate the response of different genotypes of our germplasm collections.

Materials and Methods

Plant material

The plant material used was a set of diverse genotypes of Malus and Pyrus from our germplasm collections:

• 15 Malus varieties: ancient (11) and modern (2) varieties of dessert apple, ancient varieties of cider apple (2). • 15 Pyrus varieties: ancient varieties of European pear (12), varieties of nashi (3).

Cryopreservation protocol employed at INRA UMR Genhort

The experimental protocol used in Angers was adapted from reference protocols developed in USDA NCGRP, Fort Collins [START_REF] Towill | Cryopreservation of Malus Germplasm using a winter vegetative bud method: results from 1915 accessions[END_REF][START_REF] Towill | Cryopreservation of dormant buds, in Plant Cryopreservation, a Practical Guide[END_REF]. The different steps of protocol were:

• Sampling of graftwoods: the budsticks were harvested in January 2010 in cold conditions (-3 °C). Three consecutive days of negative temperatures before harvesting were observed, which corresponds to the optimal conditions for the reference protocol. • Storage of graftwoods at -0.5 °C in airtight bag during 4 to 12 weeks.

• Preparing of samples: single nodal sections, 3.5 cm long, with the bud in central position were prepared from budsticks at -0.5 °C. • Desiccation of nodal sections at -5 °C in an incubator (Sanyo ® MIR 254). The objective was to reach 30 % moisture content (fresh weight basis). • Slow cooling: Sections were packaged in plastic tubes for slow cooling in a climatic chamber (Binder ® MK53) at 1 °C/h to -30 °C with an additional step at -30 °C for 24 h. • Storage in plastic tubes in the vapour phase over liquid nitrogen in a liquid nitrogen freezer (Taylor-Wharton ® 750 RS) for 72 h. • Slow rewarming in plastic tubes at +3 °C for 24 h and rehydration in plastic bags filled with moist peat moss at +3 °C for 8 to 15 days, depending on the variety. • Regeneration: grafting of the buds by the chip budding technique on MM106 rootstocks for apple and Kirschensaller rootstocks for pear. For each test date and variety, at least 12 buds from cryopreserved material and 6 buds from fresh material were grafted. The rootstocks were planted in pots in greenhouse for 4 to 6 weeks before grafting.

Regeneration of cryopreserved material in comparison with fresh material from different apple varieties.

Fresh material Cryopreserved material

Figure. 1: Regeneration of cryopreserved material in comparison with fresh material from different apple varieties (percentages were calculated from at least 12 buds from cryopreserved material and 6 buds from fresh material).

Results

Considering all the experiments, the average regeneration percentage after grafting of fresh material was 69.6 % for apple and 54.6 % for pear, while the average regeneration percentage after grafting of cryopreserved material was 31.8 % for apple and 26.9 % for pear.

For apple (Figure . 1), the regeneration percentages of cryopreserved material ranged from 0 to 77.8 %. Three genotypes ('Petit amer', 'Pomme Fougère', 'Rose de Benauge') did not respond to the technique. When using the same protocol in different series of tests at different dates, regeneration after cryopreservation fluctuated between 6.3 % and 66.7 % for 'Golden Delicious' and between 50.0 % and 75.0 % for 'Gala'.

Regeneration of cryopreserved material in comparison with fresh material from different pear varieties. Fresh material Cryopreserved material

Figure. 2:

Regeneration of cryopreserved material in comparison with fresh material from different pear varieties (percentages were calculated from at least 12 buds from cryopreserved material and 6 buds from fresh material).

For pear (Figure . 2), the regeneration percentages of cryopreserved material ranged from 0 to 91.7 %. Five genotypes ('La France', 'Hosui', 'Comice rouge', 'Choju', 'Epine du Mas') did not respond to the technique. When using the same protocol in different series of tests at different dates, regeneration after cryopreservation fluctuated between 11.1 % to 91.7 % for 'Williams'.

Discussion

The dormant bud cryopreservation protocol could be successfully applied in our experimental conditions and on our plant material, both with Pyrus and Malus. These results are very encouraging, especially for Pyrus which is reportedly more recalcitrant to the method. However, the current protocol does not yet guarantee a satisfactory regeneration percentage for all genotypes, nor a satisfactory reproducibility of the results for a given genotype. Several factors, which seem to significantly influence the results have been identified: bud


Madalina Olteanu

Nathalie Villa-Vialaneix

Christine Cierco-Ayrolles

Multiple Kernel Self-Organizing Maps

In a number of real-life applications, the user is interested in analyzing several sources of information together: a graph combined with the additional information known on its nodes, numerical variables measured on individuals and factors describing these individuals... The combination of all sources of information can help him to understand the dataset in its whole better. The present article focuses on such an issue, by using self-organizing maps. The use a kernel version of the algorithm allows us to combine various types of information and automatically tune the data combination. This approach is illustrated on a simulated example.

Introduction

In a number of real-life applications, the user is interested in analyzing several sources of information together: a graph combined with the additional information known on its nodes, numerical variables measured on individuals and factors describing these individuals... The combination of all sources of information can help him to understand the dataset in its whole better. The present article focuses on such an issue, by using self-organizing maps (SOM).

SOM has already been extended to the framework of non numerical data, using various approaches. Historically, median SOM, which consists in replacing the standard computation of the prototypes by an approximation in the original dataset, was introduced first. This principle was used to extend SOM to dissimilarity data in [START_REF] Kohohen | Self-organizing maps of symbol strings[END_REF]. Several variants of the original algorithm have been proposed since: [START_REF] Conan-Guez | Fast algorithm and implementation of dissimilarity self-organizing maps[END_REF][START_REF] Hammer | Topographic processing of relational data[END_REF][START_REF] Hammer | Topographic mapping of dissimilarity data[END_REF][START_REF] Olteanu | On-line relational som for dissimilarity data[END_REF] extended the framework of SOM for dissimilarity data by proposing a more flexible representation for the prototypes, either in batch or in online versions. A closely related approach is to rely on kernels to map the original data into an (implicit) Euclidean space where the standard SOM can be used [START_REF] Barreto | The kernel self organising map[END_REF][START_REF] Andras | Kernel-Kohonen networks[END_REF][START_REF] Boulet | Batch kernel SOM and related laplacian methods for social network analysis[END_REF]. Several kernels have been designed to handle complex data such as strings, nodes in a graph or several graphs [START_REF] Gärtner | Kernel for Structured Data[END_REF].

In the present paper, we propose to extend kernel SOM by using a combination of kernels, each of them dedicated to a particular aspect of the entire dataset. The methodology is described in Section 2 and illustrated in Section 3.

Method

The purpose of the present paper is to design a self-organizing map algorithm in which data coming from various sources can be used as inputs in an optimal combination. Let us suppose that we are given D sets of input data, (x d i ) i=1,...,n, d=1,...,D , all measured on the same individuals (i = 1, . . . , n), and taking values in arbitrary spaces (G d ) d . For example, (x d i ) i may represent pdimensional numerical vectors or factor (non-numerical) variables describing the n individuals. More specifically, (x d i ) i may describe the n nodes in a graph and the kind of relation between them or some text variables naming or describing the data.

Multiple kernel SOM

Suppose also that for each of these D datasets, a kernel is known, i.e., a function

K d : G d × G d → R such that K d is symmetric (∀ z, z ∈ G d , K d (z, z ) = K d (z , z)
) and positive (∀ N ∈ N, ∀ (z j ) j=1,...,N ⊂ G d and ∀ (α j ) j=1,...,N ⊂ R, j,j α j α j K d (z j , z j ) ≥ 0). For example, the linear kernel, the Gaussian kernel or the polynomial kernel, among others, may be used for numerical variables. For nodes of a graph, the heat kernel [START_REF] Kondor | Diffusion kernels on graphs and other discrete structures[END_REF], the commute time kernel [START_REF] Fouss | Random-walk computation of similarities between nodes of a graph, with application to collaborative recommendation[END_REF] or any other regularized version of the Laplacian can be used as kernels [START_REF] Smola | Kernels and regularization on graphs[END_REF], while for textual variables, string kernels, based on the number of occurences of common substrings [START_REF] Watkins | Dynamic alignment kernels[END_REF] may be considered. 1From these D kernels, a new kernel can be defined by computing a convex combination:

K(x i , x i ) = D d=1 α d K d (x d i , x d i ), where α d ∈ [0, 1], d α d = 1 and x i = (x 1 i , . . . , x d i ) ∈ G = G 1 × ... × G d .
Obviously, such an application is also symmetric and positive and thus [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] proves that there is a Hilbert space, (H, ., . ) (called the feature space) and an application φ :

G → H (called the feature map) such that K(x i , x i ) = φ(x i ), φ(x i ) .
In the SOM algorithm, the individuals i = 1, . . . , n have to be clustered into a low dimensional grid made of M neurons, {1, . . . , M }. These neurons are related to each other by a neighborhood relationship, h. Each neuron is also represented by prototypes, (p j ) j (j = 1, . . . , M ) taking values in the previously defined feature space H. Following the general framework of kernel SOM described in [START_REF] Barreto | The kernel self organising map[END_REF][START_REF] Villa | A comparison between dissimilarity SOM and kernel SOM for clustering the vertices of a graph[END_REF][START_REF] Boulet | Batch kernel SOM and related laplacian methods for social network analysis[END_REF], the prototypes are written as a convex combination of the input data in the feature space H: p j = i γ ji φ(x i ). The (squared-)distance between an individual x i and a prototype p j is then computed by

p j -φ(x i ) 2 = K(x i , x i ) -2 n l=1 γ jl K(x i , x l ) + n l,l =1 γ jl γ jl K(x l , x l ). (1)
Online multiple kernel SOM If the (α d ) d are given, the standard kernel SOM aims at optimizing (over (γ ji ) ji ) the following energy function :

E((γ ji ) ji , (α d ) d ) = n i=1 M j=1 h (f (x i ), j) φ α (x i ) -p α j (γ j ) 2 α ,
where f (x i ) ∈ {1, . . . , M } is the index of the cluster to which x i is assigned. The notations φ α , p α j and . α are used to emphasize that these quantities depend on α. The novelty of this paper is to include the optimization of the convex combination of kernels into the on-line algorithm which trains the map. More precisely, a stochastic gradient descent step is added to the original on-line SOM algorithm to optimize the energy E((γ ji ) ji , (α d ) d ) over both (γ ji ) ji and (α d ) d .

To perform the stochastic gradient descent step on the (α d ), the computation of the derivative of

E| xt = M j=1 h (f (x t ), j) φ α (x t ) -p α j (γ j ) 2
α (the contribution of the randomly chosen observation x t to the energy) with respect to α is needed.

As ∂ ∂α d [K(x i , x i )] = K d (x d i , x d i )
, we have that

D d := ∂E| xt ∂α d = M j=1 h (f (x t ), j) K d (x d t , x d t ) -2 n l=1 γ jl K d (x d t , x d l ) + n l,l =1 γ jl γ jl K d (x d l , x d l )  
Following an idea similar to that of [START_REF] Rakotomamonjy | SimpleMKL[END_REF], the SOM is trained by performing, alternatively, the standard steps of the SOM algorithm (i.e., affectation and representation steps) and a gradient descent step for the (α d ) d . The full process is given in Algorithm 1.

Algorithm 1 Multiple online kernel SOM

1: ∀ j = 1, . . . , M and ∀ i = 1, . . . , n, initialize γ 0 ji randomly in R such that γ 0 ji ≥ 0 and

n i=1 γ 0 ji = 1 and ∀ d = 1, . . . , D, initialize α 0 d in [0, 1] such that d α 0 d = 1. 2: for t = 1 → T do 3:
Randomly choose an input x l ∈ {x i } i 4:

Assignment step. Find the unit of the closest prototype

f t (x l ) ← arg min j=1,...,M φ α t-1 (x l ) -p α t-1 j (γ t-1 j ) α t-1 5:
Representation step. Update all the prototypes: ∀ j = 1, . . . , M ,

γ t ji ← γ t-1 ji + µ(t)h(f t (x l ), j) δ li -γ t-1 ji
with δ li = 1 if l = i and 0 otherwise.

6:

Gradient descent step. Update the kernel:

∀ d = 1, . . . , D, α t d ← α t-1 d + ν(t)D t d and K t ← d α t d K d .

7: end for

To ensure that the gradient step respects the constraints on α (α d ≥ 0 and d α d = 1), the following strategy can be applied: first, the gradient

D t d = ∂E t |x t ∂α d d
is reduced and it is projected so that the non-negativity of α is ensured (the reduced gradient projection procedure is described in [START_REF] Luenberger | Linear and Nonlinear Programming[END_REF][START_REF] Bonnans | Optimisation Continue[END_REF][START_REF] Rakotomamonjy | SimpleMKL[END_REF]). This leads to the descent direction D t d . The descent step is decreased with the standard rate of ν 0 /t.

In practice, to ensure similar scales for all kernels, a normalization step is performed before computing the first convex combination: all kernels are scaled to have a unit norm, according to the Frobenius norm. The tuning procedure then automatically finds out the (α d ) that give the best map, i.e., that best organizes the observations according to the global kernel K. As illustrated in [START_REF] Olteanu | Which distance use when extracting typologies in sequence analysis? An application to school to work transitions[END_REF], this can be particularly useful when the data can be described by several kernels, all briging up a different insight on the similarities between individuals, as it is often the case in social sciences (see [START_REF] Abbott | Sequence analysis and optimal matching methods in sociology[END_REF][START_REF] Wu | Some comments on "Sequence analysis and optimal matching methods in sociology, review and prospect[END_REF] for a discussion on this topic).

Finally, the additional cost of tuning the (α d ) is moderate because it is performed online: in practice, we simply multiplied by two the number of steps usually needed to train the map and this choice led to a converging algorithm.

Application

In this section, a simple example is used to test the algorithm and illustrate its behavior in the presence of complementary information. 200 observations, divided into 8 groups (numbered from 1 to 8 in the following), were generated using three different kinds of data: 1/ an unweighted graph, simulated similarly as described in [START_REF] Condon | Algorithms for graph partitioning on the planted partition model[END_REF]. The nodes of the groups 1 to 4 and the nodes of the groups 5 to 8 could not be distinguished in the graph structure: the edges within these two sets of nodes were randomly generated with a probability equal to 0.3. The edges between these two sets of nodes were randomly generated with a probability equal to 0.01; 2/ numerical data coming from two dimensional Gaussian distributions. The variables were simulated by Gaussian vectors with independent components having a variance equal to 0.3 and a mean equal to (0, 0) for the odd groups and to (1, 1) for the even groups; 3/ a factor with 2 levels. Observations of groups 1, 2, 5, and 7 were affected to the first level and observations of the other groups to the second level.

Hence, only the combined knowledge of the three datasets gave access to the eight original groups. The online multiple kernel SOM algorithm was applied to this problem with the commute time kernel for the graph and a Gaussian kernel for the numerical data and for the factor variable (re-coded in disjunctive form). The algorithm was compared with a standard kernel SOM approach using one of the three datasets only. It was also compared to kernel SOM using numerical and factor data or all the three datasets but used as if they were issued from the same dataset with a single Gaussian kernel (when the graph was added to the numerical and factor data, it was under the form of its adjacency matrix). The performance of the different approaches were compared by using the normalized mutual information (NMI) [START_REF] Danon | Comparing community structure identification[END_REF] compared to the original classes. This measure has values between 0 and 1 and is equal to 1 when the two partitions are identical. The results are given in Figure 1.

The multiple kernel SOM clearly outperformed all the other approaches, in-

Conclusion

This article proposes a way to combine multiple sources of heterogenous information in self-organizing maps. Our approach can help to select the most relevant sources of information or the most relevant ways to describe similarities between observations. The computational cost of such an approach is moderate but can however be prohibitive when run with very large datasets or a large number of kernels. A sparse approach is currently under study for improving this aspect of the algorithm.

Fig. 1 :

 1 Fig. 1: NMI (left) of online kernel SOM when using all datasets in the multiple kernel framework or only one dataset or all datasets in a single kernel framework. Example of a map resulting from the use of the algorithm (right): each node represents a neuron of the map and has an area proportional to the number of observations classified inside. The colors display the proportion of the eight original classes in the clusters.

Note that this methodology can also been applied to combine several kernels on the same data (e.g., different kinds of string kernels, encapsulating different features of the texts).


Simone Rossi
email: simone.rossi@unisi.it

Tommaso Lisini Baldi

Marco Aggravi

Monica Ulivelli

David Cioncoloni

Viola Niccolini

Lorenzo Donati

Domenico Prattichizzo

Wearable haptic anklets for gait and freezing improvement in Parkinson's disease: a proof-of-concept study

Keywords: Parkinson's disease, Freezing of gate, Wearable robotics, Gait disturbances, Vibration

Background and objective In a proof-of-concept study, we aimed to verify whether the wearable haptic anklets, a device that delivers personalized suprathreshold alternating exteroceptive stimulation at the anklets on demand, may improve the quality of walking, including the freezing of gate (FOG), in idiopathic Parkinson's disease (PD) patients. The clinical relevance of the presented device as a walking pacemaker to compensate the disturbed locomotion through the generation of a more physiological internal walking rhythm should be verified in a dedicated clinical trial. Methods We tested 15 patients diagnosed as idiopathic PD, during their regular treatment regimen. Patients were evaluated during walking with the device switched on and off, personalized at their most comfortable cadence. Stride velocity, variance, and length, as well as FOG episode duration during walking or turning of 180°, were quantified by an optical high-performance motion capture VICON system. Results The alternating, rhythmic, sensory stimulation significantly improved either walking velocity or reduced inter-stride variance. Effects were more variable on stride length. The significant reduction of FOG episodes' duration correlated with clinical severity of scales rating gate and balance. No safety problems occurred. Conclusions The WEARHAP-PD device, whose Technology Readiness Level (TRL) is 6, significantly improved some walking abilities (walking velocity and stride variance) and reduced the duration of FOG episodes in idiopathic PD patients. Unlike the traditional auditory and visual explicit cues that require the user's allocation of attention for correct functioning, the interaction of the patients with the surrounding environment was preserved, due to the likely implicit processing of haptic stimuli.

Introduction

The alleviation of continuous and episodic gait disturbances, as instability, increased cadence, reduced stride length, and reduced velocity and freezing of gait (FOG) [START_REF] Nieuwboer | Abnormalities of the spatiotemporal characteristics of gait at the onset of freezing in Parkinson's disease[END_REF], represents a major therapeutic challenge, if not an unmet need, in Parkinson's disease (PD) patients. Indeed, these symptoms often remain unresponsive to the best possible dopaminergic medical [START_REF] Nonnekes | Freezing of gait: a practical approach to management[END_REF] and surgical [START_REF] Abboud | Predictors of functional and quality of life outcomes following deep brain stimulation surgery in parkinson's disease patients: disease, patient, and surgical factors[END_REF] treatments, even when other motor symptoms improve [START_REF] Nonnekes | Freezing of gait: a practical approach to management[END_REF]. This worsens the patients' daily quality of life [START_REF] Muslimovic | Determinants of disability and quality of life in mild to moderate Parkinson disease[END_REF], increases the risk of falling and disability, and in turn increases the burden for healthcare systems.

Gait disturbances and FOG are usually positively sensitive to external rhythmic stimulations in most, but not all, patients: both auditory and visual cues can improve walking [START_REF] Spaulding | Cueing and gait improvement among people with Parkinson's disease: a meta-analysis[END_REF][START_REF] Ginis | Cueing for people with Parkinson's disease with freezing of gait: a narrative review of the state-of-the-art and novel perspectives[END_REF] and reduce FOG [START_REF] Ginis | Cueing for people with Parkinson's disease with freezing of gait: a narrative review of the state-of-the-art and novel perspectives[END_REF]. Auditory cueing demonstrated significant improvement of cadence, stride length, and velocity [START_REF] Thaut | Rhythmic auditory stimulation in gait training for Parkinson's disease patients[END_REF][START_REF] Decock | BeatPark: a new wearable device for gait auto-rehabilitation in Parkinson's disease delivering adapted musical stimulation[END_REF].

Simone Rossi and Tommaso Lisini Baldi contributed equally to this work.

In contrast, visual cueing significantly improved stride length only [START_REF] Spaulding | Cueing and gait improvement among people with Parkinson's disease: a meta-analysis[END_REF]. However, such positive effects of auditory and visual cues on walking and FOG are counterbalanced by an underestimated drawback that is inherent to explicit cues: they require the constant allocation of patients' attention to cueing (i.e., if the patient does not listen to or look at the cue, the effect cannot take place). This required active cooperation of the user to the cue's functioning obviously obliges patients to restrict contacts and interactions with the surrounding environment in everyday life, thereby representing an additive potential dangerous factor.

We reasoned that the haptic channel could represent a suitable and safer alternative for external cueing, based on these two premises: (1) vibratory stimulation is effective in transiently improving cardinal motor signs in PD, including step parameters [START_REF] Novak | Effect of step-synchronized vibration stimulation of soles on gait in Parkinson's disease: a pilot study[END_REF][START_REF] King | Short-term effects of vibration therapy on motor impairments in Parkinson's disease[END_REF][START_REF] Nunzio | Alternate rhythmic vibratory stimulation of trunk muscles affects walking cadence and velocity in Parkinson'sdisease[END_REF][START_REF] Winfree | The effect of step-synchronized vibration on patients with Parkinson's disease: case studies on subjects with freezing of gait or an implanted deep brain stimulator[END_REF]; (2) it does not require allocation of attention as auditory or visual cues to be effectively processed. Moreover, we reasoned that (3) a bilateral, alternating, rhythmic stimulation at the ankles can be tuned to a given frequency for cueing a personalized comfortable walking cadence, even in healthy subjects; [START_REF] Muslimovic | Determinants of disability and quality of life in mild to moderate Parkinson disease[END_REF] in case of unpredictable FOG, the patient (or the care giver) can even activate on demand the device by an easily accessible remote Bluetooth control (i.e., a smartphone); [START_REF] Spaulding | Cueing and gait improvement among people with Parkinson's disease: a meta-analysis[END_REF] finally, haptic stimuli applied on lower limbs might directly target the activity of spinal centers for locomotion that are particularly sensitive to proprioceptive and exteroceptive inputs [START_REF] Kuo | The relative role of feedforward and feedback in the control of rhythmic movements[END_REF][START_REF] Dimitrijevic | Evidence for a spinal central pattern generator in humans[END_REF], meanwhile being hardly responsive (if not very indirectly through polysynaptic circuitries) to auditory and visual cueing.

On these premises, we investigated whether the exploited wearable haptic vibrating anklets (WEARHAP-PD device), originally developed and patented (Patent WO2017199171) for providing a gait cadence to a subject [START_REF] Baldi | Human guidance: suggesting walking pace under manual and cognitive load[END_REF][START_REF] Baldi | Wearable haptics for remote social walking[END_REF] and to remotely guiding blind people [START_REF] Scheggi | A remote guidance system for blind and visually impaired people via vibrotactile haptic feedback[END_REF], might be safely used as a sensory cue helping to improve walking abilities-including FOG-in PD patients.

Subjects and methods

Experiment 1: Preparatory findings in healthy subjects

Twelve healthy subjects (6 males; age range 18-65 years; mean 39.7 ± 13.6, non statistically different from the PD patients group) were firstly evaluated in adapting their gait cadence to an external constant rhythm provided by the WEARHAP-PD device. This study group served to verify which could be the most appropriate cadence of alternating vibratory stimulation to be subsequently applied to Parkinsonian patients. We intended "the most appropriate" cadence of stimulation as the one at which subjects automatically align or modify their stride sequence to an external rhythm applied through alternating vibrations at the ankles.

Participants were provided with an Android phone and two haptic interfaces. A pressure sensor, connected to the master anklet, was positioned under the right heel to sample pressure data and extract stride duration, which were then transmitted via Bluetooth to the smartphone and logged. In all the trials, participants walked while wearing headphones reproducing white noise to avoid entrainment due to the sound of vibrations and external stimuli. For each participant, the experimental session started with a preliminary acquisition of selfpaced gait along a 200 m path, to record the user's most comfortable cadence. In the second trial, the haptic interfaces were activated at a random rhythm. The frequency was randomly selected in a range of 0-20% faster than the previously estimated user's comfortable stride duration. No additional instruction or hint was given to the participant. For the following 200 m, we evaluated the strides duration measured by the master haptic interface, and for each stride, we extracted the stride duration error, defined as the difference between the current stride duration and the reference gait cycle.

Experiment 2: Patients

Fifteen idiopathic PD patients regularly followed as outpatients by at least 3 years at the local clinical institution were included in the study. Each subject fulfilled the UK Brain Bank Criteria [START_REF] Hughes | Accuracy of clinical diagnosis of idiopathic Parkinson's disease. a clinicopathological study of 100 cases[END_REF] at the diagnosis. Exclusion criteria were the presence of troublesome dyskinesia (verified by objective examination and by a sum score less than 2 at the items for dyskinesia of the section IV of the revised version of the Unified Parkinson Disease Rating Scale (UPDRS)) [

Standard protocol approvals, registrations, and patient consents

The study was approved by the Local Institutional Ethics Committee, and each patient gave her/his written informed consent to participate. Authorization has been obtained for the consent to disclose any recognizable video material.

Clinical evaluation

Patients were evaluated during their regular therapeutic regimen, because we intended to test the device in a fully ecological situation, comparable with everyday life, and because previous studies failed to indicate significant differences of effects of sensory cues on walking parameters between on and off medication conditions [START_REF] Barthel | The laser shoes. A new ambulatory device to alleviate freezing of gait in Parkinson disease[END_REF]. Each patient underwent a complete neurological examination, including the part III of the revised version of the UPDRS [START_REF] Goetz | Movement Disorder Society-sponsored revision of the unified Parkinson'sdisease rating scale (MDS-UPDRS): scale presentation and clinimetric testing results[END_REF], and the compilation of the new FOG questionnaire (NFOG-Q) [START_REF] Nieuwboer | Reliability of the new freezing of gait questionnaire: agreement between patients with Parkinson'sdisease and their carers[END_REF] and of the Tinetti scale for gait and balance [START_REF] Tinetti | Tinetti balance assessment tool. Fall risk index for elderly patients based on number of chronic dis-abilities[END_REF]. The new FOG is divided in three parts: the first one identifies if the patient is a freezer or not; the second and third parts rate the severity and the impact of the FOG, respectively. Their sum produces the final score.

The Tinetti scale evaluates balance and gait: it includes 16 items, the former 9 addressing balance during sitting, arising from a chair, during standing, turning on place, and during the Romberg's test. The latter 7 items explore gait, and more specifically the height, symmetry and continuity of strides, the trunk deviation the position, and the trajectory of hells. The maximal (worst) global score is 28. Clinical-demographic characteristics of patients are summarized in Table 1.

The WERAHAP-PD device

Tactile vibratory sensitivity is influenced by the spatial location on the body, by the distance between the stimulators, by the frequency of stimulation, and by the age of the user. Previous studies have demonstrated that vibration is better sensed on hairy skin due to its thickness and nerve depth and that vibrotactile stimuli are best detected in bony areas [START_REF] Gemperle | Wearable vibro-tactile display[END_REF]. Due to the aforementioned considerations and since our aim was to design an intuitive and non-obtrusive device which could be easily worn, we concentrated on the exploitation of vibrotactile belts, to be worn on ankles.

From a technical point of view, each device is composed by two cylindrical vibro-motors, independently controlled via an external PC (or a smartphone app) using the Bluetooth communication protocol (see Fig. 1).

The communication is realized with an RN-42 Bluetooth antenna connected to a 3.3 V Arduino pro-mini microcontroller. The wireless connection baud rate is 57,600 bps. The microcontroller installed on the board is used to independently control the activation of each motor and receiving data from an external PC. Note that the proposed vibro-motors are controlled by applying a certain amount of voltage, which determines both frequency and amplitude. Thus, frequency and amplitude cannot be changed independently. Since the user's maximal sensitivity is achieved around 200-300 Hz [START_REF] Riener | Sensor-actuator supported implicit interaction in driver assistance systems[END_REF], two Precision Microdrives Pico Vibe vibration motors are used to generate a 250 Hz vibration. The motors have a vibration frequency range of 100-280 Hz, lag time of about 20 ms, and rise and stop time of 35 ms. The anklet guarantees about 4 h of battery life with one motor always turned on. Each device weighs about 90 g.

Tuning the WEARHAP-PD device

The most comfortable walking velocity in healthy subjects aged 50-70 years is on average 1.3 m/s. Under normal conditions, each stride-calculated as the time elapsed from two consecutive heel supports of the same feet-has a length corresponding to about 80% of the subject's height. For example, in a person 1.7 m tall, the stride measures 1.35 m: this means that if the individual walks at 1.3 m/s, it will take 1.05 s to complete a stride.

The left/right alternating vibrotactile stimulus at the ankles is applied with an initial frequency according to the average walking speed of the patient. After having tuned the intensity of the stimulation well above the perceptive threshold, the timing of left/right alternating stimulation is individually tuned to the most comfortable walking speed, as subjectively reported by the patient and as previously verified in the experiment 1. In our sample of patients, these most comfortable velocities ranged from 1.6 to 1 s per step.

Quantitative and qualitative walking assessment

A motion caption system (Vicon Motion Systems Ltd., UK) was used for obtaining gait data regarding the walking patterns of the patients. Each user wore 39 retroreflective markers (Fig. 2), which were labeled using the software Nexus version 2.7 (Vicon Motion Systems Ltd., UK). With Nexus, we were able to record and identify the walking pattern of each user, with a precision on the position of the markers of 0.1 mm. The trials took place in an available area of 12.20 × 5.70 m, where users were asked to walk in the middle, parallel to the long side. In our study, we considered three conditions for the evaluation: no-haptic suggestions mode, haptic-suggestions mode, and mixed mode. In the former condition, users walked from one side to the other of the room four times (two goes and two returns, ending the trial at its initial point), while no haptic suggestion was sent. In the haptic condition, the users were required to act as previously, but this time haptic stimuli, provided by the WERAHAP-PD devices, were used to give a timing information on the gait. In the latter case, the number of trials (walking from one side to another and come back) was the same, but we used haptic stimuli in the first two walking and no cueing for the remaining two. Each user performed 12 walking trials (four per condition).

For each step, we considered the (a) elevation from the ground of the mean for the three markers composing the feet (Fig. 2) and (b) its trajectory (position in time) along the space for the all trial; in this way, we obtain a quantitative measure of how much each foot was going high and how long were the steps of each patient. Velocity and inter-stride variance were also computed. FOG time was evaluated as well, either when occurring unexpectedly during walking or when patients had to turn by 180 degrees at each corner of the room.

Data availability statement

All data, including those not published within the article, are available on request from any qualified investigator to the first two authors. Data will eventually be provided in anonymized form. Data will be shared by request.

Results

Experiment 1: Preparatory findings in healthy subjects

A Spearman's rank-order correlation was run to assess the relationship between the capability of subjects to automatically align the cadence with the proposed rhythm and the distance of the suggested cadence with respect to the user's most comfortable gait. Thus, we measured the capability in synchronization as the difference between the suggested cadence and the actual one of participants.

Analyses showed the relationship to be monotonic, as assessed by visual inspection of a scatterplot. Both variables (i.e., difference between the suggested and most comfortable cadence, real user's gait) were not normally distributed, as assessed by Shapiro-Wilk test (p < .05), and there were no outliers. There was a statistically significant, strong negative correlation between the two variables, r(12) = -0.907, p < 0.0005.

Thus, results showed that subjects automatically aligned to the walking pace only if the rhythm of the vibratory stimulation was close to subjects' preferred walking cadence, while the influence of the vibration was negligible and did not affect the gait when the rhythm was far from the most comfortable one. This is a common behavior, extensively discussed in the literature for different types of cueing [START_REF] Roerdink | Walking to the beat of different drums: practical implications for the use of acoustic rhythms in gait rehabilitation[END_REF][START_REF] Ahn | Walking is not like reaching: evidence from periodic mechanical perturbations[END_REF][START_REF] Leman | Activating and relaxing music entrains the speed of beat synchronized walking[END_REF][START_REF] Damm | Why do we move to the beat? A multi-scale approach, from physical principles to brain dynamics[END_REF] and is the reason why we chose the "most comfortable cadence" in Parkinsonian patients.

Experiment 2: Patients

Table 1 summarizes clinical characteristics of patients, including current therapy. As interviewed at experimental debriefing, all the patients were comfortable with the device and subjectively reported no difficulty to use it. No particular safety concerns emerged, including the absence of any minor side effect.

Stride kinematics

Walking velocity Inspection of Vicon motion system video recordings was carried out with the WEARHAP-PD device switched on and off, for comparable epochs in the two conditions. Stride velocity in m/s increased in 11/15 patients (mean increase 4.5%; upper range 10%) (Fig. 3a) when the device was switched on. In two patients (number 8, 10), instead, the stride velocity decreased: at a closer inspection of clinical data, especially in relation to an unsatisfactory levodopa response, these patients underwent brain MRI during the study period that showed numerous ischemic lesions in the basal ganglia (that were not present at the time of diagnosis), possibly complicating the initial diagnosis of idiopathic PD. Thus, these data were excluded in the further statistical analysis.

Mean stride velocity data in the remaining 13 patients were normally distributed (Shapiro-Wilk test, p = 0.149). The increase of velocity induced by the alternating vibratory stimulation was overall significant (paired sample Student t test = -2.337, p = 0.038). Such increase was even more significant after exclusion of the two other patients (numbers 2 and 14, see Fig. 3a) who were considered non-responder.

Stride length and height As shown in (Fig. 3c), effects of the WEARHAP-PD device on the stride length were variable: increase in 7, decrease in 6 (including patients 8 and 10), and no change in two. The lack of stride length increase in 8 patients (that, however, increased walking velocity as previously said) is simply consequent to the alignment to their most comfortable cadence as tuned by the WEARHAP-PD devices. Stride height did not change significantly with the device switched on (data not shown).

The two patients (8 and 10) with vascular lesions in the basal ganglia worsened stride length and height, besides velocity, with the device on, thus suggesting that vascular lesions might impair a correct supraspinal processing of the haptic cueing.

Inter-stride variance In normal subjects, the inter-stride variance is below 1-2%, while in PD patients, it is increased due to the lack of rhythmicity [START_REF] Ginis | Cueing for people with Parkinson's disease with freezing of gait: a narrative review of the state-of-the-art and novel perspectives[END_REF]. The goal of the WEARHAP-PD device is to act as a walking pacemaker to compensate the disturbed locomotion through the generation of a new internal walking rhythm. To verify whether the device regulates inter-stride variance, in each patient, we measured the average length of left and right strides in every trial, with the device switched on and off. Figure 3b shows that the inter-stride standard deviation (STD) with the device on is reduced by 10% on average, thus suggesting that the alternating haptic stimulation regularized stride length.

To properly verify the capability of our system, we took into consideration the reduction of STD percentage with the respect of the average stride length. Individual values are graphically depicted in Fig. 3b.

Paired sample Student t tests were conducted to assess differences in walking with or without the vibrotactile cueing. Three tests were run in the sample: (1) for the entire patients' group, (2) after exclusion of subjects 8 and 10 (patients with vascular lesions), and (3) after further exclusion of patients 2 and 14 (non-responders for walking velocity improvement). It is worth pointing out that subjects 8 and 14 greatly improved their walking ability. All the tests showed a statistically significant decrease in inter-stride STD walking with the active anklets versus the off condition: test (1) t (14) = 2.319, p = 0.039;

(2) t (12) = 2.188 p = 0.048 p < 0.05; and (3) t (10) = 2.236 p = 0.049 < 0.05.

FOG time reduction

Of the 15 patients, 6 had FOG. In each patient, the duration of each FOG episode was measured during the WEARAP-PD device switched on and off: as shown in Fig. 3d, the FOG duration decreased in all patients, ranging from 2.5 to 16.2 s with the devices in off and from 1.8 to 8.2 s with the devices switched on. This corresponds to a percentage improvement ranging 19.5-49.9%. This difference was significant (p < 0.05, Student t test), after verifying the normality of distribution of FOG episode durations (Shapiro-Wilk, t = 3.009, p = 0.16 > 0.05).

Correlations There was a significant correlation between the severity of the overall Tinetti's clinical scale with the percentage reduction of FOG time (Pearson, r = -0.679; p = 0.05), as well as with the Tinetti rating gate (r = 0.843; p = 0.035).

To prove and validate results on FOG time reduction, a statistical analysis test was conducted on the data, analyzing the time in recovering from a FOG episode (dependent variable) under the two different conditions, vibrating anklets switched on or off (independent variable). Pearson's partial correlation was run to assess the relationship between rates in the NFOG questionnaire and percentage reduction of time of FOG episodes.

There was a linear relationship between the two variables. Both variables were normally distributed as assessed by Shapiro-Wilk test (p = 0.880 and p = 0.098 for questionnaire answers and time percentage reduction, respectively). A bivariate Pearson's correlation showed a statistically significant linear relationship (r = 0.858, p = 0.029).

Discussion

Traditional auditory and visual cues are defined as temporal or spatial stimuli able to regulate and facilitate repetitive movements, as walking, by providing an explicit motor target. They obviously require constant allocation of attention for their correct functioning, thereby distracting users from the surrounding environment. Previous pilot studies of alternating vibratory stimulation for walking improvement in PD used stimuli applied on the soles, time-locked with the foot stand on the ground [START_REF] Novak | Effect of step-synchronized vibration stimulation of soles on gait in Parkinson's disease: a pilot study[END_REF][START_REF] Winfree | The effect of step-synchronized vibration on patients with Parkinson's disease: case studies on subjects with freezing of gait or an implanted deep brain stimulator[END_REF], with the aim to improve the defective proprioceptive feedback usually found in PD [START_REF] Abbruzzese | Sensorimotor integration in movement disorders[END_REF]. When positioned at the ankle, the vibrating device is in very close contact with the rigid surface of the leg bones, and therefore, vibrations are more likely to be perceived than when applied to the soles [START_REF] Novak | Effect of step-synchronized vibration stimulation of soles on gait in Parkinson's disease: a pilot study[END_REF][START_REF] Winfree | The effect of step-synchronized vibration on patients with Parkinson's disease: case studies on subjects with freezing of gait or an implanted deep brain stimulator[END_REF]. Moreover, daily clinical practice indicates that sensory disturbances are usually more prominent in most distal limbs segments and that an overt or subclinical neuropathy is very common in PD, ranging from 16% if only large nervous fibers are considered to 56% if also small fibers are taken into account [START_REF] Zis | Peripheral neuropathy in idiopathic Parkinson's disease: a systematic review[END_REF]. Thus, we reasoned that positioning the vibratory device at the ankles, rather than at the soles, will render more efficiently the processing of the haptic inputs in most of subjects.

Moreover, rather than improving foot stand-related proprioceptive feedback, the novel aim of WEARHAP-PD device is to provide the haptic cueing in a feed-forward mode, as a walking pacemaker tuned at the most comfortable subjective walking speed: indeed, manipulating in real time the dynamical environment during walking is commonly used to mechanically entrain gait patterns [START_REF] Baldi | Wearable haptics for remote social walking[END_REF][START_REF] Ahn | Walking is not like reaching: evidence from periodic mechanical perturbations[END_REF]. As such, the device can be easily tuned by patients or caregivers (via a smartphone app) according to their requirements in everyday life. Such non-continuous utilization of the device, together with the possibility to easily modify stimulation parameters, will minimize habituation to cueing in everyday life [START_REF] Ginis | Cueing for people with Parkinson's disease with freezing of gait: a narrative review of the state-of-the-art and novel perspectives[END_REF]. However, this aspect has been not specifically addressed in the current laboratory investigation, and it will require a prolonged use of the WEARHAP-PD device, in a domestic context.

Locomotion is a voluntarily initiated repetitive movement and then maintained by oscillating interneuronal pools located in the zona intermedia of the gray matter of the spinal cord at L2-L4 myelomers, called central pattern generators (CPG) [START_REF] Kuo | The relative role of feedforward and feedback in the control of rhythmic movements[END_REF]. CPGs may generate their oscillatory neural activity even when the spinal cord is lesioned rostrally [START_REF] Winfree | The effect of step-synchronized vibration on patients with Parkinson's disease: case studies on subjects with freezing of gait or an implanted deep brain stimulator[END_REF]. CPGs that are silent at rest are under direct control from brainstem centers via reticulospinal neurons, which integrate afferent sensory feedback [START_REF] Zehr | Regulation of arm and leg movement during human locomotion[END_REF][START_REF] Nielsen | How we walk: central control of muscle activity during human walking[END_REF]: indeed, they are known to be extremely sensitive to large diameter sensory afferent fiber inputs, as those proprioceptive (neuromuscular spindles, Golgi tendon organs) and exteroceptive (tactile mechanoreceptors) ones continuously activated during walking [START_REF] Kuo | The relative role of feedforward and feedback in the control of rhythmic movements[END_REF]. The alternating suprathreshold haptic input generated by the WEARHAP-PD device is likely to activate most of all large diameter sensory afferents channel mediated by the posterior columns of the spinal cord: as the stimulation is tuned according to "the most comfortable walking cadence", it is likely that the artificial afferent input might favor the reactivation of more physiological pattern of oscillatory neural activity of CPGs, either acting at segmental or suprasegmental levels. This might theoretically explain either the improvement of walking parameters or the improvement/resolving of the FOG.

Walking is a complex function, in which diffuse supraspinal networks involving basal ganglia, thalamus, supplementary motor area, and cerebellum are key structures in the functional loop implementing such rhythmical automatized movements [START_REF] Kotz | Cortical speech processing unplugged: a timely subcortico-cortical framework[END_REF][START_REF] Kotz | Differential input of the supplementary motor area to a dedicated temporal processing network: functional and clinical implications[END_REF]. Basal ganglia obviously play a major and complex role in the pathophysiology of gait problems, including FOG, in PD patients [START_REF] Nutt | Freezing of gait: moving forward on a mysterious clinical phenomenon[END_REF]: one of their function is to correctly process sensory inputs [START_REF] Boecker | Sensory processing in Parkinson's and Huntington's disease: investigations with 3D H(2)(15)O-PET[END_REF][START_REF] Ulivelli | Time course of frontal somatosensory evoked potentials relation to L-dopa plasma levels and motor performance in PD[END_REF]. An inadequate integration of sensory inputs at the striatum and a defective proprioceptive feedback in PD contributes to the alteration of sensorimotor integration and movement control [START_REF] Abbruzzese | Sensorimotor integration in movement disorders[END_REF][START_REF] Lewis | Altered sensorimotor integration in Parkinson's disease[END_REF]. Functional neuroimaging investigations in PD patients have been obviously focused only on cortical and subcortical regional activity, and they point to an "imbalance between subcortical and cortical brain activation underlying the freezing episodes" [START_REF] Fasano | Neuroimaging of freezing of gait[END_REF]; moreover, they still do not offer univocal explanations of the effects of external cueing on neural processing [START_REF] Ginis | Cueing for people with Parkinson's disease with freezing of gait: a narrative review of the state-of-the-art and novel perspectives[END_REF][START_REF] Shik | Control of walking and running by means of electrical stimulation of the mid-brain[END_REF]. It cannot be excluded that the alternating vibratory input generated by the WEARHAP-PD device might additionally improve basal ganglia sensory processing and, in turn, restore a more physiological cortical descending drive toward supraspinal walking centers (i.e., mesencephalic locomotor region, locus ceruleus, and reticular formation). The lack of effects of the WEARHAP-PD device in the two patients in which vascular basal ganglia lesions were detected would indirectly support this hypothesis, as the burden of white matter subcortical hyperintensities correlates with the development of FOG and its resistance to treatment [START_REF] Chung | White matter hyperintensities as a predictor of freezing of gate in Parkinon's disease[END_REF]. Animal studies are required to definitely clarify neurophysiologically the relative weight of propriospinal or supraspinal mechanisms engaged by the WEARHAP-PD device alternating stimulation.

Whatever the neurophysiological mechanisms, current behavioral results are encouraging to look at the WEARHAP-PD device as a cheap, safe, and easy to use wearable device for the improvement of some walking parameters in PD patients under regular pharmacological treatment. Moreover, if direct effects on CPGs would really mechanistically prevail, there could be room to look at the use of WEARHAP-PD devices even in patients with some degree of cognitive impairment, a common condition in PD patients with FOG that may reduce the success of external cueing based on stimulation of specific sensory channels as hearing or sight [START_REF] Nonnekes | Comopensation strategies for gait impairment in Parkinson diseae[END_REF].

At a descriptive level, two patients in the current sample (number 8 and 10 of figures) showed vascular lesions of basal ganglia, poor response to levodopa, and highest scores to UPDRS-III items addressing bradykinesia and rigidity especially at lower limbs, as well as postural instability. Also in the other two patients who did not respond to the haptic cueing of the WEARHAP-PD device, bradykinesia and rigidity prevailed. These preliminary observations should be verified on larger samples that will help to verify which kind of Parkinsonian phenotypes are better responsive to the device.

Limitations of the study

While the small sample size of the study fits a proof-ofprinciple design addressing feasibility of the approach, safety, and acceptability of the device, it will need to be largely increased for a clinical trial. Such trial will benefit from a multicentric design, with gate and FOG measures taken outside the boundaries of the lab in a fully ecological domestic context, and possibly capitalizing by wearable sensors documenting quantitatively their occurrence.

The clinical trial should also address whether effects of the WEARHAP-PD device will change according to off and on states of patients and which patient's phenotype is the best candidate for using the device. For example, it is known that both age and PD can lead to a decline in haptic acuity [START_REF] Conte | Pathophysiology of somatosensory abnormalities in Parkinson disease[END_REF][START_REF] Konczak | Parkinson's disease accelerates age-related decline in haptic perception by altering somatosensory integration[END_REF]; therefore, a longer disease duration in elderly patients might prevent the optimal functioning of the device.

Another limitation that generalizes to all other neuroimaging investigations in PD might be represented by the difficulty to causally demonstrate in humans the hypothesized effect on spinal CPGs.

Fig. 2 a

 2 Fig. 2 a Available workspace (12.20 × 5.70 m) with cameras and trackable workspace (dashed line) of 10 × 4 meters. Tracker software workspace shows the Vicon cameras and the patient walking from point A toward point B. b The subject is equipped with 39 retro reflective

Fig. 1

 1 Fig.1The WEARHAP-PD device. The anklets are composed of two vibrating motors (1) attached to an elastic band (3). A Li-Ion battery and an Arduino board are in[START_REF] Nonnekes | Freezing of gait: a practical approach to management[END_REF]. When the device is engaged, its

Fig. 3

 3 Fig. 3 In all panels, white and black bars identify results with the anklets switched off and on, respectively. Numbers on the x axis identify single patients. a Mean stride velocity (m/s). b Stride standard deviation

Table 1

 1 Demographic and clinical characteristics of the 15 PD patients

	Age	61 ± 7.3 years
	Disease duration	9 ± 4.9 years
	LEDD	Mean, 530 ± 180; range, 350-900 mg/die
	UPDRS-I	1.7 ± 1.7
	UPDRS-II	14 ± 7.3
	UPDRS-III	26 ± 13.2
	UPDRS-IV	1.8 ± 1.2 (no troubling dyskinesia)
	FOG questionnaire	11 ± 6.1
	Tinetti (balance)	11 ± 3.4
	Tinetti (gate)	8 ± 2.8

LEDD levodopa equivalent daily dose UPDRS Unified Parkinson Disease Rating Scale FOG freezing of gate

Funding information The research leading to these results has received funding from the European Union's Horizon 2020 research and innovation program-Societal Challenge 1 (DG CONNECT/H) under grant agreement n. 643644 of the project "ACANTO: A CyberphysicAl social NeTwOrk using robot friends" and under grant agreement n°601165 of the project ``WEARHAP-WEARable HAPtics for humans and robots".

Compliance with ethical standards

Conflict of interest None.

Ethical approval None.