Cyril Raveschot

Benoit Cudennec

Barbara Deracinois

Marc Frémont

Maxime Vaeremans

Jamyan Dugersuren

Shirchin Demberel

Djamel Drider

Pascal Dhulster

François Coutte

Dr Christophe Flahaut
email: christophe.flahaut@univ-artois.fr

Proteolytic activity of Lactobacillus strains isolated from Mongolian traditional dairy products: A multiparametric analysis

Keywords: Cell envelope associated peptidases, fermented milk, Lactobacillus, peptides Angiotensin Converting Enzyme ACN: Acetonitrile CEP: Cell Envelope associated Ultra High Temperature

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Milk proteins are precursors of many different peptides, some of which display interesting bioactivities: peptides with antimicrobial, antioxidant, opioid, immunostimulating or calcium binding activities have already been described [START_REF] Nielsen | Milk bioactive peptide database: A comprehensive database of milk protein-derived bioactive peptides and novel visualization[END_REF][START_REF] Nongonierma | The scientific evidence for the role of milk protein-derived bioactive peptides in humans: A Review[END_REF][START_REF] Park | Bioactive Peptides in Milk and Dairy Products: A Review[END_REF]. Bioactive peptides can be generated from milk proteins by enzymatic hydrolysis during gastrointestinal digestion, or produced in vitro with the use of commercial enzymes. However, industrial production of peptides is limited by the lack of suitable large-scale technologies, and by the high cost of enzymes [START_REF] Wu | Simultaneous production of multifunctional peptides by pancreatic hydrolysis of bovine casein in an enzymatic membrane reactor via combinational chromatography[END_REF].

Milk fermentation by Lactobacillus bacteria represents an inexpensive way to generate bioactive peptides [START_REF] Giacometti | Peptidomics as a tool for characterizing bioactive milk peptides[END_REF]. In order to source all the amino acids required for their growth, Lactobacillus bacteria indeed hydrolyze milk proteins during fermentation, and a large variety of peptides is generated in this process. The enzymatic system responsible for milk protein degradation is composed of cell envelope associated peptidases (CEPs) and peptidases. The CEPs first hydrolyze proteins into peptides, which are then transported into the cell by a peptide transport system; intracellular peptidases finally degrade these peptides into smaller molecules and amino acids [START_REF] Savijoki | Proteolytic systems of lactic acid bacteria[END_REF]. Only one CEP (named PrtB) has been described for Lactobacillus delbrueckii subsp. bulgaricus (L. delbrueckii) [START_REF] Hou | Effect of culturing conditions on the expression of key enzymes in the proteolytic system of Lactobacillus bulgaricus[END_REF] whereas 4 CEPs (PrtH1, PrtH2, PrtH3 and PrtH4) can be expressed, in a strain dependent manner, by Lactobacillus helveticus (L. helveticus) strains [START_REF] Griffiths | Lactobacillus helveticus: The proteolytic system[END_REF]. The proteolytic activity of a given strain finally depends on its CEP and peptidase profiles, on gene variations, and on the conditions required for optimal enzyme activity [START_REF] Alcántara | Peptide and amino acid metabolism is controlled by an OmpR-family response regulator in Lactobacillus casei[END_REF].

L. helveticus LBK-16H is an example of strain that generates bioactive peptides upon milk fermentation: two tripeptides isoleucyl-prolyl-proline (IPP) and valyl-prolyl-proline (VPP) were identified, which have been shown to act as angiotensin converting enzyme (ACE) inhibitors [START_REF] Narva | Effects of Lactobacillus helveticus fermented milk on bone cells in vitro[END_REF]. Another example is the antioxidative peptide present in milk fermented by L. delbrueckii IFO13953 [START_REF] Korhonen | Bioactive peptides: Production and functionality[END_REF].

This study aimed at evaluating the peptide production capacity of Lactobacillus strains isolated in three different areas of Mongolia. Production of dairies through spontaneous milk fermentation has a long history in this country. Herders use different milks (mare, cow, yak, camel, sheep, goat) to produce a variety of fermented food such as yogurt, kefir, or airag (fermented horse milk). These products were reported to contain a wide diversity of LAB, especially Lactobacillus species [START_REF] Watanabe | Diversity of lactic acid bacteria and yeasts in Airag and Tarag, traditional fermented milk products of Mongolia[END_REF]Yu et al., 2011). Mongolian L. helveticus strains isolated from mare's fermented milk present a different CEP distribution, and generate different peptides from casein hydrolysis, when compared to European and North American strains [START_REF] Miyamoto | Distinctive proteolytic activity of cell envelope proteinase of Lactobacillus helveticus isolated from airag, a traditional Mongolian fermented mare's milk[END_REF]. This suggests that Lactobacillus species from different origins could hydrolyze milk proteins in different ways, and therefore generate peptides with original bioactivities. Actually, new peptides with ACE inhibitory activity have been identified in milk fermented by the Mongolian strain L. helveticus 130B4 (Schuangquan, Tsuda & Miyamoto, 2008).

The screening was carried on 170 Lactobacillus isolates, and involved a multiparametric analysis based on the correlation between i) fermentation characteristics (strain growth, acidification pattern), ii) CEP gene distribution and iii) peptide heterogeneity obtained after milk fermentation. The study was able to distinguish strains with specific proteolytic activities, generating distinct populations of peptides upon fermentation of milk.

Materials and methods

Sample collection

A total of 67 dairy samples were collected from three regions of Mongolia, corresponding to different ecosystems and environmental conditions: Gobi desert (16 samples), Northern Mongolia Hangai region (12 samples) and Altaï mountains (39 samples). Collected samples included mare (Equus caballus) fermented milk, camel (Camelus bactrianus) fermented milk, as well as cow (Bos taurus), goat (Capra hircus), yak (Bos grunniens) and hainakh (Bos taurus × Bos grunniens) yogurts. Samples were collected in sterile containers and kept refrigerated for a maximum of two days before analysis.

Isolation and identification of LAB

Samples were serially diluted in saline solution (NaCl 0.85% w/v), streaked on MRS plates, incubated for 48h at 37°C under anaerobic conditions [START_REF] De Man | A medium used for the cultivation of Lactobacilli[END_REF].

Morphologically distinct colonies were picked and streaked again until purification. Strains were examined with an optical microscope for Gram staining. Catalase assay was performed by application of 3 % (v/v) hydrogen peroxide solution on cultures. Bubbling indicates catalase activity. A total of 170 isolates were finally retained in the collection. Strains were stored in MRS media and glycerol (15% v/v final concentration) at -80°C.

Identification of strains was performed by 16S rDNA amplification and sequencing. Strain DNA was extracted from bacterial colonies using QIAamp DNA Mini Kit (Qiagen, Hilden, Germany). DNA amplification and sequencing were performed by GATC biotech (Konstanz, Germany). Species were identified by comparing the 16S sequence with bacterial sequences in Medline databases, a score higher than 99% was considered significant.

Screening for proteolytic activities on skim milk agar

Proteolytic activity was determined using skim milk agar medium. Strains were grown in MRS at 37°C for 48h, washed twice in the same medium and adjusted to a final optical density at 600 nm (OD600) of 1. Twenty microliters of the suspension were placed in skim milk agar (5% w/v) wells (diameter 4 mm) and inoculated plates were incubated at 37°C for 72h. Proteolysis resulted in the formation of a clear zone around the wells; diameter of this clear zone was measured for proteolytic activity quantification.

Fermentation of Ultra High Temperature (UHT) skim milk with selected strains

In order to assess the strains proteolytic activity, milk cultures were prepared in a media composed of 100% UHT skim milk (Cora, France). Strains were inoculated to 30 mL of UHT skim milk, to reach an initial OD600 of 0.3, and grown at 37°C for 48h without agitation.

Concomitantly, a non-inoculated UHT skim milk (NISM) was incubated in the same conditions and served as negative control. After incubation, the pH value of culture media was determined using a pH meter (Mettler Toledo, Greifensee, Swiss) and all milk cultures were kept at -20°C for further analysis. Experiments were performed in triplicate and the results were expressed as the mean plus standard error.

Bacterial growth determination

Bacterial concentration in skim milk cultures was determined as follows: a volume of culture was diluted 10 times in 2% ethylenediaminetetraacetic acid (EDTA), pH 12, and then centrifuged at 13,000 × g for 10 min (Abraham, De Antoni & Anon, 1993). This process was repeated and the pellet was resuspended in phosphate buffer solution (PBS) at pH 7. OD of the solution was measured at 600 nm.

Peptide extraction and purification

Proteins were precipitated from 4 mL of milk cultures by addition of a trichloroacetic acid solution (TCA, 10% w/v in water) to reach a final TCA concentration of 1%. After centrifugation at 10,000 × g for 10 min at room temperature, the peptide-containing supernatant was collected. Sugars, salt and TCA were removed from supernatant using a C18 solid phase extraction with Bond Elut C18 1000 mg microcolumns (Agilent Technologies, Santa Clara, USA). Elution of C18-retained peptides was performed with 2 mL of 80% v/v acetonitrile (ACN), 20% v/v water and 0.1% v/v trifluoroacetic acid (TFA).

Peptide quantification

After extraction and purification, the assessment of peptide concentration in milk culture was performed using the Folin-Ciocalteu reagent St Quentin Fallavier, France)) [START_REF] Mcsweeney | Chemical methods for the characterization of proteolysis in cheese during ripening[END_REF][START_REF] Lowry | Protein measurement whit the folin phenol reagent[END_REF]. Briefly, the complexation reaction was performed in a final volume of 800 µL composed of 200 µL of extracted peptides, 500 µL of sodium carbonate solution (500 mM) and 100 µL of FC reagent. After incubation at 37°C for 30 min in the dark, the OD750 was measured to evaluate the peptide amount. Standard curve was established with a peptide digest standard (Peptide digest assay standard, ThermoFisher Scientific, Waltham, USA).

Molecular weight profiles of peptides obtained from milk hydrolysates

Peptide content in milk hydrolysates was analyzed by size exclusion chromatography (SEC) based on a fast protein liquid chromatography (FPLC) using a Superdex Peptide 10/300 GL column (10 × 300-310 mm, 13 µm, GE Healthcare, Little Chalfont, UK) coupled to an AKTA Protein Purification System (GE Healthcare, Little Chalfont, UK). Twenty-five microliters of sample were loaded, eluted in isocratic conditions with 30% ACN, 0.1% TFA solution at a flow rate of 0.5 mL.min -1 for 60 minutes and monitored at 214 nm. Five substances were used as internal standards to obtain a calibration curve for calculating molecular weight distribution: albumin (60 kDa), cytochrome C (12,400 Da), aprotinin (6,500 Da), vitamin B12

