Bertoni Leaves

Jean-Baptiste Jentzer
email: jeanbaptistejentzer@ensiacetfrq.-8.jentzer

Marion Alignan
email: marion.alignan@ensiacet.fr

Carlos Vaca-Garcia
email: carlos.vacagarcia@ensiacet.fr

Luc Rigal
email: luc.rigal@ensiacet.fr

Gérard Vilarem
email: gerard.vilarem@ensiacet.fr

Response surface methodology to optimise Accelerated Solvent Extraction of steviol glycosides from Stevia rebaudiana

Keywords:

Following the approval of steviol glycosides as a food additive in Europe in December 2011, large-scale stevia cultivation will have to be developed within the EU. Thus there is a need to increase the efficiency of stevia evaluation through germplasm enhancement and agronomie improvement programs. To address the need for faster and reproducible sample throughput, conditions for automated extraction of dried stevia leaves using Accelerated Solvent Extraction were optimised. A response surface method ology was used to investigate the influence of three factors: extraction temperature, static time and cycle number on the stevioside and rebaudioside A extraction yields.

Introduction

Stevia rebaudiana Bertoni is a perennial herb of the Asteraceae family, native to the Amambay and Iguaçu districts on the borders of Brazil and Paraguay (I<inghom & Soejarto, 1985). Its principal interest is the high content of natural, dietary valuable sweeteners in its leaves (about 4-20% in leaves, dry weight) (Lemus-Mondaca, Vega-Galvez, Zura-Bravo, & Ah-Hen, 2012). These sweet compo nents are ent-kaurene diterpenoid glycosides, and today more than 30 are recognised in the scientific literature (Wèilwer-Rieck, 2012). The main compound is stevioside (4-13% of leaves dry weight) with a sweetening power of 250-300 (the sweetening power of sucrase, 'table' sugar, is 1 ), and the next most abundant is rebau dioside A (2-4%) with a sweetening power of 300-450. The other well known steviol glycosides (SG) are dulcoside A (0.3%), rebau dioside B, C (1-2%), D and F and steviolbioside [START_REF] Geuns | Stevia and steviol glycosides[END_REF][START_REF] Tavarini | Stevia rebaudiana Bertoni as a source of bioactive compounds: The effect ofharvest time, experimental site and crop age on steviol glycoside content and antioxidant properties[END_REF]. As reported in the literature, rebaudio side Ais recognised as Jess astringent than stevioside and therefore preferred by the food and beverages industry (I<ennelly, 2002;[START_REF] Tanaka | Steviol-glycosides: New natural sweeteners[END_REF]. The stevia sweeteners have functional and sensory properties superior to those of many other high-potency sweeten ers, like aspartame or cyclamate (Goya!, Samsher, & Goya!, 2010). Many studies suggest that they have antioxidant, anti-diabetic, anti-diarrheal, anti-hyperglycemic, anti-hypertensive, anti-inflam matory, anti-tumour, diuretic and immunomodulatory effects [START_REF] Chan | A double-blind placebo-controlled study of the effectiveness and tolerability of oral stevioside in human hypertension[END_REF]Gregersen, Jeppesen, Holst, & Hermansen, 2004;[START_REF] Puri | Optimisation of navel method for the extraction of steviosides from Stevia rebaudiana leaves[END_REF]Shivanna, Naika, I<hanum, & I<aul, 2013).

The Guarani Indians have long been using the leaves of stevia as a sweetening agent. The first crops outside Paraguay were grown by Sumida in Japan in 1968 [START_REF] Brandie | Stevia rebaudiana: Its agricultural, biological, and chemical properties[END_REF]. After, stevia appeared in many countries, such as Brazil, China, Indonesia, Israel, Tanzania and the United States [START_REF] Crammer | Progress in the chemistry and properties of rebaudiosides[END_REF][START_REF] Yoda | Supercritical fluid extraction from Stevia rebaudiana Bertoni using C02 and C02 + water: Extraction kinetics and identification of extracted components[END_REF]. In many of these countries the dried stevia leaves and their extracts were approved as food additives. In the United States the FDA (Food and Drug Administration) banned stevia until 1995 but highly purified steviol glycosides finally received GRAS (Generally Recog nised As Safe) status in both 2008and 2009(Wèilwer-Rieck, 2012). Introduction of steviol glycosides onto the European market as food additives was approved at the end of 2011 (European Commission, 2011 ).

Currently, the main stevia producing areas are China, especially in the north, and South Asia [START_REF] Tavarini | Stevia rebaudiana Bertoni as a source of bioactive compounds: The effect ofharvest time, experimental site and crop age on steviol glycoside content and antioxidant properties[END_REF]. Since a growing market for stevia products has developed in Europe, par ticularly following SG approval as a food additive within the EU on 2 December 2011, it would be of great interest to set up large-scale fully traceable crops to ensure the supply of stevia leaves in this area. To implement such crops, screening of different genotypes adapted to European soi! and climatic conditions is required. On a laboratory scale, SG have up to now been mostly extracted in a batch system by conventional methods: solvent extraction [START_REF] Bondarev | Peculiarities of diterpenoid steviol glycoside production in in vitro cultures of Stevia rebaudiana Bertoni[END_REF][START_REF] Liu | Subcritical fluid extraction of stevia sweeteners from Stevia rebaudiana[END_REF][START_REF] Woelwer-Rieck | Improved HPLC method for the evaluation of the major steviol glycosides in leaves of Stevia rebaudiana[END_REF], Soxhlet extrac tion [START_REF] Erkucuk | Supercritical C0 2 extraction of glycosides from Stevia rebaudiana leaves: Identification and optimization[END_REF], extraction by heating under reflux (Teo, Tan, Yang, Hew, & Ong, 2009), and cold extrac tion [START_REF] Erkucuk | Supercritical C0 2 extraction of glycosides from Stevia rebaudiana leaves: Identification and optimization[END_REF][START_REF] Jaitak | An efficient microwave-assisted extraction process of stevioside and rebaudioside-A from Stevia rebaudiana (Bertoni)[END_REF][START_REF] Teo | Validation ofgreen solvent extraction combined with chromatographie chemical fingerprint to evaluate quality of Stevia rebaudiana Bertoni[END_REF]. These methods suffer from being rather long (at least 1 h) with other drawbacks, like complicated thermoregulation without specific equipment, or use of an environmentally unfriendly solvent ( e.g. methanol). Moreover, little data concerning the extraction efficiency or reproducibility are given in the literature. Although other extraction methods are referenced, e.g. ultra sound-assisted extraction (UAE), microwave-assisted extraction (MAE) or supercritical fluid extraction (SFE) Uaitak et al., 2009;[START_REF] Liu | Subcritical fluid extraction of stevia sweeteners from Stevia rebaudiana[END_REF][START_REF] Teo | Validation ofgreen solvent extraction combined with chromatographie chemical fingerprint to evaluate quality of Stevia rebaudiana Bertoni[END_REF], these processes are also time con suming and not adapted to analysing a large number of samples.

Accelerated Solvent Extraction (ASE) is an automated extraction technique using elevated temperatures and pressures to achieve extraction in very short periods of time. The high temperature leads to better analyte solubility, faster diffusion rates, lower sol vent viscosities and weakening of the solute-matrix interactions. The high pressure allows for working with solvents as liquids above their boiling point and accelerates the overall extraction process [START_REF] Mustafa | Pressurized liquid extraction as a green approach in food and herbai plants extraction: A review[END_REF]. Finally, compared to conven tional extraction methods, higher automation, extraction yields, recovery and a shorter duration can be achieved [START_REF] Kukula-Koch | Influence of extraction procedures on phenolic content and antioxidant activity of Cretan barberry herb[END_REF]. This method can also be used to increase the through put of samples in stevia germplasm evaluation screening pro grams. Therefore ASE was the chosen method for this study. Conceming the parameters that can influence extraction yield, temperature, static time and number of cycles are the main ones and will be studied here using response surface methodology (RSM). This approach allows a reduction in the number of experi ments, and a prediction of the influence of the factors on the chosen response using a mathematical mode!. The latter can be graphically represented with response surfaces that show to what extent the influence of parameters or their interactions are signif icant or not, and can then provide the optimal conditions to improve a process [START_REF] Bharathi | Optimization of extraction of phenolic compounds from Avicennia marina (Forssk.) Vierh using response surface methodology[END_REF][START_REF] Rai | Optimization of process parameters for water extraction of stevioside using response surface methodology[END_REF]. Additionally, RSM was used on the ASE of stevia leaves to understand the impact on SG extraction yield of the main parameters ( extraction temperature, static time and number of cycles) and establish preliminary ASE conditions. An additional optimisation and experimental validation were then undertaken to set up the conditions for the ASE method retained.

Materials and methods

Plant material

For extraction method optimisation, S. rebaudiana Bertoni leaves, variety Criolla, were obtained from Paraguay. Leaves of S. rebaudiana Bert. genotypes 1-5 were obtained from Salvagnac, South West France, genotype 6 from Riom, Central France and genotype 7 from Nicaragua. Plant material was sun dried and stored at room temperature. Before extraction, dried leaves were ground in an electric grinder P19 Pulverisette (Fritsch, Germany) to pass through a 1 or 0.5 mm diameter circular mesh sieve.