(1,355 Da) and glutathione (307.3 Da). Peptide content in each sample was estimated by peak area integration. Two molecular weight ranges were distinguished, one corresponding to the high molecular weight proteins/peptides (HMWPs) with an apparent molecular weight higher than 1,700 Da, and one corresponding to the low molecular weight peptides (LMWPs) with an apparent molecular weight lower than 1,700 Da. Peptide content estimation for each range was expressed either as percentage of total profile area (%HMWPs and %LMWPs) or as percentage compared to non-inoculated milk sample (%HMWPs/Milk and %LMWPs/Milk).

Polymerase chain reaction (PCR)-based detection of CEP genes

Bacteria from milk cultures were harvested as described in the 2.4.1. section. Cellular DNA was isolated using a Wizard Genomic DNA Purification Kit (Promega, Madison, USA). Each strain was tested with five different primer sets corresponding to different CEP genes. The primers described by [START_REF] Hou | Effect of culturing conditions on the expression of key enzymes in the proteolytic system of Lactobacillus bulgaricus[END_REF] were used for detection of prtB gene in L. delbrueckii strains. The primers used for the detection of prtH1 and prtH2 genes of L. helveticus were from [START_REF] Genay | PrtH2, Not prtH, Is the ubiquitous cell wall proteinase gene in Lactobacillus helveticus[END_REF] and the primers for prtH3 and prtH4 (same species) were from [START_REF] Broadbent | Genetic diversity in proteolytic enzymes and amino acid metabolism among Lactobacillus helveticus strains[END_REF]. PCR products were detected with GelRed staining after electrophoresis on a 1% w/v agarose gel. An O'GeneRuler DNA ladder Mix (ThermoFisher Scientific, Waltham, USA) was used as molecular weight marker. The strain CNRZ32 CIRM-BIA 103 (provided by CIRM-BIA, INRA, Rennes, France) was used as control for L. helveticus strains.

Peptides identification by UPLC-ESI-qTOF-MS/MS

The peptides purified on Bond Elut C18 microcolumns (section 2.4.2.) were dried by centrifugal evaporation (miVac Centrifugal Vacuum Concentrators, Gene Vac, Ipswich, UK) for 2h at 40°C. Dried peptides were redissolved in 100 µL H2O, 0.1 % TFA and centrifuged for 10 min at 8,000 × g. Ten µL were chromatographed on an ACQUITY UPLC system (Waters, Manchester, UK) using a C18-Kinetex column (150 x 4.6 mm, 2.6µ 100Å, Phenomenex, Torrance, USA), a linear gradient of ACN containing 0.1% FA (5-15% ACN over 30 min, 15-30% ACN over 60 min, 30-50% ACN over 10 min and 50-95% ACN over 10 min) and a flow rate of 500 µL.min -1 . The HPLC eluent was directly electrosprayed from the column end at an applied voltage of 3 kV, using a desolvation gas (N2) flow of 600 L.h -1 , a nebulizer gas flow of 2.5 bar and desolvation temperature of 300°C, respectively. The chromatography device was coupled to SYNAPT-G2-Si mass spectrometer (Waters, Manchester, UK). MS analysis was performed in sensitivity, positive ions and data dependent analysis (DDA) modes and MS data were collected in the 100-2000 m/z range with a scan time of 0.2 s. A maximum of 15 precursor ions with an intensity threshold of 10,000 counts were selected for ion trap collision-induced dissociation (CID) fragmentation and subjected to a collision energy ramping from 8V to 9V for low mass and 40V to 90V for high mass. The MS/MS spectra were recorded on the 100 to 2,000 m/z range with a scan time of 0.1 s. Peaks were analyzed using Mass Lynx software (ver.4.1, Waters, Manchester, UK).

Database searches via PEAKS Studio 7.0 (Bioinformatics Solutions) were performed using the UniProt databases (May 15, 2017) restricted to Bos taurus organism. A mass tolerance of 35 ppm and 3 missing cleavage sites and an MS/MS tolerance of 0.2 Da were allowed.

Variable methionine oxidation was also considered. The relevance of protein and peptide identities was judged according to their identification generated by PEAKS Studio 7.0 (p<0.05) and a false discovery rate < to 1%.

β-casein heat map construction from MS analysis

Heat maps were constructed from β-casein peptides. The frequency of appearance of each amino acid was calculated from all identified peptides. Apparition frequency was established from in silico approach [START_REF] Lacroix | Evaluation of the potential of dietary proteins as precursors of dipeptidyl peptidase (DPP)-IV inhibitors by an in silico approach[END_REF] using Microsoft Excel software ( 2013version) with the following formula: A = a/N, where A is the apparition frequency of an amino acid, a the amino acid occurrence and N the number of peptide sequences. According to calculated frequency, a color was associated to each amino acid from white (no apparition frequency) to dark grey (high apparition frequency).

Statistical analysis and principal component analysis (PCA)

All statistical differences between strains were measured using one way ANOVA and Tukey's post-hoc test. Differences between means were considered significant when p value < 0.05. In order to analyze all measured variables at the same time, a first principal component analysis (PCA) was performed on data corresponding to results from proteolytic activity on milk agar, growth determination, pH determination and estimation of peptide contents. A second PCA was conducted to compare data from MS analysis (β-casein heat maps). PCA results were represented both as a variables map (correlation between studied variables) and as individuals map (similarity between individuals). All statistical analyses and PCA were performed on R software (FactorMineR package) (R core team, 2016, Vienna, Austria).

Results

Identification of Lactobacillus species in the collected dairy samples

One hundred and seventy (170) strains of LAB were isolated from 67 dairy products (fermented milks and yogurts) collected in three different regions of Mongolia: Gobi desert, Hangai region, Altaï mountains. Twelve different bacterial species were represented among the isolates (Figure 1A). Three dominant species accounted for more than 75% of the isolates: L. helveticus (73 strains), L. delbrueckii (39 strains) and L. kefiranofaciens (29 strains). There were however differences between the three regions: the dominant species in samples from Gobi area was L. kefiranofaciens (half of the isolates, retrieved in all types of samples including yogurt, airag, camel fermented milk) whereas this species was virtually absent in the Altaï samples (only one isolate). In the Altaï samples, the dominant species was L. helveticus, which for instance accounted for 76% of the isolates obtained from airag. In airag samples collected in the Gobi area, L. helveticus accounted for only 42% of all isolates.

Selection of highly proteolytic L. helveticus and L. delbrueckii strains

A standardized skim milk agar-well diffusion test was used to evaluate the proteolytic activity of all 170 strains. Strains showing the highest proteolytic activity (clear zone diameter larger than 28 mm, Table 1) belonged to the species L. helveticus and L. delbrueckii. Fifteen strains (nine L. helveticus, six L. delbrueckii) were selected for the further steps of the screening strategy presented in Figure 1B. Briefly, this strategy relies on 48h skim milk fermentations with the selected strains, in order to i) follow milk acidification, ii) harvest cells to monitor growth and perform PCR-based detection of peptidase genes and iii) purify and extract peptides for quantification, evaluation of their apparent molecular weight distribution and identification of sequences by peptidomics approach.

Multiparametric analysis of fermented cow milk

Correlation between strain growth and acidification capacity in UHT skim milk

Skim milk was inoculated with the selected strains; bacterial growth and pH were measured after 48 hours incubation. Two clearly distinct growth behaviors were observed (Table 1, upper and lower parts): L. helveticus strains showed a much faster growth, reaching an OD600 mean of 3.03 ± 0.322 whereas L. delbrueckii strains only reached a mean of 0.856 ± 0.234.

Acidification capacity of the selected strains was evaluated by pH determination after 48h of culture. The initial pH of NISM (control) was 6.5 and reached 6.45 ± 0.05 after 48h of incubation (Table 1). L. helveticus and L. delbrueckii strains displayed two distinct acidification capabilities. In almost all cases L. helveticus strains led to a pH decrease of minimum 3 pH units, reaching a mean of 3.43 ± 0.08, whereas the pH of milk inoculated with L. delbrueckii decreased by a maximum of 2 pH units, reaching a mean of 4.57 ± 0.07. L. helveticus strains were therefore more efficient than L. delbrueckii strains for both bacterial growth and acidification capacity. As expected, all milk fermentations led to milk protein coagulation because of casein precipitation at pH lower than 4.6 [START_REF] Raak | Enzymatic Cross-Linking of Casein Facilitates Gel Structure Weakening Induced by Overacidification[END_REF].

Increase of peptide concentration in inoculated UHT skim milk

Peptides in fermented milk were extracted by protein precipitation and purified by C18 solid phase extraction. Assessment of peptide concentrations using the FC reagent showed, for all tested strains, an enrichment after 48h fermentation (Table 1). Peptide concentration in NISM was 2.06 g.L -1 ± 0.62 after 48h incubation while inoculated milks showed concentrations higher than 4 g.L -1 in all cases but one. On average, L. helveticus strains generate more peptides than L. delbrueckii strains (mean 5.17 g.L -1 ± 0.92 vs. 4.96 g.L -1 ± 0.74). The two best strains, with a final peptide concentration higher than 6 g.L -1 , are two L. helveticus (25c, 45a).

The size of the peptides generated during fermentation was analyzed by SEC. The chromatogram presented (Figure 1C) depicts the peptide apparent molecular weight distribution. Roughly, the SEC profiles showed two distinct regions: a region of high molecular weight proteins/peptides (HMWPs: fermentation products with an apparent molecular weight higher than 1,700 Da) and a region with low molecular weight peptides (LMWPs: products with an apparent molecular weight lower than 1,700 Da). Profiles were analyzed by integration of the two distinct molecular weight ranges (HMWPs, > 1700 Da and LMWPs , < 1700 Da) and were expressed (Table 1) as % of total area or as % of control (NISM). Proportions for NISM were 85.62% ± 2.72 and 14.39% ± 2.72 of HMWPs and LMWPs respectively. As expected, inoculation and fermentation of milk with lactic acid bacteria lead to an increase of LMWPs, and a decrease of HMWPs compared to NISM (Figure 1C), reflecting the proteolytic activity of the strains against protein and large peptides. L. helveticus and L. delbrueckii strains displayed two distinct peptide production capabilities: the increase of LMWPs proportion compared to NISM was, on average, 1.5 fold for L. delbrueckii and 4.2 fold for L. helveticus (Table 1, last column). The highest activity was observed for strain 25c (% of LMWPs increased 6.7 times compared to NISM).

Strain growth and acidification capacities correlate with peptide production

A first PCA was applied to evaluate correlations between the following parameters: proteolysis diameter on milk agar, cell growth, acidification and peptide content evaluation (FC reagent-based peptide quantification, total area of SEC profiles, %HMWPs, %HMWPs/NISM, %LMWPs and %LMWPs/NISM). Figure 2A shows the obtained variable map. The principal components Dim1 and Dim2 explained 60.82% and 23.56% of the total variance, respectively. A positive correlation was observed between acidification and HMWP proportion (%HMWPs). Moreover strain growth was correlated with a set of variables linked to peptide production (FC reagent-based peptide quantification, %LMWPs and %LMWPs/NISM) and with proteolysis diameter on milk agar plates (Table 1). Acidification and HMWP content were negatively correlated with variables such as growth, LMWP proportion and skim milk agar test results. Regarding the individuals map (Figure 2B), all tested strains were clearly separated from the NISM control. L. delbrueckii strains formed an homogeneous group, clustered together on the negative sides of Dim1 and Dim2 whereas L. helveticus strains were scattered along the positive sides of Dim1 and Dim2.

CEP gene distribution

The distribution of CEP genes was analyzed by PCR in the selected strains and in a control L. helveticus strain (CNRZ32) known to possess all four prtH genes (prtH1 to prtH4). Results are shown in Table 2. The prtB gene was detected in all L. delbrueckii strains but not in L. helveticus strains. The prtH1 gene was detected in 6 L. helveticus strains (16b, 25b, 25c, 49d, 58b and 60b) and, as expected, in CNRZ32 control strain. Interestingly, all L. helveticus strains presenting the prtH1 gene were isolated from mare fermented milk. prtH2 and prtH3 genes were detected in 3 L. helveticus strains (45a, 46b and 47d) and in the CNRZ32 control strain. The prtH4 gene was detected in the CNRZ32 control strain but was not found in any of the Mongolian strains.

Peptide heterogeneity in fermented milk

Peptides were extracted from fermented milk by C18 SPE and analyzed by peptidomics approach. Analysis of the control sample (NISM) allowed defining the base level as 596 peptides belonging to 22 proteins. Nevertheless, after filtering compared to the cow milk protein library (http://www.dallaslab.org/resources -CowMilkProteinDatabase.fasta), 12 cow milk proteins (αS1-casein, β-casein, αS2-casein, GLYCAM1, κ-casein, BTN1A1, βlactoglobulin, FGFBP, lactoperoxidase, osteopontin, perilipin, TRIP6) were finally identified in the NISM control sample. 32% of all peptides were from β-casein (data not shown). In comparison, between 150 and 722 peptides were identified in milks fermented with the selected strains. In these fermented milks, β-casein peptides accounted for 44 to 67% of all sequences identified.

Because of the large proportion of β-casein peptides identified in milk cultures, and the natural high abundance of β-casein in milk, we decided to focus the analysis on this protein.

Figure 3 shows the heat maps obtained for two strains with different CEP profiles, 60b (prtH1) and 47d (prtH2/prtH3). In such type of graphical representation, all peptides identified are repositioned (multiple horizontal grey bars in the middle of Figure 3) on the βcasein backbone, clearly illustrating the β-casein proteolysis patterns. Very few peptides came from the β-casein N-terminal region (signal peptide excluded); the vast majority of them came from a part of β-casein covering amino acids 56 to 224. The repartition of peptides on the βcasein backbone defines a specific proteolysis signature. For the two strains shown as examples (60b, 47d) these signatures were clearly different: the preferential areas of cleavage (dark grey regions on the heat maps) were not situated in the same regions of the protein.

β-casein heat maps are correlated with CEP gene distribution

In order to compare the β-casein proteolysis signatures obtained for each of the 15 strains, a second PCA was performed which included the heat map data. Figure 4A represents the individuals map of this PCA. The principal components Dim1 and Dim2 explained 43.11% and 31.43% of the total variance, respectively. The CEP gene profile of the strains is reported on the figure. Interestingly, the PCA clearly grouped the strains according to their CEP gene profile. As expected, NISM condition (no strain = no prt gene) stands apart from the other points and only depends on Dim2. L. delbrueckii strains, which possess the prtB gene, were grouped on the positive sides of Dim1 and Dim2. L. helveticus strains were scattered along Dim1 and on the negative side of Dim2. Interestingly, strains harboring prtH1 gene were separated from strains harboring prtH2 and prtH3 genes. prtH1 strains showed a higher dispersion, one half on the negative side of Dim1 (strains 16b, 25b and 60b) and the second half on its positive side (strain 25c, 49d and 58b). This dispersion may reflect the presence of different prtH1 gene variants in our isolates, which cannot be distinguished by our PCR approach.

β-casein is composed (before physiological maturation) of 209 amino acids. The variables map of the PCA (Figure 4B) was constructed using 209 variables corresponding to the occurrence frequency of these 209 amino acids. Interestingly, the most implicated variables on Dim1 and Dim2 (defined by cos² > 0.85) were not corresponding to amino acids randomly distributed in the protein. Instead, different specific regions on the β-casein sequence were identified. This observation allowed to distinguish the β-casein areas that were differentially hydrolyzed by the strains. Dim2 corresponded to the Nt region of β-casein (hydrophilic region), whereas Dim1 corresponded to different areas spread along β-casein until the Ct region (hydrophobic parts of the protein) (Figure 4C). In fine, this representation illustrates the protein regions preferentially hydrolyzed during fermentation.

Discussion

Strain identification and screening for proteolytic activity

170 strains were isolated from various Mongolian dairy products, and identified by amplification and 16S rDNA sequencing. Strain identification showed a predominance of 3 Lactobacillus species: L. helveticus, L. delbrueckii and L. kefiranofaciens. This observation is in agreement with previous reports. Several authors have reported that L. helveticus was the dominant LAB in mare fermented milk (airag) [START_REF] Miyamoto | Distinctive proteolytic activity of cell envelope proteinase of Lactobacillus helveticus isolated from airag, a traditional Mongolian fermented mare's milk[END_REF][START_REF] Watanabe | Diversity of lactic acid bacteria and yeasts in Airag and Tarag, traditional fermented milk products of Mongolia[END_REF]. Yu et al. (2011) obtained 668 isolates from different fermented milk products and reported a majority of L. delbrueckii and L. helveticus strains. L. kefiranofaciens was also a species predominantly recovered from Mongolian dairy products [START_REF] Watanabe | Diversity of lactic acid bacteria and yeasts in Airag and Tarag, traditional fermented milk products of Mongolia[END_REF]. However, our findings show that the prevalence of these species depends on the geographical area: L. kefiranofaciens, predominant in dairy samples collected in the Gobi desert, was not found in samples from the Altaï mountains. L. helveticus was more represented in the Altaï samples than in samples collected in the two other regions.

The proteolytic activity of the strains was screened on milk agar. L. helveticus strains showed the highest proteolytic activity of all tested strains. Actually, L. helveticus is considered as the species displaying the highest proteolytic activity among all Lactobacillus [START_REF] Griffiths | Lactobacillus helveticus: The proteolytic system[END_REF]. The other species showing a significant proteolytic activity in our study was L. delbrueckii. Again this is in agreement with previous reports; L. delbrueckii is used as fermentation starter for initiating proteolysis in yogurt production [START_REF] Pailin | Detection of extracellular bound proteinase in EPS-producing lactic acid bacteria cultures on skim milk agar[END_REF]. L. kefiranofaciens strains in our collection only displayed moderate proteolytic activities (data not shown).

Dairy products collected in the Altaï region are therefore characterized by a more proteolytic microbiota (mostly L. helveticus) than products from other areas. These products may therefore present specific properties and potential health benefits (more bioactive peptides).

On the other hand, products collected in the Gobi desert are hosting a large variety of L. kefiranofaciens strains, known to produce kefiran, an exopolysaccharide with beneficial effects on immune system function and cardiovascular health (Maeda, Zhu, Suzuki, Suzuki & Kitamura, 2004). Studies investigating the link between consumption of these products and protection of the local populations against specific pathologies could be undertaken.

Differences in fermentation capacities between L. helveticus and L. delbrueckii strains

The fermentation capacity of L. helveticus and L. delbrueckii strains was investigated by measuring growth, milk acidification and peptide production 48 hours after UHT skim milk inoculation. Overall, L. helveticus strains reached a higher cell density after fermentation than L. delbrueckii strains. As expected, this high bacterial density resulted in a strong acidification of the medium, higher than for samples inoculated with L. delbrueckii, which is in agreement with previous reports [START_REF] Gandhi | Cell growth and proteolytic activity of Lactobacillus acidophilus, Lactobacillus helveticus, Lactobacillus delbrueckii ssp. bulgaricus , and Streptococcus thermophilus in milk as affected by supplementation with peptide fractions[END_REF].

A PCA was used to determine the correlations between all analyzed parameters. PCA are useful to synthesize the information carried by a large number of quantitative variables [START_REF] Lacou | An in silico approach to highlight relationships between a techno-functional property of a dairy matrix and a peptide profile[END_REF]. [START_REF] Bezerra | Proteolysis in goat "coalho" cheese supplemented with probiotic lactic acid bacteria[END_REF], for instance, have used a PCA to study the proteolysis occurring in goat "coalho" cheese. In our study, the variables In fine, the L. helveticus species is characterized by a high intraspecific diversity of CEP profiles. Fermentation of milk with different L. helveticus strains may therefore lead to an important heterogeneity of peptides (Sadat-Mekmene, Genay, Atlan, Lortal, & Gagnaire, 2011a).

Correlation between CEP profiles and β-casein cleavage patterns

β-casein is the main substrate of Lactobacillus proteases [START_REF] Jensen | Variation in aminopeptidase and aminotransferase activities of six cheese related Lactobacillus helveticus strains[END_REF].

It is thought that different CEPs can have very different cleavage specificities [START_REF] Oberg | Diversity in specificity of the extracellular proteinases in Lactobacillus helveticus and Lactobacillus delbrueckii subsp. bulgaricus[END_REF]. However it is somewhat difficult to compare the results of separate publications: depending on the study, caseins are indeed used either as fragments, as purified forms, or as micelles in milk, and this of course may influence the ability of proteases to access specific cleavage sites in the protein (Sadat-Mekmene et al., 2011a). [START_REF] Oberg | Diversity in specificity of the extracellular proteinases in Lactobacillus helveticus and Lactobacillus delbrueckii subsp. bulgaricus[END_REF], who assessed the proteolytic activity of 14 L. delbrueckii strains based on the cleavage of an αs1-casein fragment, reported several distinct digestion patterns.

In contrast, in our study, the β-casein cleavage patterns of all L. delbrueckii strains revealed some common cleavage sites. Our strains may be genetically closer to each other; but our observation may also result from the fact that we worked with native milk β-casein (micelle configuration), and not with a fragment of this protein.

The PCA separated L. delbrueckii strains from L. helveticus strains, indicating dissimilarities on β-casein hydrolysis patterns. Sadat-Mekmene et al. (2011b) obtained similar results when incubating purified β-casein with either L. helveticus CP790 or L. delbrueckii CRL581 strains (Sadat-Mekmene et al., 2011b). It is notable however that three of our L. helveticus strains (25c, 49d, 58b, all prtH1) were located very close to the prtB strains (L. delbrueckii). The three other prtH1 strains (16b, 25b, 60b) were located close to the prtH2/prtH3 strains (Figure 4A). Therefore, it appears that prtB and prtH2/prtH3 genotypes are associated with clearly distinct proteolytic activities; the situation for prtH1 is however more complex, which may be related to the existence of different alleles for this gene, not differentiated by our PCR [START_REF] Genay | PrtH2, Not prtH, Is the ubiquitous cell wall proteinase gene in Lactobacillus helveticus[END_REF]. Interestingly, while PrtH2 and PrtH3 share 58% of amino acid strict identity, PrtH1 only shares 32 and 34% of amino acid strict identities with PrtH2 and PrtH3, respectively [START_REF] Broadbent | Genetic diversity in proteolytic enzymes and amino acid metabolism among Lactobacillus helveticus strains[END_REF]. It is remarkable that despite this high genetic variation, some prtH1 strains still display a proteolytic activity which appears very similar to the activity of prtH2/prtH3 strains.

As reported in section 3.6, the Dim2 of PCA depicted the hydrophilic region of β-casein while Dim1 depicted the hydrophobic region of β-casein (Figure 4). The selection of the most implicated variables from PCA allowed to highlight the main important differences between conditions. The NISM condition was used to define the basal proteolytic pattern due to endogenous milk peptidases such as plasmin and cathepsin D [START_REF] Dallas | Peptidomic analysis reveals proteolytic activity of kefir microorganisms on bovine milk proteins[END_REF]Sadat-Mekmene et al., 2011b). As previously reported for human milk (Khaldi et al., 2014), our results indicate that these endogenous milk peptidases hydrolyzed the hydrophilic part (Nt region) of β-casein during the our experimental conditions (37°C and 48h): the point corresponding to NISM was indeed placed by the PCA right and high on the Dim2 axis (Figure 4A).