Chemicals

Acetonitrile (analytical grade) and formic acid were provided by Sigma-Aldrich (France). Deionised water (resistivity > 18.2 MO cm) was made in situ using a Milli-Q Gradient AlO system purchased from Millipore Corporation (Bedford, MA, USA). Stevioside (purity ~95%) and rebaudioside A (purity ~98%), used as standards in this study, were purchased from ChromaDex (Molsheim, France).

Determination of moisture content

To determine the dry matter content for each sample, 1 g of leaves was oven-dried at 103 ac in an Air Concept AC 120 system (Firlabo, France) to constant weight.

Extraction methods

Conventional extraction method

Conventional extraction was performed with 15 g of stevia leaves, by adding 600 ml of deionised water in a 1 1 bottle and maintaining at 60 ac for 120 min with magnetic stirring at 250 rpm. Extracts were diluted if necessary and the amount of ste vioside and rebaudioside A were determined by HPLC, method described below, to find the extraction yield of SG from extracts.

Pressurised solvent extraction procedures

An ASE 350 system Dionex Corporation (Sunnyvale, CA, USA) was used for the pressurised liquid extraction. The cells were equipped with a stainless steel frit and a cellulose filter at the bot tom to avoid collection of suspended particles in the collection via!. Five hundred milligrams of crushed leaves were mixed with Fon tainebleau sand (VWR, France) and loaded into a 10 ml stainless steel extraction ce!!. Following the addition of purified water, the cell was pressurised, heated and extracted statically under condi tions obtained from RSM (Table 1) with a rinse volume (100%) of 5 ml and a 100 s nitrogen purge at the end of each extraction. The pressure in the cell was set at 10.3 MPa. The aqueous extract was completed to exactly 50 ml in a volumetric flask with deion ized water. Extracts were diluted if necessary and the amount of stevioside and rebaudioside A were determined by HPLC, method described below, to determine the extraction yield of SG in extracts.

Determination of maximal extractable mass of steviol glycosides in stevia leaves

In order to determine the maximal extractable SG content in stevia Ieaves, three successive Accelerated Solvent Extractions (sol vent = water, temperature = 60 °C, static time = 15 min and 1 cycle) were performed on the same cell to exhaust the Ieaves. The exper iment was made in triplicate and each fraction analysed by HPLC.

Experimental design

This study was made using an experimental design based on a three-factor Doehlert's type uniform network as described else where [START_REF] Doehlert | Uniform shell designs[END_REF]. This consisted of a three-variable (k = 3) Doehlert's design and needed thirteen experiments (plus five extra repetitions at the centre of the experimental domain). In the pres ent study, the ranges of experimental parameters were selected, based on Iiterature and extraction kinetics' results presented below. The experimental conditions are shown in Table 1.

The response values Ywere expressed as the main SG extraction yield obtained by ASE, i.e. stevioside and rebaudioside A. y = mass of extracted stevioside and rebaudioside A in leaves maximal extractabie mass of stevioside and rebaudioside A in leaves A full second-order polynomial mode) of the design was used to evaluate the extraction yield as the response (Y), as a function of the three independent variables (xi), namely extraction tempera ture (x 1 ), static time (x 2 ) and number of cycles (x 3 ), and their interactions.

The behaviour of the system can be described by the following second-order polynomial equation:

k k� k k Y= ao + L)ixi + L L aijxixi + LaiiXf i=l i=l j=i+l i=l
where Y is the response value, a 0 the response value at the centre of the experimental domain, ai are the linear coefficients, aii are the quadratic coefficients, a i i are the interaction coefficients, and x i the coded values of the three independent variables. 6%

7. HPLC analysis of extracts

Ali the extracts were diluted with deionized water and filtered through a 0.45 µm cellulose membrane filter (Xilab, Belgium) before being subjected to HPLC analysis. The HPLC system con sisted of a Dionex P680 pump, a Dionex ASI-100 automated sample injector, a Dionex Ultimate 3000 Column Compartment (Thermo Scientific, France), a diode array detector (DAD) (Hewlett Packard, France) and piloted by a Chromelon chromatography data system (Thermo Scientific, France) to obtain chromatographie profiles of the extracts. The chromatographie method was adapted for the JECFA method of SG assay published in FAO JECFA Monographs 10 (2010). The mobile phase consisted of0.1% aqueous formic acid (solvent A) and acetonitrile (solvent B). The elution mode was iso cratic at 1 ml/min for 30 min; 69-31% solvent B. Samples were injected (20 µI) onto an apolar Luna reversed-phase C18 column (250 x 4.6 mm ID, 5 µ, Phenomenex, Le Pecq, France) at 40 °C. Chromatograms were recorded at 200 nm. The standards used in the experiments were weighed and dissolved in deionized water. Calibration curves were generated with concentrations ranging from 50 to 500 mg/1 of stevioside and rebaudioside A.

Results and discussion

Selection of experimental conditions and ranges

In general, the extraction efficiency of a compound is influenced by multiple parameters such as time, temperature and solvent polarity, and their effects may be either independent or interactive. Before developing the study using RSM, extraction kinetics were determined by conventional extraction methods to select experi mental conditions and ranges. Steviol glycosides' extraction kinet ics were undertaken in triplicate in a batch system, at a temperature of 60 °C and with water as a solvent [START_REF] Giovanetto | Method for the recovery of steviosides from plant raw material[END_REF]. In order to measure the influence of particle size, Ieaves were used either whole or ground to pass through a 1 mm sieve. Results were plotted Fig. 1. Firstly, as can be seen, the kinetics of the two solutes are the same, so the total yield of stevio side + rebaudioside A can be chosen as the response for the model. Secondly, for ground Ieaves the speed of the extraction was faster and levelled off after 20 min, meaning that an equilibrium had been reached and that the extraction has ended. From this it was concluded that 20 min were sufficient to obtain the maximal SG concentration with ground Ieaves under these conditions. Thirdly, and finally, as stevioside is stable in solution up to 100 °c (Kinghom, 2002), this temperature value was chosen as the upper limit for our experiments.

� l: .!!!' 5% •• I ••••••••• •••• I ••••••••••• •• :C• •••••••• ••• • I• . ••••••• .... 3: • •••••••••••• ;i •• •• • �--..i 4% 3% � ..... ::i 2% 1% 0% o I• •• --r--- _ _.
.,,,,,, r -

Maximal extractable mass of Steviol glycosides in stevia leaves

Results from exhausting the Ieaves showed that about 86% of the SG were removed with the first extraction and that there was no more SG quantifiable (Limit Of Quantification < 0.01 g/1 of stevia extract, as defined by the European Pharmacopoeia ( 2007)) in the third extract (data not shown). After the second extraction, the rebaudioside A total content was 3.16 (±0.08)g/100 g of Ieaves dry weight and stevioside was 6.32 (±0.22)g/100 g of leaves dry weight, giving a standard profile of 67% stevioside and 33% rebau dioside A in agreement with [START_REF] Lemus-Mondaca | Stevia rebaudiana Bertoni, source of a high-potency natural sweetener: A comprehensive review on the biochemical, nutritional and functional aspects[END_REF] and [START_REF] Tavarini | Stevia rebaudiana Bertoni as a source of bioactive compounds: The effect ofharvest time, experimental site and crop age on steviol glycoside content and antioxidant properties[END_REF]. An example of a first extract chro matogram is given in Fig. S1 (see Supplemental material).

A mathematical mode! to describe the processing of ASE of steviol glycosides from S. rebaudiana Bertoni leaves

The Doehlert experimental design for three variables, the thir teen experiments and the experimental yields of the extractions are shown in Y= 93.78(±0.62) + 6.71 (±0.76)x 1 + 5.58(±0.76)x2 + 5.57(±0.76)X3 -1.02(±1.76)X1X2 -4.39(±1.97)X1X3 + 4.45(± 1.97)x 2 x 3 -2.28(± 1.24)x� -3.92(± 1.24)x� -2.87(± 1.17)x� For a detailed calculation of the coefficients of the quadratic equation, see [START_REF] Rossi | Application of the Doehlert experimental design to molecularly imprinted polymers: Surface response optimization of specific template recognition as a function of the type and degree of cross-linking[END_REF]. An F-test and a Student's T-test were used to check the statistical significance and suitability of the model. The F value (0.95) was Jess than the Fisher F parameter at the 95% confidence level ( 5.41 ), so the mathematical mode) fits with the data set. [START_REF] Joglekar | Product excellence through design of experiments[END_REF] have suggested that for a good fit of a mode), the regression coefficient R 2 should be at least 80%. Thus the proposed model had a sufficiently high correlation coefficient value (R 2 = 0.9642) to indicate that the extraction data was adequately explained, i.e. 96% of the data can be explained by the mode). Therefore, the mode) can be used as an estimate of a trend, to be able to predict the influence of the extraction parameters on the extraction yield. Fig. 2 shows that extraction temperature was the most influential parameter, and interaction between extraction temperature and static time, the Ieast. Moreover, only the linear terms were significant (p < 0.05), meaning that the response variables depended more upon the individual change of the independent variables, rather than their interactions.

Effects of extraction variables on extraction yield of steviol glycosides