All the L. delbrueckii strains were closely grouped together, reflecting a high homogeneity in terms of proteolytic capacities; they were placed high on Dim1 and close to 0 on the Dim2 axis, suggesting that hydrophobic areas along β-casein were preferentially hydrolyzed by these strains. In contrast, strains of L. helveticus were scattered along the Dim1 axis, reflecting a much higher heterogeneity; strains harboring prtH2 and prtH3 were located far in the negative side of Dim1, showing a preferential hydrolysis of hydrophobic β-casein regions.

Our PCA therefore revealed the existence of different proteolytic activities, resulting in different cleavage patterns of β-casein. The capacity to cleave β-casein in its more hydrophobic regions, which are more difficult to access when the protein is under micellar form [START_REF] Dumas | Structure primaire de la caseine bovine: Sequence complete[END_REF], is a significant species advantage.

Conclusion

In this study, we used an original approach to screen a collection of 170 Lactobacillus strains for their ability to produce peptides derived from milk proteins. All strains had been isolated from dairy products collected in three different regions of Mongolia. Fifteen strains were first selected by milk agar-well diffusion test; bacterial growth, milk acidification and production of peptides were then evaluated after 48h fermentation of skim milk. The CEP gene profile of each strain was also determined. A multiparametric analysis demonstrated that growth and milk acidification capacity correlated closely with peptide production. The CEP gene distribution was different in Mongolian L. helveticus strains than in European and North American strains, and PCA revealed that presence of specific CEPs was linked to specific βcasein hydrolysis patterns. The tested strains display, at various levels, the capacity to cleave β-casein on specific sites buried inside the micellar form of caseins (hydrophobic regions). This multiparametric analysis allowed to select strains harboring i) distinct proteolytic patterns, ii) different CEP gene profiles, iii) high growth rate in milk and iv) high acidification capacity. The selected strains can be used as new fermentation starters, able to generate milk peptides with potential bioactivities. The analysis of such bioactivities, using in silico approach by comparison with milk bioactive peptide database [START_REF] Nielsen | Milk bioactive peptide database: A comprehensive database of milk protein-derived bioactive peptides and novel visualization[END_REF], is ongoing. Diameters obtained on skim milk agar during the first phase of screening for the 15 selected strains.

Results from the second phase of screening on the selected strains obtained after UHT skim milk fermentation at 37°C for 48h: growth, acidification, peptide quantification with FC reagent and integration results from chromatogram obtained by SEC ((high molecular weight proteins/peptides (HMWPs); low molecular weight peptides (LMWPs)). NISM: Non-Inoculated UHT Skim Milk. *Statistical analysis using one way ANOVA and Tukey's post-hoc test. Mean values in the same column with different uppercase subscripts differ significantly (p < 0.05). Distribution of CEP genes detected by PCR (prtH1, prtH2, prtH3, prtH4 and prtB) in the 15 selected strains and in the reference strain CNRZ32. Five primer set were used for the PCR.

(

  2011). Diversity of lactic acid bacteria associated with traditional fermented dairy products in Mongolia. Journal of Dairy Science, 94(7), 3229-3241. https://doi.org/10.3168/jds.2010-3727 Highlights • Mongolian strains show specific skim milk fermentation capacities • Cell envelope associated peptidase genes distribution is different in Mongolian Lactobacillus helveticus strains than in Western strains • Peptide production is closely correlated with growth and acidification patterns • Patterns of β-casein hydrolysis are correlated with cell envelope associated peptidase genes profiles

Figure 1 :Fig. 1 .

 11 Figure 1 :

Figure 2 :Fig. 2 .

 22 Figure 2 :

Figure 3 :

 3 Figure 3:

Figure 4 :Fig. 4 .

 44 Figure 4:

Table 1 .

 1 Results of the strain screenings .46 ad 225.08 ± 118.3 a 58.56 ± 14.51 ef 42.19 ± 14.63 b 43.87 ± 29.24 bc 263.68 ± 77.95 ab 58b L. helveticus 37 ± 1.7 ab 3.403 ± 0.394 ef 3.44 ± 0.01 ab 5.95 ± 1.36 cd 263.89 ± 21.7 a 59.24 ± 8.93 df 40.79 ± 8.95 bc 49.48 ± 4.74 c 322.98 ± 93.51 ab 60b L. helveticus 35 ± 1 ac 2.81 ± 0.257 de 3.49 ± 0.11 ab 5.09 ± 0.91 ad 303.15 ± 165.9 a 46.4 ± 6.48 bd 59.41 ± 8.02 cd 29.10 ± 7.85 ac 526.87 ± 274.97 ab

	Strain	Species	Ø milk agar (mm)	Growth (OD)	Acidification (pH)	FC reagent (g.L -1 )	Total area (mAU)	%HMWPs	%LMWPs %HMWPs/NISM %LMWPs/NISM
	NISM	×	×	×	6.45 ± 0.05 e	2.06 ± 0.62 a	347.51 ± 42.8 a 85.62 ± 2.72 g	14.39 ± 2.72 a	100 d	100 a
	23b	L. delbrueckii 33 ± 3 ac	0.739 ± 0.01 ac	4.52 ± 0.02 c	4.86 ± 0.21 ad	101.48 ± 31.8 a 61.37 ± 3.44 df 40.06 ± 5.13 bc	19.74 ± 5.66 ab	123.32 ± 52.81 a
	33c	L. delbrueckii 32 ± 2 ac 0.825 ± 0.126 ac	4.64 ± 0.02 c	5.54 ± 0.61 bcd 141.65 ± 27.2 a 54.83 ± 1.28 cde 45.98 ± 2.43 bc	24.69 ± 4.3 ab	194.64 ± 47.3 a
	34b	L. delbrueckii 39 ± 3.5 c 0.987 ± 0.03 bc	4.52 ± 0.02 c	5.32 ± 0.21 ad	115.96 ± 4.35 a 59.07 ± 3.17 df 41.08 ± 3.35 bc	21.79 ± 0.37 ab	141.83 ± 16.75 a
	48b	L. delbrueckii 30 ± 0 ac 0.473 ± 0.118 ab	4.68 ± 0.01 c	3.59 ± 0.19 ad	123.56 ± 12.3 a 57.26 ± 1.76 de 44.14 ± 2.52 bc	22.50 ± 1.65 ab	162.65 ± 24.69 a
	50b	L. delbrueckii 30 ± 1.5 ab 0.973 ± 0.199 bc	4.51 ± 0 c	4.89 ± 0.3 bcd	114.70 ± 12.5 a 58.12 ± 6.22 de 40.90 ± 4.26 bc	21.12 ± 1.81 ab	140.37 ± 28.22 a
	56a	L. delbrueckii 28 ± 1.1 a	1.14 ± 0.124 c	4.52 ± 0.01 c	5.58 ± 0.29 bcd 131.81 ± 35.8 a 58.35 ± 3.93 df 42.32 ± 5.03 bc	24.21 ± 4.98 ab	169.25 ± 65.55 a
	Average L. delbrueckii 32 ± 3.8	0.856 ± 0.234	4.57 ± 0.07	4.96 ± 0.74	121.53 ± 14.1	58.17 ± 2.1	42.4 ± 2.3	22.34 ± 1.9	155.34 ± 25.39
	16b	L. helveticus 31 ± 1.5 ab	2.79 ± 0.03 de	3.56 ± 0.05 b	4.49 ± 0.02 ac	184.11 ± 46.7 a 34.14 ± 1.98 ab 73.24 ± 2.06 de	19.89 ± 4.18 ab	400.81 ± 103.87 ab
	25b	L. helveticus 32 ± 1 ac	3.18 ± 0.194 df	3.48 ± 0.1 ab	4.17 ± 0.68 ad	155.52 ± 4.92 a	38.83 ± 2.2 ab	64.62 ± 1.89 de	15.39 ± 0.6 a	503.42 ± 297.66 ab
	25c	L. helveticus 32 ± 2 ac	2.62 ± 0.138 d	3.45 ± 0.17 ab	6.78 ± 1.63 d	268.7 ± 121.8 a	20.13 ± 3.6 a	85.7 ± 7.47 e	16.91 ± 6.77 a	678.54 ± 309.89 b
	45a	L. helveticus 32 ± 4 ac	3.57 ± 0.12 f	3.32 ± 0.04 a	6.16 ± 1.32 bcd	289.95 ± 52 a	40 ± 3.74 ad	60.1 ± 3.77 cd	37.35 ± 11.95 ac	514.7 ± 86.67 ab
	46b	L. helveticus 35 ± 2.6 ac	3.03 ± 0.02 df	3.39 ± 0.1 ab	4.44 ± 0.17 ad	139.24 ± 6.51 a	37 ± 3.1 abc	74.2 ± 3.67 de	16.37 ± 0.74 a	307.41 ± 27.56 ab
	47d	L. helveticus 32 ± 1 ac	3.12 ± 0.326 df	3.34 ± 0.04 ab	5.01 ± 0.8 ad	170.7 ± 52.9 a 44.73 ± 12.24 bd 58.91 ± 16.94 cd	23.93 ± 7.31 ab	297.41 ± 114.8 ab
	49d 4.46 ± 0Average L. helveticus 33.4 ± 2.1 L. helveticus 35 ± 4.3 ac 2.72 ± 0.05 d 3.43 ± 0.01 ab 3.03 ± 0.322 3.43 ± 0.08 5.17 ± 0.92	222.26 ± 61.7	42.11 ± 12.1	62.13 ± 14.6	28.03 ± 12.8	423.98 ± 139.63

Table 2 .

 2 CEP gene distribution

	Strain	Specie prtH1 prtH2 prtH3 prtH4 prtB Isolated from
	23b L. delbrueckii			+	Cow milk
	33c L. delbrueckii			+	Cow milk
	34b L. delbrueckii			+	Mare milk
	48b L. delbrueckii			+	Cow milk
	50b L. delbrueckii			+	Yak milk
	56a L. delbrueckii			+	Yak milk
	CNRZ32 L. helveticus +	+	+	+	CIRM (INRA)
	16b	L. helveticus +				Mare milk
	25b	L. helveticus +				Mare milk
	25c	L. helveticus +				Mare milk
	49d	L. helveticus +				Mare milk
	58b	L. helveticus +				Mare milk
	60b	L. helveticus +				Mare milk
	45a	L. helveticus	+	+		Yak milk
	46b	L. helveticus	+	+		Yak+Cow milk
	47d	L. helveticus	+	+		Camel milk

Acknowledgements

This research was supported by a grant from BPI France/Métropole Européenne de Lille/Région Hauts de France (DOS0020745) and by a grant from Association Nationale Recherche et Technologie (ANRT). This work was supported by the Hauts-de-France region funding through the ALIBIOTECH research program. The HPLC-MS experiments were performed on the REALCAT platform funded by a French governmental subsidy managed by the French National Research Agency (ANR) within the frame of the "Future Investments'program (ANR-11-EQPX-0037). The Hauts-de-France region and the FEDER, the Ecole Centrale de Lille and the Centrale Initiatives Foundation are also warmly acknowledged for their financial contributions to the acquisition of REALCAT platform equipment. We are grateful to Prof. Anil Kulkarni (University of Texas Medical School, Houston, USA), Dr Anarkhuu Bold-Erdene and Dr Amarbat Baatarnum (National Cancer

Center, Ulaanbaatar, Mongolia) for their support during the initiation of the project. The inhouse tool used for heat maps and peptide mappings was designed by Yohan Guegaden. We thank Aurélie Mateos for helpful discussions.