The three-dimensional response surfaces shown in Fig. 3 for the independent variables (static time, temperature and number of cycles) were obtained by keeping the number of cycles constant, which indicated the changes in SG extraction yield under different ASE conditions. Fig. 3a shows the evolution of the SG extraction yield according to the extraction temperature and static time with the number of cycles kept constant at 1. As can be seen on the three-dimensional response surface, the extraction yield values increased principally with increasing extraction temperature. This means that tempera ture had a strong effect on the yield, whereas static time had no significant effect. The contact time between the solute and solvent was rather short, therefore the main step influencing the extrac tion efficiency is the solute dissolving in the solvent. When the temperature was increased, the viscosity of the water was reduced, thereby increasing its ability to wet the matrix and solubilise the solutes. There is also more energy to break analyte-matrix bonds, thus the diffusion of these analytes is facilitated (Teo, Tan, Yang, Hew, & Ong, 2010). This increase of SG extraction yield with tem perature has also been observed by [START_REF] Puri | Optimisation of navel method for the extraction of steviosides from Stevia rebaudiana leaves[END_REF] and [START_REF] Erkucuk | Supercritical C0 2 extraction of glycosides from Stevia rebaudiana leaves: Identification and optimization[END_REF].

Fig. 3b illustrates the evolution of the SG extraction yield according to extraction temperature and static time, with the num ber of cycles kept constant at 2. As can be seen from the graph, the extraction yield values increased with increasing extraction tem perature and static time. Thus temperature and static time have a significant effect, contact time between the two phases is signif icantly longer, and higher extraction yields are obtained. On the one hand, the addition of an extraction cycle allowed for renewal

Experimental validation

Three repetitions were made under the above conditions ( extrac tion temperature = 100 °C, static time = 4 min and 1 cycle) and the following extraction yields were obtained: 88.3%, 95.2% and 91.9% with an average of 91.8% ± 3.4%. The relative difference with the 95.3% expected was 3.7%. In addition, the verification experiments clearly demonstrated a good fit for the curve, and good reproducibil ity of results for an extraction using these parameters.

In order to further increase extraction yield while minimising extraction time, the particle size was reduced, thus increasing the exchange surface between the matrix and solute. The chosen conditions were as follows: static time 4 min, temperature 100 °c, one extraction cycle and leaves ground in order to pass through a 0.5 mm sieve. Extraction experiments were performed three times and the following extraction yields were obtained: 98.6%, 104.6% and 99.3% with an average of 100.8% ± 3.3%, thus these extraction conditions were retained. As stated earlier, the SG Accelerated Solvent Extraction method was optimised using leaves from S. rebaudiana Bertoni var. Criolla, which is usually called a variety but is in fact a population. This explains why yields slightly greater than 100% were obtained.

Table 2 (see Supplemental material) summarises the different stevioside extraction methods ( exclu ding conventional methods) from S. rebaudiana leaves reported in the literature (Pol, Ostra, I<arasek, Roth, Benesova, I<otlafikova & Caslavsky, 2007;[START_REF] Erkucuk | Supercritical C0 2 extraction of glycosides from Stevia rebaudiana leaves: Identification and optimization[END_REF][START_REF] Jaitak | An efficient microwave-assisted extraction process of stevioside and rebaudioside-A from Stevia rebaudiana (Bertoni)[END_REF][START_REF] Liu | Ultrasonically assisted extraction of total carbohydrates from Stevia rebaudiana Bertoni and identification of extracts[END_REF][START_REF] Teo | Pressurized hot water extraction (PHWE)[END_REF][START_REF] Puri | Optimisation of navel method for the extraction of steviosides from Stevia rebaudiana leaves[END_REF][START_REF] Rai | Optimization of process parameters for water extraction of stevioside using response surface methodology[END_REF]. The ASE method we have developed is faster than the other methods except for Microwave Assisted Extraction developed by [START_REF] Jaitak | An efficient microwave-assisted extraction process of stevioside and rebaudioside-A from Stevia rebaudiana (Bertoni)[END_REF] that lasted 1 min. However, this uses mainly methanol as a solvent, which is not eco-friendly. Moreover, ASE provides filtered extracts, a big time-saver with many samples, and in addition gives a maximal extraction yield for stevioside and rebaudioside A as developed in this study. Yield data is often absent in the literature.

Experimental validation on different genotypes

Optimisation of the SG extraction method was made using the same Stevia sample (var. Criolla). In order to confirm that the cho sen conditions were suitable for any Stevia genotype, seven of these with different stevioside and rebaudioside A profiles were tested. Two extractions were performed on the same cell to confirm that there was no SG quantifiable in the second extract. Experiments were made in triplicate and SG contents determined (Table 3, in Supplemental material).

As expected, ail the SG were removed in the first extraction (i.e. no more SG quantifiable in the second extract), and the conditions set up were suitable for determining the SG content in the leaves of any stevia genotype. Therefore this ASE method can be used and is well-adapted for stevia genotype screening.

Conclusion

During ASE of steviol glycosides from stevia leaves, tempera ture, static time and number of cycles significantly influenced the extraction yield, and response surface analysis did not demonstrate interactions between these independent variables. For the condi tions retained, the second-order polynomial equation predicted the extraction conditions for a yield of 95.3% at 100 °C, 4 min static time and 1 cycle, and experimental results gave an average extrac tion yield of 91.8 ± 3.4%. An additional optimisation was made by reducing the particle size to under 0.5 mm and the experimental results gave an extraction yield of 100.8 ± 3.3%. Finally, the chosen conditions were validated on leaves from seven stevia genotypes with different stevioside and rebaudioside A contents, and a maximal extraction yield was obtained for each. Therefore, the ASE method developed was very fast (13 min 30 s), very efficient, reproducible and can be used for screening S. rebaudiana geno types. Compared to conventional and non-conventional extraction methods it is faster, eco-friendly, more convenient and does not require as much energy.

Fig. 1 .

 1 Fig.1. SG extraction kinetics for stevia whole leaves and ground stevia leaves (particle size = 1 mm) at 60 'C with liquid/solid ratio = 20.

Fig. 2 .

 2 Fig. 2. Standardized Pareto chart showing the effect of different factor terms on SG extraction yield values. Bars exceeding the vertical line on the graph indicate that the corresponding factor tenns are significant (p < 0.05 ).

Fig. 3 .

 3 Fig. 3. Response surfaces plot for SG extraction yield showing the effect of extraction temperature and static time with number of cycles kept constant at 1 (a), 2 (b) and 3 (c).

Table 1

 1 Doehlert design of factors in real values for optimisation of process variables in SG content (stevioside + rebaudioside A) of Stevia rebaudiana Bertoni leaves and exper-

	imental results.
	Experiments Extraction	Static	Number Experimental
		temperature extraction	of cycles extraction yield of
		x1 ( °C)	time x2	X3	stevioside +
			(min)		rebaudioside A y
					{%)
		70	4	2	91.3
	2	100	4	2	99.0
	3	85	7	2	97.5
	4	55	7	2	91.8
	5	40	4	2	84.0
	6	55	1	2	82.1
	7	85		2	89.6
	8	85	5	3	99.3
	9	55	5	3	98.0
	10	70	2	3	88.9
	11	85	3		90.0
	12	55	3		80.9
	13	70	6		88.0
	14	70	4	2	93.3
	15	70	4	2	95.1
	16	70	4	2	93.5
	17	70	4	2	95.6
	18	70	4	2	94.0

Table 1 .

 1 The maximum yield of extractable stevio side + rebaudioside A (99.3%) was recorded in experiment No. 8, with parameters of 85 °c extraction temperature, 5 min static extraction time and 3 extraction cycles. The Iowest steviol glyco sides' yield (82.1%) was found with 55 °C extraction temperature, 1 min static extraction time and 2 extraction cycles. The data yielded the following second-order polynomial equation (results are given with standard deviations):

Acknowledgements

We would like to thank the Midi-Pyrénées Region and BPI France (French govemment) for financially supporting this work (IUF 1107019M), the Epi Salvagnacois (Salvagnac, France) for lead ing the project,Joël Perret and Mathieu Penin (Stevia Natura, Riom, France) for providing us with Stevia leaves and good advice, and the Agronomy Laboratory and Agrophysiology Laboratory of the Purpan Engineering School (Toulouse, France) for their collaboration.

Appendix A. Supplementary data

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.foodchem.2014. 06.078.


Nicola Accettura
email: nicola.accettura@laas.fr

Eric Alata

Pascal Berthou

Daniela Dragomirescu

Thierry Monteil

Addressing scalable, optimal and secure communications over LoRa networks: challenges and research directions

Keywords: LoRa, LPWAN, Scalability, QoS, Security

The Internet of Things (IoT) enables large scale deployments of very low power devices connected through wireless lossy links and able to interact with the surrounding environment (sensing and actuation). Two main challenges are then present: make them communicate; handle their energy consumption while respecting some cost constraints. Low Power Wide Area Networks (LPWANs) tackle these challenges by offering long-distance coverage while guaranteeing the use of a very little amount of energy for communications. Among many LPWAN technologies, Long Range (LoRa) networks provide a very promising but incomplete basis for satisfying the needs expressed by the applications running on low power devices. This paper describes the LoRa technology from the architectural point of view and points out those aspects that permit its seamless integration into the IoT. As major contribution, a focus on the current and future research on LoRa networks is provided by inspecting three facets: scalability, Quality of Service and security.

INTRODUCTION

The availability of miniaturized and cheap circuitry is still contributing to the growth of the Internet of Things (IoT). In that, the market of IoT networks has been recently extended with long-range technologies in the form of Low Power Wide Area Networks (LPWANs). Several technologies have been developed, including Long Range (LoRa) networks, SIGFOX, WI-SUN and Narrow Band IoT. In order to set up novel monitoring applications, e.g., targeting optimization of industrial processes [START_REF] Thubert | Low-Power Wide-Area Networks at the IETF[END_REF] , low power radiolocation, and geofencing 2 , it is mandatory to build interoperable networks, ready to be plugged and played through the Internet. In facts, such technologies are featured by some common aspects including: optimized radio modulation, star topology, tens of bytes long frames, a few frames transmitted mostly uplink per day with very low datarates and variable Maximum Transmission Units (MTUs). Among these technologies, LoRa has triggered a huge research interest given the inherent availability of different communication schemes, each mapping to a trade-off between device energy efficiency and network scalability.

In such a novel technological context, this contribution aims at surveying the recent research performed by the authors on LoRa networks. In doing that, it has been recognized that there are three main facets to be followed in future investigations, as also highlighted in the title of this paper. Indeed, the increasing amount of low power devices connected to the Internet requires communication protocols able to handle scalable deployments. From a Quality of Service (QoS) point of view, the consequent spring of novel applications, each having antagonist resource requirements, gives impulse to an anticipated study of optimal resource sharing mechanisms able to accommodate differentiated traffic patterns. Moreover, such networks will require robustness through reliable implementation of security aspects. In details, Sec. 2 pictures the LoRa protocol stack with a bottomup approach, paving the way to a complete communication architecture. Then, to browse the main research challenges in LoRa networks and to promote their seamless integration into the IoT, Sec. 3 addresses the description of the current and future investigations from the point of views of scalability, QoS and security. Eventually, Sec. 4 draws conclusions.

LORA COMMUNICATIONS

Currently, LoRa networks are mainly defined through some specifications released by the LoRa Alliance. At the same time, the Internet Engineering Task Force (IETF) is already working on the definition of primitives that will allow IPv6-enabled LPWANs, including LoRa networks. The layered protocol architecture defined by both the LoRa Alliance and the IETF is a foundation on the top of which a service layer can support interoperable Machine-To-Machine (M2M) communications. This concept is illustrated in Fig. 1 and detailed in the following subsections. . Any datarate has an associated maximum transmission range and a maximum time-on-air . Decreasing the datarate, both the time-on-air of frames and the transmission range increase [START_REF] Petäjäjärvi | On the coverage of LPWANs: range evaluation and channel attenuation model for LoRa technology[END_REF] . End-devices choose the highest possible datarate, depending on their distance from the gateway. For example, end-devices whose distance from the gateway belongs to ] , -1 ] choose as transmitting datarate. The surface affected by is then the annulus around the gateway, with the radius of the outer ring being -1 and the radius of the inner ring being . The area of the annulus increases as decreasing the datarate.

LoRa physical layer

LoRaWAN medium access

From an architectural point of view, a Long Range Wide Area Network (LoRaWAN) is handled by a central controller, namely a network server. Such a server coordinates several gateways by exploiting reliable communications on cabled or wireless technologies. At the same time, very low power constrained end-devices can communicate with the gateways over LoRa links. Generally speaking, the LoRaWAN specification [START_REF]LoRaWAN TM 1.1 Specification[END_REF] defines the Medium Access Control (MAC) layer of LoRa networks.

Whenever an information is available in its output buffer, an end-device immediately turns on its radio to transmit a linklayer frame over one among the available channels. If a duty cycle policy [START_REF]LoRaWAN TM 1.1 Regional Parameters[END_REF] applies, such a transmission is properly delayed to comply with the related limitation. Meanwhile, since gateways stay overhearing on all available channels, such a frame is received by all the gateways falling into the transmission range of the end-device. However, only a portion of such gateways do not observe collisions, thus correctly decoding the frame, extracting the included payload, and forwarding it to the server jointly with the perceived signal strength. Hence, multiple copies of the frame are received by the server, that in turn discards replicas. Then, the server selects the most appropriate gateway that will acknowledge the frame reception. The server piggybacks in the acknowledgment (ACK) the information it may wish to deliver toward the considered end-device. In this sense, the LoRaWAN scheme relies on gateways and server to handle the communication complexity, while easing the related implementation on end-devices. Indeed, the needed circuitry is not complex, thus making LoRaWAN end-devices cheap and easily available.

As a major need for LoRaWAN operators, it is very important to get reliability. This is achieved by properly bounding communication pitfalls, mainly represented by very frequent frame collisions due to the ALOHA-like nature of the access scheme. However, collisions can also be avoided through both a proper understanding of communication patterns and a robust design of network architectures. To offer a trade-off between power-efficiency and the availability of communication resources, there are three modes of operations for end-devices. Mode A is the main communication pattern, with bidirectional data exchanges asynchronously and asymmetrically started by end-devices. With Mode B, end-devices listen to beacons broadcast by gateways to get synchronization on a common understanding of time. Being synchronized, end-devices can timely open reception windows, to listen to incoming frames from the network. Hence, this mode enables synchronous, bidirectional and symmetric communications. A synchronized access also permits Time Division Multiple Access scheme, which can be exploited either with a Slotted ALOHA access or with collision-free resource scheduling. Although some energy is wasted for synchronization purposes, the network capacity is increased, with time-critical applications benefiting from such a scheme. Finally, with Mode C, any communication exchange is bidirectional and asynchronously started by either an end-device or a gateway with a symmetric pattern. This is achieved by having all LoRaWAN devices continuously overhearing on the radio when they have nothing to transmit. This mode of operation fits actuators, or mains-powered end-devices that could perform relaying.

Internet-connected LoRa networks

Given the inherent capability of addressing up to ∼ 3.4 ⋅ 10 38 entities, IPv6 is the key enabler protocol of IoT-aided automated houses, buildings, factories, and cities. However, wireless low power networks use data-link frames whose maximum length is not enough to contain a 1280-bytes long IPv6 MTU. On this regard, within the Internet Engineering Task Force (IETF) standardization body, the "IPv6 over Low power WPAN" (6LoWPAN) working group (WG) first and the "IPv6 over Networks of Resource-constrained Nodes" (6lo) WG later have already been dealing with the design of compression and fragmentation schemes able to let IPv6 packets fit 127 bytes long IEEE802.15.4 link-layer frames.

Herein, the spring of long range low power networks has triggered a renewed interest into adaptation layers similar to 6LoWPAN and able to guarantee RESTful CoAP-enabled communications and IPv6 connectivity [START_REF] Haa | Comparison of 6LoWPAN and LPWAN for the Internet of Things[END_REF] . Indeed, the wireless communication of such networks, by involving sporadic transmissions of very small frames at low datarates, puts an additional design constraints on IPv6 compression and fragmentation schemes. To this aim, the "IPv6 over Low Power Wide-Area Networks" (lpwan) IETF WG has been created within the IETF to target IPv6 connectivity on the top of some LPWAN technologies. In details, an informational document has been produced to describe the common characteristics and the actual needs of LPWAN technologies supporting IPv6. In addition, a specification has been being elaborated to enable the compression and fragmentation of CoAP messages encapsulated into IPv6 packet traveling over LPWANs. In details, a Static Context Compression Scheme is adopted to let CoAP/UDP/IPv6 packets be encapsulated in very small link-layer frames. A context defines how the packet header format is built. The context is static because it is pre-installed on end-devices and not learned through packet exchanges, since this in not compatible with LPWAN characteristics (e.g., duty cycle limitations or very sporadic traffic). Such a compression mechanism is independent of the specific LPWAN technology.

oneM2M-interoperable LoRa deployments

The characteristics and capacities of LoRa networks address some of the problems that are found in connected object systems such as smart cities and even Industry 4.0. Nevertheless, the manipulation of a single communication protocol in these systems remains marginal. In general, it must be possible to interact with several protocols and be able to deal with the specific features of each of them. The solution was to propose standards for IoT able to take into account this challenge either in a specific field (e.g., in industrial contexts through the Industrial Internet Consortium) or in a multi-domain vision with the oneM2M standard. Indeed, the latter defines a common M2M service platform which can be implemented as an horizontal solution interconnecting all network devices. oneM2M enables multiple communication protocols binding, reuse of existing remote devices management mechanisms, and interworking with existing legacy devices. A specific entity called Interworking Proxy Entity (IPE) allows to connect specific protocols of sensors like those connected on a specific LoRa network. This standard is based on a RESTful approach with open interfaces to enable developing services and applications independently of the underlying network, thus easing the deployment of vertical applications. The oneM2M standard defines: (i) a set of services, i.e., registration, discovery, security, group management, data management, device management, network service exposure, communication, etc.; (ii) a structured data model based on resources; (iii) a standard Application Programming Interface.

It is worth remarking that the oneM2M standard does not define by now how to integrate and take care of specific LoRa features (limited payload sizes, limited bandwidth, energy management, selection of the most fitting LoRaWAN mode).