map showed that high cellular growth and high acidification capacities were linked to a more abundant peptide production. Such correlation makes sense, especially in the case of L. helveticus species, where the efficient proteolytic system fulfill their nitrogen requirement leading to their high growth in pure milk. LAB is a heterogeneous microorganism group with multiple auxotrophies for amino acids and the capacity to rapidly degrade milk proteins, and thereby release these amino acids, is a major asset [START_REF] Courtin | Cell-wall proteinases PrtS and PrtB have a different role in Streptococcus thermophilus/Lactobacillus bulgaricus mixed cultures in milk[END_REF][START_REF] Griffiths | Lactobacillus helveticus: The proteolytic system[END_REF][START_REF] Savijoki | Proteolytic systems of lactic acid bacteria[END_REF].

CEP distribution in the fifteen selected strains

All L. delbrueckii strains were found to possess the prtB gene, which is so far the only CEP described for this species [START_REF] Hou | Effect of culturing conditions on the expression of key enzymes in the proteolytic system of Lactobacillus bulgaricus[END_REF]. L. helveticus strains showed more heterogeneity, being divided into two groups: 6 strains with prtH1 gene, and 3 strains with the combination prtH2/prtH3. The CEP profile does not seem to be related to the geographical origin of the strains: among the six prtH1 strains, one was from Gobi, two from the Hangai region, and three from Altaï. Rather, the CEP profiles of the strains seem to be related to the type of sample from which they were isolated: all prtH1 strains were isolated from airag samples (horse milk) whereas the three strains harboring prtH2/prtH3 were isolated from other sources (yak and camel milks). The presence of prtH1 in L. helveticus strains isolated from airag had already been reported in another study [START_REF] Miyamoto | Distinctive proteolytic activity of cell envelope proteinase of Lactobacillus helveticus isolated from airag, a traditional Mongolian fermented mare's milk[END_REF]. However, the combination of the two CEP genes prtH2 and prtH3 is unusual. Investigating CEP genes distribution in 50 L. helveticus strains of European and North American origin, [START_REF] Broadbent | Genetic diversity in proteolytic enzymes and amino acid metabolism among Lactobacillus helveticus strains[END_REF] reported only three strains with the same profile. This seems to confirm previous observations that Mongolian L. helveticus strains differ from Western strains in terms of CEP gene profiles [START_REF] Miyamoto | Distinctive proteolytic activity of cell envelope proteinase of Lactobacillus helveticus isolated from airag, a traditional Mongolian fermented mare's milk[END_REF].

Conflict of interest

CR, MF and MV are affiliated with VF Bioscience SAS, France, a company commercializing lactic acid bacteria strains for use as probiotics or dairy starters.


Rabie Aloui
email: rabie.aloui@enim.rnu.tn

Walid Larbi
email: walid.larbi@cnam.fr

Mnaouar Chouchane
email: mnaouar.chouchane@enim.rnu.tn

Global Sensitivity Analysis of Piezoelectric Energy Harvesters

Keywords: Global Sensitivity Analysis, Elementary Effects, Piezoelectric Energy Harvesting, Finite Element Method

Vibration energy harvesting using the direct effect of piezoelectricity has attracted increasing attention during the last two decades. Different modeling techniques have been applied to describe the electromechanical coupling effect of a piezoelectric harvester and to predict its electrical output.

This study aims to identify the most important properties of both harvester substrate material and piezoelectric material that cause uncertainty in the predicted performances of the harvester. Global sensitivity analysis, applied in this paper, is a promising method used to identify systems parameters which have significant impact on the system output. In this paper, the Elementary Effects method (EEs), a particular implementation of the global sensitivity method, is used to identify the impact of substrate and piezoelectric material properties on the voltage frequency response function of a typical bimorph piezoelectric energy harvester with fixed geometry. With a small number of model evaluations at selected ranges of material properties,

Introduction

Energy harvesting using piezoelectric material has attracted considerable interest in the last decade. Piezoelectric material displays converting capabilities due to its direct and indirect electromechanical effect. Several approaches have been developed to model Piezoelectric Energy Harvesters (PEHs). In particular, two methods have been accepted and are the most widely used in the scientific community: (i) the analytical distributed parameter model introduced by Erturk and Inman [START_REF] Erturk | Piezoelectric energy harvesting[END_REF][START_REF] Erturk | A Distributed Parameter Electromechanical Model for Cantilevered Piezoelectric Energy Harvesters[END_REF] is applied for unimorph and bimorph harvesters to obtain a modal solution of second-order ordinary differential equations, (ii) the finite element method, adapted by De Marqui Junior et al. [START_REF] Marqui Junior | An electromechanical finite element model for piezoelectric energy harvester plates[END_REF][START_REF] Rosa | Modeling and Analysis of a Piezoelectric Energy Harvester with Varying Cross-Sectional Area[END_REF], uses standard discretization to provide discrete models with less restrictive assumptions on the global electrical variables. The two modeling approaches are particularly useful as they have been amply tested showing close correlation with experimental data.

In the frame of numerical approaches, sensitivity analysis is a very useful tool for modeling and design analysis of mechatronic systems. It is applied in various engineering problems such as structural analysis, model updating, design optimization of structures, system control and uncertainty propagation [START_REF] Lasecka Plura | Design sensitivity analysis of frequency response functions and steady state response for structures with viscoelastic dampers[END_REF]. Sensitivity analysis enables the evaluation of the degree of influence of input parameters on the output responses of the modeled system. Furthermore, this method process is widely considered as the most important step in design optimization since significant number of optimization algorithms require the prediction of the evolution of the output responses of the system for different values of design variables [START_REF] Lima | Sensitivity analysis of frequency response functions of composite sandwich plates containing viscoelastic layers[END_REF].

Sensitivity analysis techniques can be classified into local and global methods. Local sensitivity analysis, often referred to as a one factor at time analysis, is based on the approximation of partial derivatives to assess how uncertainty in one factor affects the model response keeping other factors fixed at their nominal values. Global Sensitivity Analysis (GSA), on the other hand, offers a comprehensive approach to the studied model, since it evaluates the effect of factors that are varying within the considered multidimensional input space [START_REF] Campolongo | From screening to quantitative sensitivity analysis. A unified approach[END_REF].

In recent published research, Sharp and Brooks [START_REF] Sharp | Sensitivities of frequency response functions of linear dynamic systems to variations in design parameter values[END_REF] introduced expressions for the sensitivities of the Frequency Response Functions (FRFs) of linear constant coefficient systems with distinct eigenvalues to perturbations of the design parameters suggesting several engineering applications of the developed theory. Lasecka-Plura et al [START_REF] Lasecka Plura | Design sensitivity analysis of frequency response functions and steady state response for structures with viscoelastic dampers[END_REF] studied the sensitivity of the FRF and steady-state vibration response of planar frames with viscoelastic dampers. A comparison was made between the first-order, the second-order and the exact solution sensitivities. Analysis show that the second order sensitivity gives results which are very close to the exact solution, when the perturbation of the parameter is smaller than 20%. Lima et al [START_REF] Lima | Sensitivity analysis of frequency response functions of composite sandwich plates containing viscoelastic layers[END_REF][START_REF] De Lima | Sensitivity analysis of viscoelastic structures[END_REF] investigated the first order sensitivity of complex frequency response functions for composite sandwich plates composed of a combination of fiber-reinforced and elastomeric viscoelastic layers, frequently used for noise and vibration attenuation. The study aims to find the sensitivity of viscoelastic struc-tures to variation of temperature and structural parameters. Ng and Liao [START_REF] Ng | Sensitivity Analysis and Energy Harvesting for a Self-Powered Piezoelectric Sensor[END_REF] used the first order sensitivity to evaluate the performances of sensing and power generation of the piezoelectric energy harvesters to variation of elastic substrate properties for three type of sensors. The reported research results show that for the series bimorph piezoelectric sensor, the voltage sensitivity is the highest and the charge sensitivity is the lowest. For the parallel bimorph piezoelectric sensor, the voltage sensitivity is the lowest and charge sensitivity is highest. A unimorph sensor, however, has relatively medium sensitivities of charge and voltage responses. Ruiz et al [START_REF] Ruiz | Uncertainties propagation and global sensitivity analysis of the frequency response function of piezoelectric energy harvesters[END_REF] analyzed the uncertainty propagation and global sensitivity of the FRFs of the output voltage of piezoelectric energy harvesters. Their approach is based on the identification of Sobol' indices to assess the robustness of the stochastic prediction. Peralta et al [START_REF] Peralta | Bayesian Framework to Quantify Uncertainties in Piezoelectric Energy Harvesters[END_REF] introduced a procedure to update the electromechanical properties of the PEHs based on Bayesian updating techniques and global sensitivity analysis. The introduced approach constitutes a powerful method for the robust design and prediction of the performances of piezoelectric energy harvesters.