ISSUES AND CHALLENGES

The presented overview on LoRa networks has been focused on the innovative features that make it the preferred technology supporting very low traffic and reliable monitoring applications. However, there are still some issues that need to be addressed to make LoRa tailored for huge IoT deployments supporting differentiated types of traffic while preserving secure communications.

Scalable deployments

LoRa deployments will connect a huge number of end-devices performing sensing operations. Thus, it is mandatory to study the capacity of such technology, in order to understand its scalability. The capacity of LoRaWAN has also been studied for single gateway deployments [START_REF] Mikhaylov | Analysis of Capacity and Scalability of the LoRa Low Power Wide Area Network Technology[END_REF] and city scaled ones [START_REF] Centenaro | Long-range communications in unlicensed bands: the rising stars in the IoT and smart city scenarios[END_REF] . Based on the consideration about datarates and transmission ranges presented above in Sec. 2.1, it comes out that city scaled LoRa deployments should prefer the use of the highest datarate. Indeed, assuming that the end-device density per area stays constant, for a datarate lower than the fastest one, (i) more end-devices contend the access to the radio, and (ii) transmitting the same amount of information takes more time (i.e., longer time-on-air). As a consequence, the collision rate is higher for lower datarates, thus making worth using just the fastest datarate in dense deployments of city scaled scenarios. Remarkably, it has been quickly recognized that the throughput of Pure ALOHA based LoRaWANs had not been investigated yet, to target multi-channel multi-gateway deployments [START_REF] Accettura | Simulating scalable Long Range Wide Area Networks for very low power monitoring applications[END_REF] . In such environment, the classical ALOHA probability models do not suffice in picturing the increased capacity due to gateway redundancy. Moreover, some form of low power tracking and geofencing can be enabled through signal multilateration 2 . The system model sprang out from the assumption that gateways are placed according to regular patterns. In addition, the gateway redundancy was achieved by insuring that any enddevice falls into the coverage range of at least three gateways. The throughput formulas found in those works [START_REF] Accettura | Low power radiolocation through long range wide area networks: A performance study[END_REF][START_REF] Accettura | Simulating scalable Long Range Wide Area Networks for very low power monitoring applications[END_REF] were designed to reflect the variation in the density of end-devices per unit area. Such probabilistic models were also validated through wide simulation campaigns. The simulation of realistic duty cycle limitations [START_REF]LoRaWAN TM 1.1 Regional Parameters[END_REF] showed also that the throughput trend can slightly differ when each device transmits very frequently, even though such discrepancy can be easily predicted. As matter of facts, such an investigation has been the first milestone in the study of scalable LoRaWANs, with several research directions open ahead.

Benefit of multiple gateways:

To increase the number of gateways able to capture frames from end-devices, it is sufficient to reduce the relative distance among them. On this regard, the network capacity must be assessed when varying the ratio between the density of gateways and the density of end-devices.

Impact of ACKs: It has to be understood the trend in the throughput when uplink frame transmissions are confirmed (i.e., ACKs are required to the aim of reliability). Indeed, duty-cycle limitation put constraints also on the downlink transmissions, thus dwarfing the capacity of gateways to acknowledge correct frame receptions.

From LoRaWAN to 5g and beyond: A question to be answered is about which features of LoRaWAN will be imported in the future IoT standard technologies enabled in 5+ generation of cellular systems.

LoRa-equipped low-orbit satellite networks: With the recent availability of low-orbit satellites equipped with LoRa chipsets, there have been some envisaged scenarios [START_REF] Qu | LEO Satellite Constellation for Internet of Things[END_REF] where some inaccessible areas could be easily monitored. In that, a thorough investigation should be done to the aim of understanding how many end devices can be managed over very long ranges and how to handle intermittent communications.

Quality of Service

LoRaWANs have been evaluated in several research works in terms of coverage, throughput offered to end-devices, network capacity and scalability [START_REF] Yi | A Study of LoRa: Long Range & Low Power Networks for the Internet of Things[END_REF][START_REF] Mikhaylov | Analysis of Capacity and Scalability of the LoRa Low Power Wide Area Network Technology[END_REF][START_REF] Bor | Do LoRa Low-Power Wide-Area Networks Scale?[END_REF] . Evaluation results show that LoRaWANs can offer good performances in terms of coverage with limited amount of exchanged traffic 8 . However, the scalability of the network can be limited in terms of exchanged data due to the joint utilization of ALOHA access scheme with duty cycle limitations 2 . In fact, in a dense network of relatively high traffic devices, applications experience high delays and low reliability [START_REF] Bor | Do LoRa Low-Power Wide-Area Networks Scale?[END_REF] . In addition, a high amount of confirmed end-device traffic (requesting ACK from the network) leads to severe degradation of network performances [START_REF] Pop | Does Bidirectional Traffic Do More Harm Than Good in LoRaWAN Based LPWA Networks?[END_REF] . Accordingly, modes B and C (enabling a symmetric communication), may even trigger additional uplink communication by pulling data from the end-devices. Hence, to allow LoRaWAN technology to broaden its scope to applications with strict latency and/or reliability requirements, it is essential to make changes to offer better quality of communications [START_REF] El Fehri | LoRa technology: MAC layer operations and Research issues[END_REF] . Several ideas have been pursued to optimize the capacity of the LoRaWAN network and achieve better performances. For example, the use of mode B enables the use of resource scheduling that can be handled in a centralized way, in order to reduce the collision rate and to provide time determinism. In this sense, the following research lines are envisaged to be addressed to achieve optimal traffic management on LoRa networks.

Scheduling LoRaWAN access: The LoRaWAN performance shift when employing mode B has to be thoroughly investigated. On this regard, scheduling algorithms for LoRaWANs have to be properly designed to account for multi-gateway deployments.

Coexistence of mode A and B: Mode A and B can coexist on the same deployments. It is of utmost importance finding a threshold mechanism able to let end-devices switch among the two modes of operation depending on the traffic conditions.

Defining new LoRaWAN modes: End-devices with higher requirements in terms of data exchange will need to use other communication patterns, possibly half-way in the middle between mode B and mode C. Specific channels can be assigned to devices operating with such mode. The channel access has to be designed to respect duty cycle limitation, while spreading the end-device transmissions over the time based on pseudo-random calculations and soft synchronization tools.

Listen Before Talk (LBT): Such a channel access should be investigated as replacement for the current policy used for attaining the management of the duty cycle limitation.

LoRaWAN gateways connected through geostationary satellite links: Satellite links represent a backup solution in case of malfunctioning of backbone infrastructure among gateways. They can also be used to replace cable links onto inaccessible areas. The performance evaluation of such networks is indeed finalized to combine resource scheduling in beacon-enabled LoRaWANs with the communication scheduling on geostationary satellite infrastructures.

Security

LoRa provides confidentiality with encryption based on multiple keys. A good level of confidentiality is achieved if keys are different among LoRa devices of the same type and if the AES implementation works properly. However, in a critical context, integrity and availability are also important properties to be addressed on LoRa devices [START_REF] Powell | MAFTIA (Malicious and Accidental-Fault Tolerance for Internet Applications)[END_REF] . Even though the LoRaWAN MAC layer is publicly available, the LoRa PHY is not completely open, with the only documents available being a patent and a specification [START_REF] Semtech | LoRa TM Modulation Basics[END_REF] . It is important to assess the security of the LoRa physical layer.

The protocol has been assessed in several studies. A work by Aras et al. [START_REF] Aras | Selective Jamming of LoRaWAN using Commodity Hardware[END_REF] proposed to selectively jam one device while leaving other devices unaffected. This weakness impacts availability. Two other vulnerabilities were found [START_REF] Aras | Exploring the Security Vulnerabilities of LoRa[END_REF] . The first one is related to the attacker's ability to open a device to retrieve the keys. The second one is linked to the reply of an already transmitted message. The reverse engineering of the LoRa PHY [START_REF] Frédéric | Security in the Internet of Things -a LoRa study[END_REF] has been done using an Universal Software Radio Peripheral (USRP) and GNU-radio software. This reverse engineering led to the implementation of LoRa modulation/demodulation as a GNU-radio block. This result allows to retrieve the content of messages that respects the format described in the MAC specification [START_REF]LoRaWAN TM 1.1 Specification[END_REF] , and to confirm the previous vulnerabilities. Based on these studies, two future directions have been identified.

End-to-end security: To deal with the identified security issues, devices should be designed with end-to-end security, from the PHY to the application layer. In particular, a poor reception quality must cause a software exception. Nevertheless, this strategy must be implemented with care, otherwise, an excess of corrupted messages may overwhelm the software with exceptions.

Identification of the emitter: Investigating the behavior of chips may be used to establish a "behavioral" signature and to identify the emitter. This solution is based on the assumption that a chip and its behavior cannot be cloned.

CONCLUSIONS

LoRa networks will play a great role in shaping the access infrastructure for the Internet of Things. Belonging to the wider class of LPWANs, they will ease the connection of billions of low power devices to the Internet, thus increasing the M2M traffic offered for novel IoT applications. However, there is still some ongoing debate on the real scope of such networks. Investigating the scalability is then necessary for an anticipated management of large LoRa network sizes. At the same time, different applications will need different treatments, calling for built-in mechanisms, like resource scheduling, able to insure an acceptable level of QoS. Also a thorough investigation on the security of such network is needed, given the sensitive data generated and relayed over these networks. This contribution has detailed these open issues and future research directions, by fitting them to an interoperable communication stack ready to be plugged and played through the IoT. As additional outcome, the paper has also given some pitch on future integration of LoRa networks coordinated by satellite links, for a totally reliable network infrastructure.


Ald Éric

Solym Mawaki Manou-Abi

A NOTE ON CONVEX ORDERING FOR STABLE STOCHASTIC INTEGRALS

Keywords: 2010 Mathematics Subject Classification. 60G52, 60H05, 60E15, 60G44 Stable process, stable stochastic integral, convex ordering, forwardbackward stochastic calculus

We establish a convex ordering between stochastic integrals driven by strictly α-stable processes with index α ∈ (1, 2). Our approach is based on the forward-backward stochastic calculus for martingales together with a suitable decomposition of stable stochastic integrals.

Introduction

It is by now well-known that theoretical development of stochastic orders provides elegant and practical tools for comparison of random phenomena, motivated by various applications in financial mathematics, risk management and stochastic networks. One of the most relevant stochastic orders is the convex order. Given X and Y two integrable random variables, we say that X is less than Y in the convex order, and we write X ≤ cx Y , if we have the inequality

E [φ(X)] ≤ E [φ(Y )] ,
for all convex functions φ such that the expectations exist. Such a stochastic order has some nice properties, among others the stability under mixtures and under convolution. In particular, if X and Y are two integrable random variables such that X ≤ cx Y , then they have the same mean and moreover,

E[(X -K) + ] ≤ E[(Y -K) + ], K ∈ R,
where F + denotes the positive part of a given random variable F . The latter inequality reveals to be convenient for financial applications (option pricing, hedging) and is actually equivalent to another stochastic order commonly denoted ≤ icx , the so-called increasing convex order, for which φ is additionally assumed to be non-decreasing. Another interesting consequence of the convex ordering is to facilitate the computation of distances between distributions such as the stop-loss distances (including the famous L 1 -Wasserstein distance) and the Zolotarev distances, see for instance [START_REF] Boutsikas | On the distance between convex-ordered random variables, with applications[END_REF]. Historically, the notion of convex comparison was first introduced by Hoeffding in [START_REF] Hoeffding | On the distribution of the number of successes in independent trials[END_REF][START_REF] Hoeffding | Probability Inequalities for sums of bounded random variables[END_REF] to obtain quantitative tail estimates. We refer to the monographs [START_REF] Müller | Comparison methods for stochastic models and risks[END_REF][START_REF] Shaked | Stochastic Orders[END_REF] for an extensive analysis of many stochastic orders and their potential applications, with precise credit and references for this large body of work. From a dynamical point of view, the convex ordering has been investigated in the context of Markov processes and semimartingales. See for instance the paper [START_REF] Karoui | Robustness of the Black and Scholes formula[END_REF] dealing with diffusions, with a special emphasis on the Black-Scholes model, and also [START_REF] Bergenthum | Comparison of option prices in semimartingale models[END_REF][START_REF] Bergenthum | Comparison of semimartingales and Lévy processes[END_REF] for option pricing involving general semimartingales and Lévy processes. In all these papers, the key point is the so-called propagation of convexity property. Other references such as [START_REF] Rüschendorf | On a comparison result for Markov processes[END_REF][START_REF] Rüschendorf | Comparison of Markov processes via infinitesimal generators[END_REF] use techniques based on the comparability of the infinitesimal generators of Markov processes together with the stochastic monotonicity property. Recently, another fruitful method was considered in [START_REF] Klein | Convex concentration inequalities and forward-backward stochastic calculus[END_REF], namely the so-called forward-backward stochastic calculus for martingales, which was followed by a series of papers based on this technology, cf. [START_REF] Arnaudon | Convex ordering for random vectors using predictable representation[END_REF][START_REF] Breton | Convex comparison inequalities for exponential jump-diffusion processes[END_REF][START_REF] Breton | Bounds on option prices in point process diffusion models[END_REF][START_REF] Breton | Convex comparison inequalities for non-Markovian stochastics integrals[END_REF][START_REF] Laquerrière | Deviation inequalities for exponential jump-diffusion processes[END_REF]. This approach reveals to be convenient in the non-Markovian framework, when dealing with Brownian stochastic integrals or stochastic integrals driven by a compensated Poisson random measure, and might be extended without any additional effort to the case of stochastic integrals driven by a Lévy process. However, the underlying assumptions required to obtain these convex comparison results enforce the driving Lévy process to have bounded jumps, excluding for instance the important class of stable processes. Following these observations, the aim of this note is to establish a convex ordering for the stochastic integral T 0 K t dZ t driven by a real strictly α-stable process (Z t ) t∈[0,T ] , where T > 0 is a finite time horizon. We restrict our attention to the case α ∈ (1, 2) since the convex ordering only makes sense for random variables with finite mean. Although this stable process has sample paths of infinite variation and unbounded jumps, the forward-backward stochastic calculus combined with a suitable representation of stable stochastic integrals allow us to show the convex comparison

T 0 K t dZ t ≤ cx T 0 k(t)dZ t ,
where k is some deterministic function bounding, in a particular sense, the predictable process

(K t ) t∈[0,T ] .
The paper is organized as follows. In Section 2, we recall some basic material on stable processes and stable stochastic integrals. Then we state our main contribution of the present work, Theorem 2.1, in which we derive the convex ordering mentioned above. Our result is a special instance of Theorem 2.2, which is more general but less tractable for practical issues since it compares in the convex order forward and backward stable stochastic integrals. Our last part, Section 3, is devoted to the proofs of these two results, on the basis of the forward-backward stochastic calculus introduced in [START_REF] Klein | Convex concentration inequalities and forward-backward stochastic calculus[END_REF] together with the representation of stable stochastic integrals developed by Kallenberg in [START_REF] Kallenberg | Some time change representations of stable integrals, via predictable transformations of local martingales[END_REF].

Preliminaries and main results

Consider on a probability space (Ω, A, P) a strictly α-stable process Z := (Z t ) t∈[0,T ] where T is a finite time horizon and α ∈ (1, 2). Here, strict-stability means the following self-similarity property

(Z kt ) t∈[0,T ] d = (k 1/α Z t ) t∈[0,T ] ,
where k > 0 and the equality d = is understood in the sense of finite dimensional distributions. The characteristic function is of the form

ϕ Zt (u) = exp t +∞ -∞ (e iuy -1 -iuy)ν(dy) , t ∈ [0, T ],
where ν stands for the Lévy measure defined on R \ {0} by

ν(dx) := dx |x| α+1 c + 1 {x>0} + c -1 {x<0} .
The parameters c + , c -above are non-negative with furthermore c + + c -> 0. In the case c + = c -> 0, the process is said to be symmetric. The centered process Z has sample paths of infinite variation and unbounded jumps, and it is a pure-jump martingale with respect to its natural filtration

F Z t := σ(Z s : s ∈ [0, t]), t ∈ [0, T ].
Denoting F Z this filtration, we assume that it satisfies the usual hypothesis, that is, completeness and right-continuity. The Lévy-Itô decomposition is given by

Z t = b R t + t 0 |x|≤R x (µ -σ)(ds, dx) + t 0 |x|>R x µ(ds, dx), t ∈ [0, T ],
where R is some arbitrary positive truncation level (classically chosen to be 1) and µ is a Poisson random measure on R + × R with intensity measure σ(dt, dx) = dt ⊗ ν(dx). Here b R is the drift parameter given by b

R := - |x|>R x ν(dx) = - (c + -c -) R 1-α α -1 .
Denote respectively by Z R-and Z R+ the two independent Lévy processes defined by the integrals above. The first one has a compactly supported Lévy measure and is a square-integrable martingale with infinitely many jumps bounded by R on each compact time interval, whereas the second one is an integrable compound Poisson process with jumps larger than R. Let K be a

F Z -predictable process belonging to L 2 (Ω × [0, T ]), i.e., T 0 E[K 2 t ] dt < +∞.
Then the stable integral X T := T 0 K s dZ s is well-defined as a stochastic integral with respect to the Lévy-Itô decomposition of Z, that is

X T = A R T + X R- T + X R+ T with A R T := b R T 0 K t dt, X R- T := T 0 K t dZ R- t and X R+ T := T 0 K t dZ R+ t ,
the second aforementioned process being a square-integrable martingale whereas the two other integrals are constructed in the Lebesgue-Stieltjes sense. In particular X T is integrable and therefore the process X given by

X t := E[X T |F Z t ] = t 0 K s dZ s , t ∈ [0, T ], is a F Z -martingale.
Let us turn to the question we are interested in, namely to establish a convex ordering between stable stochastic integrals. When dealing with a symmetric stable process Z and deterministic integrands, we can proceed directly as follows. Let

k, k ∈ L 2 ([0, T ]) be two functions such that |k(t)| ≤ | k(t)| for any t ∈ [0, T ]. Then the stable integrals X T := T 0 k(t)dZ t and XT := T 0 k(t)dZ t ,
have the same distribution as the random variables Z σ T and Z σT respectively, where Z is a copy of the process Z and the quantities σ T and σT are defined by