In the literature, several GSA methods are applied [START_REF] Pianosi | A Matlab toolbox for Global Sensitivity Analysis[END_REF][START_REF] Chen | Analytical Variance-Based Global Sensitivity Analysis in Simulation-Based Design Under Uncertainty[END_REF][START_REF] Kucherenko | Monte Carlo evaluation of derivative-based global sensitivity measures[END_REF] ranging from qualitative screening techniques [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF][START_REF] Sanchez | Application of sensitivity analysis in building energy simulations: Combining firstand second-order elementary effects methods[END_REF] to quantitative methods based on variance decomposition [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF][START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF]. The Morris method [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF] is a global sensitivity screening technique used to derive sensitivity measures from a set of local derivatives sampled on a grid throughout the parameter space [START_REF] Herman | Technical Note: Method of Morris effectively reduces the computational demands of global sensitivity analysis for distributed watershed models[END_REF].

This technique investigates the model response to a change in the inputs by varying One-At-a-Time (OAT), while keeping all other variables fixed. The local sensitivity measure associated to each factor is called elementary effect and is defined as the ratio between the variation in the model response output and the variation in the input factor. In order to obtain a global sensitivity measure, different elementary effects (local derivatives taken at points sampled throughout the parameter space) for each parameter are estimated and averaged. Two sensitivity measures are carried out, the first one is the mean of the EEs which is refined by Campolongo et al [START_REF] Campolongo | From screening to quantitative sensitivity analysis. A unified approach[END_REF][START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF] consisting in averaging the absolute values of the elementary effects. The second measure is the standard deviation of the EEs. By plotting both sensitivity measures, the Morris method identifies the hierarchical order of influence of the input parameters on the model output.

Most of the previous researches were limited to local sensitivities for variation of system parameters in the vicinity of the considered nominal configuration. Generally, an interval variation of system parameter is more acceptable than a nominal parameter, hence the importance of applying global sensitivity. Morris method computes global sensitivity in an efficient way. Previously published research work does not seem to have applied this method in sensitivity analysis of piezoelectric harvesters. In this paper, we propose combining the implementation of the PEH model with Morris elementary effect estimation method to determine the effect of uncertain electromechanical parameters of a bimorph piezoelectric energy harvester on the output voltage for short-circuit and open-circuit conditions. Furthermore, the effective electromechanical coupling factor is monitored for the model evaluations used to compute elementary effects. In order to predict the output voltage of the studied harvester, a finite element formulation of a laminated piezoelectric harvester is applied. The mathematical development of the model is given in Section 2. This model is used in Section 3 to derive the complex FRFs of a bimorph piezoelectric energy harvester considered as a case study in this paper. In Section 4, the Morris method is applied using the modulus of the voltage FRF of the harvester as an output and a set of uncertain physical parameters. The considered numerical example presented in Section 5 shows the robustness of the EEs method compared to others GSA methods. Furthemore, the global sensitivities of the voltage FRF of a bimorph piezoelectric harvester based on elementary effects have been compared with those computed using the variance decomposition based sensitivity analysis method described in reference [START_REF] Ruiz | Uncertainties propagation and global sensitivity analysis of the frequency response function of piezoelectric energy harvesters[END_REF] for different configurations of the PEHs.

Finite Element Modeling of a Piezoelectric Composite Beam

Considering a laminated piezoelastic beam composed of K layers including P piezoelectric layers excited under base motion at the clamped end as shown in Figure 1. In the absence of mechanical dissipative effects, the extended Hamilton's principle for an electromechanical body, between the two known times t 1 and t 2 , is [START_REF] Amini | Finite element modeling of functionally graded piezoelectric harvesters[END_REF]:

t 2 t 1 (δT -δU + δW nc )dt = 0 ( 1
)
where δT is the first variation of the kinetic energy, δU is the first variation of the internal energy including the piezoelectric effect contribution and δW nc is the virtual work done by the non-conservative electric charge components.

The thin piezoelectric layers of the cantilever beam are poled in the thickness direction with an electrical field applied parallel to this polarization direction. Such a configuration is characterized in particular by the electromechanical coupling between the axial strain and the transverse electrical field [START_REF] Thomas | Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients[END_REF]. The reduced constitutive equations corresponding to the p th piezoelectric layer is written as follows, p ∈ {1, • • • , P }: 3 is normal to the electrodes of the piezoelectric layer and its intensity is assumed to be uniform in the piezoelectric thickness direction. The superscript (p) refers to the piezoelectric material. Moreover, the thickness of the electrodes is much smaller than the thickness of piezoelectric layer so that the equipotentiality assumption in the piezoelectric electrodes is well satisfied. Thus the electric field can be expressed as:

σ (p) 1 = c (p) 11 ε (p) 1 -e (p) 31 E (p) 3 (2)
E (p) 3 = - V (p) h (p) (4)
where V (p) is the potential difference between the upper and the lower electrode surfaces and h (p) = z p -z p-1 is the thickness of the p th piezoelectric layer.

The equations of the coupled electromechanical variational formulation can be rewritten in the following form [START_REF] Thomas | Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients[END_REF][START_REF] Larbi | Finite element formulation of smart piezoelectric composite plates coupled with acoustic fluid[END_REF]:

K k=1 Ω k ρ k .. u x δu x + .. u z δu z dΩ + K k=1 Ω k c k 11 ε 1 δε 1 dΩ + P p=1 V (p) h (p) Ω (p) e (p)
31 δε 1 dΩ = 0

(5)

- P p=1 δV (p) h (p) Ω (p) e (p) 31 ε 1 dΩ + P p=1 δV (p) h (p) Ω (p) (p) 33
V (p) h (p) dΩ + P p=1 δV (p) Q (p) = 0 (6)
where ρ k and Ω k are respectively the mass density and the domain occupied by the k th layer, and Q (p) is the electrical charge of the p th piezoelectric layer. Furthermore, the mechanical displacements fields u x and u z of the laminated harvester under base excitation, based on Euler-Bernoulli beam assumptions, are defined as follow:

u x (x, z, t) = u(x, t) -zθ(x, t) (7)
u z (x, z, t) = w(x, t) = w b (t) + w rel (x, t) ( 8
)
where u is the axial displacement of the center line of the beam, θ = ∂w rel /∂x is the beam section rotation, w is the transverse displacement, w b (t) = w b e jωt is the base displacement and w rel is the relative transverse displacement of clamped-free beam.

Applying the finite element method, the previous formulations uses a standard discretization of N mechanical Degrees Of Freedom (DOFs) and provides less restrictive assumptions on the electrical DOFs: only one electrical DOF for each piezoelectric layer is needed due to the equipotentiality in the electrodes. In this paper, three mechanical degrees of freedom per node are assumed (u, w, θ). Details of the stratified beam finite element discretization are provided in [START_REF] Thomas | Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients[END_REF] for the master element shape functions.

The various terms appearing in the variational formulation of Equations ( 5) and ( 6) in terms of Hamilton's principle are now discussed. The discrete form for each term is provided [START_REF] Amini | Finite element modeling of functionally graded piezoelectric harvesters[END_REF]:

(i)
The kinetic energy variation is:

δT = - K k=1 Ω k ρ k .. u x δu x + .. u z δu z dΩ ⇒ -δU T M m .. U (t) + δU T F(t) (9)
where M m is the global (N × N ) mass matrix, U(t) = Ue jωt is the global (N × 1) vector of mechanical coordinates displacements, F(t) = Fe jωt is the global (N × 1) forcing vector due to base excitation, which can be expressed as an (N × 1) effective mass vector m * multiplied by the base acceleration [START_REF] Marqui Junior | An electromechanical finite element model for piezoelectric energy harvester plates[END_REF][START_REF] Amini | Finite element modeling of functionally graded piezoelectric harvesters[END_REF] as follows:

-

K k=1 Ω k ρ k δw rel ẅb dΩ ⇒ δU T m * ẅb (t) = δU T F(t) (10)
(ii) The internal energy variation is:

δU = K k=1 Ω k c k 11 ε 1 δε 1 dΩ + P p=1 V (p) h (p) Ω (p) e (p)
31 δε 1 dΩ

+ P k=1 δV (p) h (p) Ω (p) e (p) 31 ε 1 dΩ - K k=1 δV (p) h (p) Ω (p) (p) 33
V (p) h (p) dΩ (11)
• The first term of Equation ( 11) is the mechanical contribution to the internal energy variation:

K k=1 Ω k c k 11 ε 1 δε 1 dΩ ⇒ δU T K m U(t) (12)
where K m is the global (N × N ) stiffness matrix.

• The second and third terms in Equation ( 11) are the piezoelectric contributions to the internal energy variation related to the direct and inverse effects:

P p=1 V (p) h (p) Ω (p) e (p) 31 δε 1 dΩ ⇒ δU T K c V(t) (13)
P k=1 δV (p) h (p) Ω (p) e (p) 31 ε 1 dΩ ⇒ δV T K T c U(t) (14)
where

K c = (K (1) c • • • K (P )
c ) is the global (N × P ) electromechanical coupling matrix, V(t) = Ve jωt is the global (P × 1) vector of voltage outputs.

• The last term in Equation ( 11) is the electrostatic contribution to the internal energy variation is:

K k=1 δV (p) h (p) Ω (p) (p) 33
V (p) h (p) dΩ ⇒ δV T K e V(t) (15)
where

K e = diag(C (1) • • • C (P )
) is the diagonal global (P × P ) capacitance matrix, and

C (p) = b × L h (p) (p)
33 is the capacitance of the p th piezoelectric layer.

(iii) The variation of work due to the charges is:

δW nc = P p=1 δV (p) Q (p) ⇒ δV T Q(t) (16)
where Q(t) = Qe jωt is the global (P × 1) vector of electric charge outputs.

The harvested energy is dissipated through a resistive electrical load R.