σ T := T 0 |k(t)| α dt and σT := T 0 | k(t)| α dt.
Denoting F Z the natural filtration of Z , Jensen's inequality entails

E [φ (X T )] = E[φ(Z σ T )] = E[φ(E[Z σT | F Z σ T ])] ≤ E[E[φ(Z σT ) | F Z σ T ]] = E[φ(Z σT )] = E[φ( XT )],
where φ is a convex function such that the expectations exist. In other words we obtain the convex ordering X T ≤ cx XT . For instance, such a result might be used to obtain a convex ordering between stable driven Ornstein-Uhlenbeck processes. Indeed consider the following random variables

X T := e -θT X 0 + T 0 e -θ(T -t) dZ t and Y T := e -κT X 0 + T 0 e -κ(T -t) dZ t ,
where Z is a symmetric α-stable process independent of the centered initial value X 0 , and θ, κ ∈ R are two parameters such that θ ≥ κ. If we denote

k(t) := e -θ(T -t) and k(t) := e -κ(T -t) , t ∈ [0, T ],
then we have k(t) ≤ k(t) for any t ∈ [0, T ] so that the previous convex ordering yields

T 0 e -θ(T -t) dZ t ≤ cx T 0 e -κ(T -t) dZ t .
Applying now the following classical property on the increasing convex order ≤ icx , cf. [START_REF] Müller | Comparison methods for stochastic models and risks[END_REF]:

if 0 ≤ a ≤ b and E[X 0 ] = 0 then aX 0 ≤ icx bX 0 ,
then with the choices of a := e -θT and b := e -κT , we obtain the increasing convex ordering X T ≤ icx Y T (after a slight modification of the argument emphasized to derive the stability by convolution of the order ≤ icx ). Finally, since

E[X T ] = E[Y T ], the order X T ≤ icx Y T thus implies the convex order X T ≤ cx Y T .
In the case of non-deterministic integrands, the story is more delicate, as it is already the case for Brownian stochastic integrals, cf. [START_REF] Karoui | Robustness of the Black and Scholes formula[END_REF] and also the recent papers [START_REF] Breton | Convex comparison inequalities for non-Markovian stochastics integrals[END_REF][START_REF] Laquerrière | Deviation inequalities for exponential jump-diffusion processes[END_REF] in which the authors overcome this difficulty by using the forward-backward stochastic calculus combined with the Malliavin calculus. In the stable setting we focus on in this note, we adapt the forward-backward stochastic calculus initially developed in [START_REF] Klein | Convex concentration inequalities and forward-backward stochastic calculus[END_REF] and use a convenient representation of stable stochastic integrals emphasized by Kallenberg in [START_REF] Kallenberg | Some time change representations of stable integrals, via predictable transformations of local martingales[END_REF] to obtain the main result of this note. Since it is a particular case of the more abstract Theorem 2.2 which will be stated in a moment, we postpone the proof to Section 3. Recall that the notation F + and F - stand respectively for the positive and negative parts of a given random variable or deterministic function F . Theorem 2.1. Let Z be a strictly stable process with index α ∈ (1, 2). Let K ∈ L 2 (Ω × [0, T ]) be a predictable process with respect to the natural filtration of Z and let k ∈ L 2 ([0, T ]). Denote for any t ∈ [0, T ] the rates

γ t := c 1/α + K + t + c 1/α -K - t and γt := c 1/α + k + (t) + c 1/α -k -(t), λ t := c 1/α + K - t + c 1/α -K + t and λt := c 1/α + k -(t) + c 1/α -k + (t)
, and assume that they satisfy the a.s. conditions:

γ t ≤ γt and λ t ≤ λt , t ∈ [0, T ].
Then the following convex ordering between stable stochastic integrals holds:

T 0 K t dZ t ≤ cx T 0 k(t)dZ t .
Certainly, if the two inequalities in the assumptions are reversed, then the convex ordering is also reversed, as expected. Moreover, note that our result might be applied, among others, in the following two classical situations:

• the process Z is symmetric and a.s. |K t | ≤ |k(t)| for any t ∈ [0, T ];

• both K and k are of the same constant sign and satisfy a.s.

K t ≤ k(t) for any t ∈ [0, T ].
For practical issues, we mention that the convex comparison result above is available for all convex functions φ whose growth at infinity is at most polynomial of degree p, with 1 ≤ p < α. For instance, letting the convex function φ(x) := |x| p where 1 ≤ p < α, we need to verify that

E T 0 K t dZ t p < +∞,
so that we will obtain, under the assumptions of Theorem 2.1, the following inequality:

E T 0 K t dZ t p ≤ E T 0 k(t)dZ t p .
Let us proceed as in [START_REF] Joulin | On maximal inequalities for stable stochastic integrals[END_REF] by using a convenient truncation level R in the Lévy-Itô decomposition. Using the notation introduced previously with X T := T 0 K t dZ t , we have for any fixed x > 0,

P (|X T | > x) ≤ P |A R- T | > x/2 + P |X R- T | > x/2 + P X R+ T = 0 . Denote T R 1 the first jump time of the Poisson process (µ ({y ∈ R : |y| > R} × [0, t])) t∈[0,T ] on the set {y ∈ R : |y| > R}. Since T R
1 is exponentially distributed with parameter ν ({y ∈ R : |y| > R}), we have

P(X R+ T = 0) ≤ P(T R 1 ≤ T ) ≤ T ν ({y ∈ R : |y| > R}) = (c + + c -) T α R α .
Now Chebyshev's inequality and the isometry formula for Poisson stochastic integrals entail

P |X R- T | > x/2 ≤ 4 x 2 E   T 0 |y|≤R K t y (µ -σ)(dt, dy) 2   = 4 x 2 |y|≤R y 2 ν(dy) T 0 E[K 2 t ] dt = 4 (c + + c -) R 2-α (2 -α) x 2 T 0 E[K 2 t ] dt,
whereas by Markov's inequality,

P |A R T | > x/2 ≤ 2 |b R | x T 0 E[|K t |] dt = 2 |c + -c -| R 1-α (α -1) x T 0 E[|K t |] dt.
Therefore choosing the truncation level R = x and rearranging the terms above entail the existence of some positive constant, say M , depending on all the parameters except x, such that P (|X T | > x) ≤ M x -α . Finally, we obtain by integrating

E [|X T | p ] = +∞ 0 p x p-1 P (|X T | > x) dx ≤ 1 + M +∞ 1 dx x α-p+1 , which is finite since p < α.
In order to state our second main contribution of the present note, let us recall the notion of backward processes. A backward strictly α-stable process Z is a strictly α-stable process running backwards in time. It shares the same properties as the (forward) process Z and in particular it is a backward martingale with respect to its own filtration. Indeed if F Z denotes the complete and left-continuous filtration F Z t := σ( Zs : t ≤ s ≤ T ), then the backward martingale property is satisfied, that is

E[ Zs |F Z t ] = Zt , 0 ≤ s ≤ t ≤ T. For instance the process Z defined by Zt := Z T -t , with F Z t := F Z T -t , t ∈ [0, T ]
, is a backward strictly α-stable process. For the sake of briefness, the (forward) processes considered in the sequel are assumed to be right-continuous with left limits, whereas the backward processes are naturally supposed to be left-continuous with right limits. Now let K be a F Z -predictable process belonging to L 2 (Ω × [0, T ]). Then the backward stable integral X0 := T 0 Ks d Zs is defined similarly to the forward case, the backward Itô differential d meaning that the limit in probability is taken over sums of the type

pn i=1 Ht n i ( Zt n i-1 -Zt n i ),
for refining subdivisions (t n i ) of [0, T ]. The random variable X0 is integrable and thus the backward process X given by

Xt := E[ X0 |F Z t ] = T t Ks d Zs , t ∈ [0, T ],
is a F Z -martingale. Now we are in position to state our second result on which Theorem 2.1 is based. In order to relieve the presentation, the proof is also postponed to Section 3. Theorem 2.2. Let Z, Z be two independent forward and backward strictly stable processes, respectively, both with the same index α ∈ (1, 2). Denote c + , c -and c+ , c-the rates of the respective Lévy measures ν and ν. Let K, K ∈ L 2 (Ω × [0, T ]) be two predictable processes with respect to the natural filtrations of Z and Z, respectively. Denote for any t ∈ [0, T ] the rates

γ t := c 1/α + K + t + c 1/α -K - t
and γt := c1/α

+ K+ t + c1/α - K- t , λ t := c 1/α + K - t + c 1/α -K + t and λt := c1/α + K- t + c1/α - K+
t , and assume that they satisfy the a.s. conditions:

γ t ≤ γt and λ t ≤ λt , t ∈ [0, T ].

Then we have the following convex ordering between forward and backward stable stochastic integrals:

T 0 K t dZ t ≤ cx T 0 Kt d Zt .
Our result remains in the spirit of the following well-known fact: if Z is an independent copy of Z and Z is the backward process Zt = Z T -t , t ∈ [0, T ], then the distribution of the random variable Z t + Zt does not depend on time t ∈ [0, T ] and thus has the same distribution as Z T and Z0 .