Using Ohm's law, the following additional equation relates the voltage vector V to the time derivative of the charge vector Q:

V(t) = R . Q (t) (17)
To account for system damping, a proportional viscous damping is assumed and included in the model. The global (N × N ) damping matrix C m is assumed to be a linear combination of the mass and stiffness matrices (the Rayleigh damping), so that:

C m = αM m + βK m ( 18
)
where α and β are real proportional coefficients which are typically determined from the measurement of modal damping.

The general damped electromechanical equations obtained from the previous finite element formulation can be expressed as [START_REF] Shu | An improved analysis of the SSHI interface in piezoelectric energy harvesting[END_REF][START_REF] Larbi | Coupled FEM/BEM for control of noise radiation and sound transmission using piezoelectric shunt damping[END_REF][START_REF] Moita | Active-passive damping in functionally graded sandwich plate/shell structures[END_REF]:

M m .. U (t) + C m . U (t) + K m U(t) + K c V(t) = F(t) (19
)

K e V(t) -K T c U(t) + Q(t) = 0 ( 20
)
where Equation ( 19) corresponding to the mechanical equation of motion includes the inverse effect electromechanical coupling term, whereas Equation [START_REF] Herman | Technical Note: Method of Morris effectively reduces the computational demands of global sensitivity analysis for distributed watershed models[END_REF] corresponding to the electrical equation includes the direct effect electromechanical coupling term introduced by the piezoelectric behavior.

Natural frequencies of the harvester for the short-circuit and the opencircuit conditions

The natural frequencies of the harvester are calculated by solving the free undamped system corresponding to Equation [START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF]. For piezoelectric energy harvesters, two natural frequencies are defined as a characteristic limit bounds for each vibration mode: (i) the short-circuit natural frequencies are calculated by assuming that no potential difference exists across the piezoelectric layers in free vibration. As a result, the electromechanical coupling in Equation ( 19) is omitted, (ii) the open circuit natural frequencies, however, are calculated by assuming that no charge flows in the electrical circuit [START_REF] Larbi | Reduced order finite element formulations for vibration reduction using piezoelectric shunt damping[END_REF]. In practice, these natural frequencies are determined in terms of the load resistance value. The short-circuit condition corresponds to a low resistive load connected to the electrodes of the piezoelectric layers (R → 0), whereas the open-circuit condition is determined for a high load resistance (R → ∞). The main motivation of computing these particular limits is to determine the electromechanical coupling factor, using a classical elastic mechanical formulation, which corresponds to the electromechanical energy ratio.

The normal modes of the short-circuit system are obtained assuming harmonic solutions of the undamped equations of motion [START_REF] Saltelli | Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF] for V = 0

and

F = 0. K m -ω 2 i M m Φ i = 0 (21)
where ω i is the natural frequency and Φ i is eigenvector for the i th mode for short-circuit condition. These modes verify the following orthogonality properties:

Φ T i M m Φ j = δ ij and Φ T i K m Φ j = ω 2 i δ ij ( 22
)
where δ ij is the Kronecker symbol and the mode Φ i is normalized with respect to the structural mass matrix.

The eigenvalue problem for the open-circuit condition is derived by coupling Equations ( 19) and ( 20) substituting Q = 0 and F = 0:

K m -ω 2 i M m Φ i = 0 (23)
where

K m = K m + K c K -1 e K T c
is the stiffness matrix for the open-circuit condition which depends on the piezoelectric system properties, ω i is the natural frequency and Φ i is eigenvector for the i th vibration mode. These modes verify the following orthogonality properties:

Φ T i M m Φ j = δ ij and Φ T i K m Φ j = ω 2 i δ ij (24)

Effective Electromechanical Coupling Factor

The dimensionless Effective Electromechanical Coupling Factor (EEMCF) of the piezoelectric laminated beam for the i th mode can be expressed in terms of the natural frequencies of the short and open-circuit conditions.

This factor defines the exchange of energies of the vibrating structure at the i th mode. It is usually defined, as follows [START_REF] Trindade | Effective Electromechanical Coupling Coefficients of Piezoelectric Adaptive Structures: Critical Evaluation and Optimization[END_REF][START_REF] Larbi | Structural-Acoustic Vibration Reduction Using Switched Shunt Piezoelectric Patches: A Finite Element Analysis[END_REF][START_REF] Pereira Da Silva | Topology optimization of shunted piezoelectric elements for structural vibration reduction[END_REF]:

k 2 ef f,i = ω 2 i -ω i 2 ω i 2 (25)

The Bimorph Piezoelectric Energy Harvester

The finite element formulation, presented in Section 2, can be easily adapted to the case of a bimorph piezoelectric energy harvester considered as a case study in this paper. The harvester is composed of two identical piezoelectric layers, so that P = 2 and C (1) = C (2) = C, bonded to an elastic substrate beam as shown in Figure 2. Thus, the harvester has a total of three layers (K = 3). The harvester is assumed to be a linear inputoutput system which is characterized by its Frequency Response Functions. The FRFs are defined here as the output responses of the harvester (displacement or voltage) per base acceleration versus the excitation frequency.

Transforming Equations ( 19) and ( 20) to the frequency domain, the equivalent expression for nodal displacements FRFs relative to the electrical configurations (load resistance mounted in series or parallel connection) is given in the following equation in which the electrical DOFs are eliminated [START_REF] Anton | Multifunctional Unmanned Aerial Vehicle Wing Spar for Low-Power Generation and Storage[END_REF][START_REF] Aloui | Sensitivity analysis of frequency response functions for load resistance of piezoelectric energy harvesters[END_REF].

U -ω 2 w b =   -ω 2 M m + jωC m + K m + jω K eq c K eqT c 1 R + jωC eq   -1 m * ( 26
)
where K eq c and C eq are respectively the equivalent electromechanical coupling vector and the equivalent capacitance for a series and parallel connections of the load resistance with piezoelectric layers [START_REF] Aloui | Sensitivity analysis of frequency response functions for load resistance of piezoelectric energy harvesters[END_REF], which are given in Table 1. For the mechanical vibration response, only the transverse tip dis-

placement FRF w n -ω 2 w b
is considered in this study. The index n = N/3 is the total number of nodes used to discretize the beam harvester with linear elements.

Table 1: Equivalent terms of a bimorph energy harvester for the series and the parallel connections of the piezoceramic layers

Terms Series connection Parallel connection

K eq c K (1) c + K (2) c 2 K (1) c + K (2) c C eq C 2 2C
Similarly, for the electrical outputs, the resultant voltage FRF of the harvester is obtained as a function of the nodal displacements FRFs as follows [START_REF] Anton | Multifunctional Unmanned Aerial Vehicle Wing Spar for Low-Power Generation and Storage[END_REF]:

V -ω 2 w b = jω K eqT c 1 R + jωC eq p U -ω 2 w b ( 27
)
The current I and the electrical power P FRFs are derived from the voltage FRF:

I -ω 2 w b = 1 R V -ω 2 w b ( 28
) P (-ω 2 w b ) 2 = 1 R V -ω 2 w b 2 (29)
The global finite element matrices appearing in the frequency response functions Equations ( 26) and ( 27) establish the dependence of the response of the system on a set of parameters representing the physical characteristics of the harvester structure. In general, we assume that any complex frequency response function H of the BPEH can be expressed in the following form:

H = H (ω, θ) ( 30
)
where θ is a vector of physical parameters of the harvester.

Global Sensitivity Analysis of the Frequency Response Functions

The Elementary Effects Method

The Elementary Effects Method constitutes a simple way for screening the effect on the output responses of a few important input factors among the many factors that can be contained in the studied model. The fundamental idea behind this method is given by Morris [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF][START_REF] Morris | Exploratory designs for computational experiments[END_REF], who introduced the concept of elementary effects, proposing the construction of two sensitivity measures: the mean and the standard deviation of the EEs for each input parameter. In practice, this method is conducted by: (a) defining the model and its input parameters, (b) assigning Probability Density Function (PDF)

to each input parameter, (c) generating an input matrix using an appropriate random sampling method, (d) calculating the corresponding output response vector, and (e) computing the EEs and then the sensitivity measures of each input/output relationship [START_REF] Hamby | A review of techniques for parameter sensitivity analysis of environmental models[END_REF]. The average and standard deviation are then computed for a set of elementary effects for each factor. The main goal is to determine which input factors could be considered to have effects which are [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF][START_REF] Sanchez | Application of sensitivity analysis in building energy simulations: Combining firstand second-order elementary effects methods[END_REF]:

1. Negligible (low mean, low standard deviation).

2. Linear and additive (high mean, low standard deviation).

Nonlinear or involved in interactions with other factors (high standard deviation).

This method computes the sensitivities of the outputs considering one factor at time among the input factors through the concept of EEs, which are approximations of the first order partial derivatives of the response. These elementary effects are estimated at a set of sampled points in the input space.

The finite element model of the bimorph piezoelectric energy harvester defined in Equations ( 26) and ( 27), is the model-parameters considered. The frequency response function H in Equation ( 30) is the model-output and the k independent parameters in vector θ are the model-inputs. Each parameter θ i is assumed to be scaled to take on values in the interval [0, 1] as follows:

θ i = θ i -θ i,min θ i,max -θ i,min (31)
where θ i,max and θ i,min are respectively the maximum and minimum bounds of the i th parameter. The bar is omitted below to simplifying the notation. Thus, these non-dimensional parameters are randomly selected in a k-dimensional unit hypercube.

A discretized approach of the input space is proposed. In fact, the possible input values will be restricted to an experimental design, which is a regular k-dimensional p-levels grid, where each θ i may take a value from the set {0, 1/(p -1), 2/(p -1), • • • , 1}. For a given vector θ, the elementary effect of the i th input parameter is defined as follows [START_REF] Arruda | Mechanical Joint Parameter Estimation Using Frequency Response Functions and Component Mode Synthesis[END_REF][START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF]:

EE i (θ) = |H(θ i + ∆, θ ∼i )| -|H(θ)| ∆ ( 32
)
where |.| is the modulus, θ i varies between 0 and (1 -∆); ∆ is the variation size which is a predetermined multiple of 1/(p -1), and θ ∼i is the set of all parameters except the i th parameter. The originality of the elementary effects method is based on selecting a set of r trajectories where parameters are changed one at a time on the design grid [START_REF] Campolongo | From screening to quantitative sensitivity analysis. A unified approach[END_REF]. A trajectory enables the calculation of an elementary effect for each input parameter i between two points of the trajectory (θ (j) and θ (j+1) , where j ∈ {1, .., k}) [START_REF] Sanchez | Application of sensitivity analysis in building energy simulations: Combining firstand second-order elementary effects methods[END_REF][START_REF]Global sensitivity analysis: the primer[END_REF].

Therefore, each trajectory requires (k + 1) model evaluations (simulations)

to calculate k elementary effects (one EE per parameter).

Figure 3 presents a trajectory in a 3-dimensional space which is composed of a four sampling points {θ (1) , • • • , θ (4) } and ∆ = 0.25. Along a trajectory each input parameter is increased or decreased by the same step ∆. Considering, for example, the third component θ 3 . This components differs only in the two consecutive sample points θ (3) and θ (4) , thus:

EE 3 = |H(θ 4 )| -|H(θ 3 )| 0.25 (33)
A set of r different trajectories (with index t, where t = 1, .., r) provides r estimates of elementary effects related to each input parameter i, at the cost of r × (k + 1) simulations.

The first measure of sensitivity is the average of the elementary effects which is computed for each input parameter i:

µ i = 1 r r t=1 EE * (t) i ( 34
)
where EE * (t) i is the normalized elementary effects of EE (t) i . This measure has been rectified by Campolongo et al. [START_REF] Campolongo | An effective screening design for sensitivity analysis of large models[END_REF], who has recommended the use of the mean of absolute elementary effects, defined in Equation [START_REF] Anton | Multifunctional Unmanned Aerial Vehicle Wing Spar for Low-Power Generation and Storage[END_REF], rather than the usual average, since some EEs can eliminate each other in non-monotonic models. Thus, the mean of the absolute values of EEs is:

µ * i = 1 r r t=1 |EE * (t) i | ( 35
)
For a given parameter i, if the r estimates of EEs have the same sign, that means that this parameter has a monotonic effect on the output response, increasing or decreasing, depending on the sign of the EEs. In this case, µ * i is equal to the absolute value of µ i . The reverse is also true; if

µ * i = abs(µ i ),
the effects of the i th parameter are monotonic.

The second proposed measure is the standard deviation of the EEs σ i , which assesses the extent of interactions and non-linear effects of each input parameter i. This measure is defined as follows:

σ i = 1 r -1 r t=1 (|EE * (t) i | -µ * i ) 2 (36)
A graphical representation of the (µ * i , σ i ) points in the (µ * , σ) plane allows a comparison of the sensitivities of different parameters taking into account at the same time the two sensitivity measures.

Latin Hypercube Sampling for Elementary Effects Method

Latin Hypercube Sampling (LHS) has become a popular technique for generating designs to computational experiments [START_REF] Andres | Sampling methods and sensitivity analysis for large parameter sets[END_REF][START_REF] Stein | Large sample properties of simulations using latin hypercube sampling[END_REF]. To do that for the elementary effects method, three steps are essential:

Numerical Example

This section presents a numerical example of the harvester introduced in Section 4, which is composed of two identical piezoceramic layers of PZT-5A bonded on the top and the bottom surfaces of an aluminum substrate.

The geometrical dimensions and material properties of the harvester used in this example are given in Table 2 (adapted from Erturk and Inman [START_REF] Erturk | Piezoelectric energy harvesting[END_REF]).

The ratio of the overhang length to the total thickness of the harvester is about 85.7, which makes it reasonable to neglect the shear deformation and the rotary inertia effects for the first few vibration modes. For the numerical simulation, we used modal mechanical damping ratios ζ 1 = 0.01 and ζ 2 = 0.012 for the first two modes. The proportional coefficients α and β given in Equation ( 18) are computed using these two damping ratios. The load resistance R is mounted in series with the piezoceramic layers as shown in Figure 2. For the finite element discretization, the harvester is modeled in 1D using 30 linear elements with three degrees of freedom per node.

The model is validated for the frequency range from 0 to 5000 Hz by comparing the first three resonance frequencies computed by the finite element model to those given by the analytical solution computed by Erturk et Inman [START_REF] Erturk | Piezoelectric energy harvesting[END_REF] for the same harvester for the short-circuit and open circuit conditions. The first three natural frequencies are listed in Table 3 with the effective electromechanical coupling factor k ef f defined in Equation [START_REF] Shu | An improved analysis of the SSHI interface in piezoelectric energy harvesting[END_REF] which is computed for the FE-model. We can notice that the natural fre-

quencies for the open-circuit condition are higher than those obtained in the short-circuited case due to the electromagnetical coupling effects and the added stiffness in the open-circuit condition. Therefore, the natural frequencies for the two approaches are in excellent agreement for all modes.

An infinitesimal error for the third mode can be due to the numerical computational errors caused by the approximation solution of the finite element method, which is so sensitive to the assigned discretization. These results confirm the validation of the finite element harvester model and allow for pursuing the study. . Mode (r) f cc r (FE) f cc r (Anal. [START_REF] Erturk | Piezoelectric energy harvesting[END_REF]) f co r (FE) f co r (Anal. [START_REF] Erturk | Piezoelectric energy harvesting[END_REF]) k ef f,r (FE) The selected uncertain parameters are the Young modulus and the material densities for both the substrate and the piezoelectric layers, also the premittivity and piezoelectric constants for the piezoceramic layers. Table 4 lists the tolerances provided by the manufacturers as percent varation from the nominal values [START_REF] Ruiz | Uncertainties propagation and global sensitivity analysis of the frequency response function of piezoelectric energy harvesters[END_REF]. A tolerance of ±20% is reported for the piezoceramic PZT-5A electromechanical properties and ±10% for the Young's modulus and density of the aluminum substrate material. An uniform distribution is assumed for all materials properties of the harvester. Parameters Bounds

E p ±20% E s ±10% ρ p ±20% ρ s ±10% e 31 ±20% 33

±20%

The application of the Morris method for the voltage FRF of the bimorph piezoelectric energy harvester consists in computing r = 10 elementary effects for each parameter listed in table 4. Thus, 70 simulations are required for a fixed excitation frequency and load resistance. The aim of the study is to identify the effect of the excitation frequency on the importance order of the parameters particularly in the vicinity of the first mode. Furthermore, the system is be studied for both a short-circuit condition simulated using a load resistance R = 100 Ω and an open-circuit condition simulated using for a load resistance R = 10 7 Ω. [START_REF] Larbi | Structural-Acoustic Vibration Reduction Using Switched Shunt Piezoelectric Patches: A Finite Element Analysis[END_REF] as the ratio of the exchanged mechanical energy and electrical energies is shown in Figure 10 for each of the 10 electromechanical characteristic sample vectors used in this GSA study. One can observe that the EEMCF can reach a maximum value of 0.34 for the 42 th sample vector. This results can be very useful for the research of optimal harvester parameters that maximize this factor.

Conclusion

In this paper, a finite element formulation of the electromechanical coupling problem for a laminated piezoelectric cantilever beam introduced. The finite element formulation is applied in the case of a symmetric bimorph piezoelectric energy harvester. The FRFs of the harvester are used as the performance predictors. The Morris global sensitivity method using the el- Future research may consider second order sensitivity analysis using the elementary effects method to minutely study the interactions between pa-

Figure 1 :

 1 Figure 1: A multilayer piezoelectric cantilever beam: k ∈ {1, • • • , K}

Figure 2 :

 2 Figure 2: A Clamped free bimorph piezoelectric energy harvester under base excitation with a load resistance mounted in series

Figure 3 :

 3 Figure 3: An example of a trajectory in the input domain for k = 3 and p = 5

( 1 )

 1 The LHS generates p = 2×r random samples of k uncorrelated parameters in the k-dimensional unit hypercube; where each component of θ has an uniform distribution in the interval [0, 1]. (2) Returning to the specified distribution of the uncertain parameters by inverting the Cumulative Distribution Function (CDF). (3)Applying the One factor at a time strategy, which consists in adding the r ×(k -1) points in order to establish r trajectories, each trajectory contains k + 1 points. Therefore, the computational experiments matrix is composed of r × (k + 1) random points of k uncorrelated parameters. Each row of the matrix is a point in the input space. The rows of the matrix are sorted in r blocks, each block includes k + 1 rows (each block presents a trajectory).Inside each block, points differ only in one component at a time. Thus, each block allows to compute one Elementary Effect per input parameter.

Figure 4 Figure 4 :

 44 Figure 4 shows an example of two parameters uniformly distributed in [0 1] 2 space. In order to determine 4 elementary effects per parameter, the LHS method generates 8 random points given in Figure 4.a. Adding 4 points whose one parameter changes at the time, provides 4 trajectories as shown in Figure 4.b. The LHS with the one factor at the time strategy ensure two levels per input parameter for one trajectory of the EEs method.

Figure 5 Figure 5 :Figure 6 33 Figure 6 :

 556336 Figure 5 shows the voltage and tip displacement FRFs responses of the harvester versus the excitation frequencies in the vicinity of the first mode for 70 samples. One can be observe the variability in the first natural frequency and in amplitude peaks. The voltage FRF peak varies between 0.02032 and 0.01152 V /g whereas the tip displacement FRF peak varies between 442.46

Figure 7 33 Figure 7 :Figure 8 :Figure 9 Figure 9 :

 7337899 Figure7shows the (µ * , σ) scatter graph of the FRF voltage for the six parameters to the variation of the electromechanical characteristics for a load resistance R = 100 Ω and and excitation frequency equal 185 Hz. The limits of the boxes correspond to a confidence interval of 95%. Clearly, we can see that the sensitivities of the absolute mean and standard deviation of the piezoelectric density ρ p and the elastic modulus E p are relatively significant compared to those of the piezoelectric constant, permittivity constant and substrate properties. Therefore, the density and elastic modulus of piezoelectric layers are the most important parameters that affect the output voltage FRF. Furthermore, the average and the standard deviation are of the same order of values. This means that the elastic modulus and

Figure 10 :

 10 Figure 10: The EEMCF k ef f for the 70 simulations used in the GSA

Table 2 :

 2 Physical parameters values of the harvester structure of Fig.2

	Parameters Descriptions	Units 1 PZT-5A Aluminum
	L	Beam length	mm	30	30
	b	Beam width	mm	5	5
	h p , h s	Layers thickness	mm	0.15	0.05
	E p , E s	Young's modulus	GP a	61	70
	ρ p , ρ s	Mass density	kg/m 3	7750	2700
	e 31	Piezoelectric constant C/m 2	-10.4	-
	33	Permittivity constant nF/m	13.3	-

Table 3 :

 3 The first three natural frequencies for the short circuit and open circuit conditions computed from the Finite Element Method (FE) and the Analytical solution (Anal.)[START_REF] Erturk | Piezoelectric energy harvesting[END_REF]

Table 4 :

 4 Bounds for the uniform PDF distribution associated to the electromechanical parameters; bounds are expressed in terms of the nominal values[START_REF] Ruiz | Uncertainties propagation and global sensitivity analysis of the frequency response function of piezoelectric energy harvesters[END_REF]

Unit of parameter (PU), it will be used to define sensitivities units

rameters properties of the harvester.