Before turning to the proofs, let us mention that Theorem 2.2 might be potentially extended. Indeed, up to some slight changes in the proof below, we are able to consider stable stochastic integrals for which a drift term is added, at the price of weakening the convex order ≤ cx by the increasing convex order ≤ icx . Under the assumptions and notation of Theorem 2.2, consider the drifted integrals

T 0 K t dZ t + T 0 L t dt and T 0 Kt d Zt + T 0 Lt dt,
where the processes L, L are respectively F Z -and F Z -adapted and satisfy the a.s. inequality L t ≤ Lt for any t ∈ [0, T ]. Then one needs in the forthcoming proof of Theorem 2.2 to assume that the convex functions of interest are also nondecreasing to preserve the order of the drift parts, so that we obtain the increasing convex ordering

T 0 K t dZ t + T 0 L t dt ≤ icx T 0 Kt d Zt + T 0 Lt dt.
In particular under appropriate assumptions, the strict stability of the processes Z, Z might be relaxed since any integrable stable process can be transformed into a strictly stable process by removing its expectation (thus corresponding to a drift term).

Proofs

The core of the present part is devoted to the proof of Theorem 2.2 (the proof of Theorem 2.1 will follow immediately). Before entering into the details, let us concentrate first our attention on two important tools that will be used in the remainder of the paper, namely the forward-backward stochastic calculus developed in [START_REF] Klein | Convex concentration inequalities and forward-backward stochastic calculus[END_REF] and the decomposition of stable stochastic integrals due to Kallenberg [START_REF] Kallenberg | Some time change representations of stable integrals, via predictable transformations of local martingales[END_REF]. In the sequel we work under the assumptions and notation of Theorem 2.2.

Define the enlarged forward and backward filtrations F := (F t ) t∈[0,T ] and F := ( Ft ) t∈[0,T ] as

F t := F Z t ∨ F Z 0 and Ft := F Z T ∨ F Z t , t ∈ [0, T ],
where we recall that F Z and F Z are the natural filtrations generated by the processes Z and Z, respectively. Then Z is still a F-martingale which is F-adapted whereas Z is a F-martingale which is F-adapted. In particular the processes X and X defined by

X t := t 0 K s dZ s and Xt := T t Ks d Zs , t ∈ [0, T ],
inherit the properties of Z and Z, respectively. Let ∆X u := X u -X u-be the jump of X at time u ∈ (0, T ] and define similarly the backward jump of X at time u ∈ [0, T ) by ∆ Xu := Xu -Xu+ . Following [START_REF] Klein | Convex concentration inequalities and forward-backward stochastic calculus[END_REF], we recall Itô's formula for forwardbackward martingales, applied in our context to the pure-jump stable stochastic integrals X, X: for any real-valued function f ∈ C 1 (R 2 ) and any 0

≤ s ≤ t ≤ T , f (X t , Xt ) = f (X s , Xs ) + t s + ∂f ∂x 1 (X u-, Xu ) dX u - t - s ∂f ∂x 2 (X u , Xu+ ) d Xu + s<u≤t f (X u-+ ∆X u , Xu ) -f (X u-, Xu ) -∆X u ∂f ∂x 1 (X u-, Xu ) - s≤u<t f (X u , Xu+ + ∆ Xu ) -f (X u , Xu+ ) -∆ Xu ∂f ∂x 2 (X u , Xu+ ) ,
where ∂f ∂x i is the partial derivative of the function f with respect to the i th coordinate, i ∈ {1, 2}. Now let us introduce the decomposition of stable stochastic integrals emphasized by Kallenberg in [START_REF] Kallenberg | Some time change representations of stable integrals, via predictable transformations of local martingales[END_REF]. Since the following argument can be easily generalized to the backward case, we will only focus our attention on the forward situation. First, note that the Lévy measure ν of the strictly α-stable process Z is of the form

ν = c + ν + + c -ν -,
where ν + and ν -are the measures on R + and R -, respectively, with Lebesgue density |x| -α-1 . Hence the self-similarity entails trivially the following representation

Z = c 1/α + Z -c 1/α -Z ,
where Z , Z are two independent strictly α-stable processes with Lévy measure ν + . In other words, this decomposition is obtained by separating the positive and negative jumps of Z. Therefore the process X can be rewritten as

X t = t 0 c 1/α + K s dZ s - t 0 c 1/α -K s dZ s = t 0 c 1/α + K + s dZ s + t 0 c 1/α -K - s dZ s - t 0 c 1/α + K - s dZ s - t 0 c 1/α -K + s dZ s .
Then Kallenberg's decomposition for the stable stochastic integral X stands as follows:

X t = t 0 γ s d Ẑ s - t 0 λ s d Ẑ s ,
where the processes Ẑ , Ẑ given by the stable stochastic integrals

Ẑ t = t 0 1 {Ks≥0} dZ s + t 0 1 {Ks<0} dZ s and Ẑ t = t 0 1 {Ks<0} dZ s + t 0 1 {Ks≥0} dZ s ,
are independent and have the same distribution as Z , Z .

As announced above, we are now in position to prove Theorem 2.2.

Proof of Theorem 2.2. Our aim is to establish the inequality

E φ T 0 K t dZ t ≤ E φ T 0 Kt d Zt ,
for any convex function φ such that the expectations exist. Since this inequality rewrites as E φ(X T + XT ) ≤ E φ(X 0 + X0 ) , the proof will be achieved once we have established that the function t ∈ [0, T ] → E[φ(X t + Xt )] is non-increasing. At the price of an unessential regularizing procedure, we can assume without loss of generality that the convex function φ is also C 1 and Lipschitz. Using then the Itô formula above with f (x, y) = φ(x + y), we have

φ(X t + Xt ) = φ(X s + Xs ) + t s + φ (X u-+ Xu ) dX u - t - s φ (X u + Xu+ ) d Xu +J X s,t -J X s,t ,
where J X s,t and J X s,t are defined by

J X s,t := s<u≤t φ(X u-+ ∆X u + Xu ) -φ(X u-+ Xu ) -∆X u φ (X u-+ Xu ) , J X s,t := s≤u<t φ(X u + Xu+ + ∆ Xu ) -φ(X u + Xu+ ) -∆ Xu φ (X u + Xu+ ) .
By Kallenberg's decomposition of the stable stochastic integral X, we can rewrite J X s,t in the following way:

J X s,t = s<u≤t φ(X u-+ γ u ∆ Ẑ u -λ u ∆ Ẑ u + Xu ) -φ(X u-+ Xu )) -(γ u ∆ Ẑ u -λ u ∆ Ẑ u )φ (X u-+ Xu ) .
Taking then expectation, we obtain

E J X s,t = E t s +∞ -∞ +∞ -∞ φ(X u + Xu + γ u x -λ u y) -φ(X u + Xu ) -(γ u x -λ u y) φ (X u + Xu ) σ Ẑ , Ẑ (du, dx, dy) ,
where σ Ẑ , Ẑ is the F-dual predictable projection of the bivariate process ( Ẑ , Ẑ ). Since its coordinates are independent and have the same distribution as Z , Z , the measure σ Ẑ , Ẑ is the tensor product of the Lebesgue measure on R + and the two-dimensional Lévy measure concentrated on the axis and with two marginals ν + . Therefore we get from the above computations, To achieve this short note, it remains to prove our main result Theorem 2.1 on the basis of Theorem 2.2.

Proof of Theorem 2.1. The proof is a straightforward consequence of Theorem 2.2 since under the conditions c+ = c + and c-= c -and when the process K is actually deterministic (equal to the given function k), the backward stable stochastic integral T 0 Kt d Zt has the same distribution as the forward integral T 0 Kt dZ t .

  u + Xu + γ u x) -φ(X u + Xu ) -γ u x φ (X u + Xu ) du ν + (dx) u + Xu -λ u y) -φ(X u + Xu ) + λ u y φ (X u + Xu ) du ν + (dy) u + Xu + a) -φ(X u + Xu ) -a φ (X u + Xu ) u + Xu + b) -φ(X u + Xu ) -b φ (X u + Xu ) λ α u du db |b| α+1 ,where in the last line we used two obvious changes of variables. Hence we obtainu + Xu + x) -φ(X u + Xu ) -x φ (X u + Xu ) ×(γ α u 1 {x>0} + λ α u 1 {x<0} ) du dx |x| α+1 .Certainly, a similar method applied in the backward case allows us to get Xu +X u + x) -φ( Xu + X u ) -x φ ( Xu + X u ) ×(γ α u 1 {x>0} + λα u 1 {x<0} ) du dx |x| α+1 .Finally taking expectation in the Itô formula above and plugging the previous calculations, we haveE φ(X t + Xt ) -E φ(X s + Xs ) u + Xu + x) -φ(X u + Xu ) -xφ (X u + Xu ) (γ α u -γα u ) u + Xu + x) -φ(X u + Xu ) -xφ (X u + Xu ) (λ α u -λα u ) du dx |x| α+1 ,which is non-positive by the convexity of φ and our assumptions. The proof of Theorem 2.2 is now complete.

Acknowledgments. The authors thank Yutao Ma for useful discussions. They are also grateful to the French Agency Campus France for financial support.