Fabio Candotto Carniel

Davide Gorelli

Emmanuel Flahaut

Lorenzo Fortuna

Cecilia Del Casino

Giampiero Cai

Massimo Nepi

Maurizio Prato

Cecilia Del Casino

Maurizio Prato De

Mauro Tretiach

Graphene oxide impairs the pollen performance of Nicotiana tabacum and Corylus avellana

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

Introduction

Graphene is "a two-dimensional crystal composed of monolayers of carbon atoms arranged in a honeycombed network with six-membered rings" 1 and is the first 2D material made available. 2 Graphene based materials (GBMs) have attracted wide interest in the last decade because they display extraordinary physical and chemical properties, such as supreme electrical and thermal conductivities and high surface area, clearly distinguishing them from graphite and carbon nanotubes. Therefore, GBMs are expected to be used more often in emerging applications in sectors such as optoelectronics, nanocomposites, biomedicine, energy generation and storage assessing the impact of a wide variety of chemicals on plant metabolism. 20 It was also used as an indicator of air pollution because the pollen performance, when described in terms of pollen germination and pollen tube elongation, can be affected (depending on species) by gaseous airborne pollutants such as NO x , O 3 and SO 2 . [21][22][23] More recently, it was demonstrated that also airborne particulate matter affects pollen performance: Speranza et al. 24,25 found that Ag-and Pd-nanoparticles impair the in vitro pollen germination and the elongation of the pollen tube in the kiwifruit vine, Actinidia deliciosa C. F. Liang & A. R. Ferguson. Interestingly, this effect is derived from the nanoparticles themselves rather than from the ions carried on their surface. 24,25 Given these premises, our hypothesis is that the interaction of aerodispersed GBM flakes with pollen might impair the pollen performances affecting negatively the reproductive process of seed plants. To verify this, the effect of GO, the most reactive GBM tested so far, was verified on the pollen performance of two species, the entomophilous tobacco plant (Nicotiana tabacum L.) and the anemophilous common hazel (Corylus avellana L.), through in vitro germination experiments.

Materials and methods

Graphene oxide preparation and characterization

Graphene oxide (GO) was provided by Grupo Antolin (Burgos, S) within the Graphene Flagship project (http://grapheneflagship.eu/) and it was produced by oxidation of carbon fibers (GANF Helical-Ribbon Carbon Nanofibres, GANF®) in sulfuric acid with sodium nitrate at 0 °C. Dry GO was resuspended in water and sonicated in an ultrasonic bath prior to the experiments.

The surface chemical composition of GO was analyzed by X-ray photoelectron spectroscopy (XPS) using a Thermoelectron K-α apparatus (ThermoFisher Scientific, MA, USA). The photoelectron emission spectra were recorded using Al-Kα radiation (hυ = 1486.6 eV) from a monochromatized source. The X-ray spot diameter on the sample surface was 400 μm. The pass energy was fixed at 30 eV. The spectrometer energy calibration was performed using the Au 4f 7/2 (83.9 ± 0.1 eV) and Cu 2p 3/2 (932.8 ± 0.1 eV) photoelectron lines. The background signal was removed using the Shirley method. 26 Atomic concentrations were determined from photoelectron peak areas using the atomic sensitivity factors reported by Scofield, 27 taking into account the transmission function of the analyzer. The morphology of the dried particles and their sizes were characterized using a JEM 1400 (Jeol USA, Inc., MA, USA) Transmission Electronic Microscope (TEM) operated at 120 kV, and the GO samples for the observations were prepared according to Garacci et al. 28 The GO average flake thickness was assessed by X-ray powder diffraction (XRD) analysis using a D4 ENDEAVOR diffractometer (Bruker Corp., MA, USA) at a wavelength of Cu Kα1 (λ = 0.15418 nm). The GO structure was assessed by Raman spectroscopy using a Labram-HR800 (Horiba Group, J) equipped with a thermoelectrically cooled CCD, using and high sensitivity sensors. 2 The estimated annual production of graphene is around 120 tons 3 but is supposed to increase to 1000 tons by 2019. 4 Estimates for the GBM market foresee investments of almost $400 million by 2025, [START_REF] Ghaffarzadeh | Graphene 2D Materials and Carbon Nanotubes: Markets Technologies and Opportunities 2016-2026[END_REF] and industrial patent applications have promptly increased in recent years. 4 Hence, an important part of the applications exploiting GBMs will likely end up in devices and objects of common use increasing the GBM demand. In this regard, GBM inclusion in the rubber of bike tires is already a reality (https://www.vittoria.com/graphene-evidence/), [START_REF] Vittoria | Graphene, Join the revolution[END_REF] whereas the inclusion in tires for motor vehicles is a patented application. [START_REF] Zhamu | [END_REF] Furthermore, new graphene oxide-Ag composites are being developed to fight crop plant diseases, such as the bacterial leaf blight (Xanthomonas oryzae pv. oryzae) affecting rice, 8 through aerial applications. The large scale production of GBMs, wearing out of graphene-enriched products, and poor disposal of the derived wastes might result in a significant release and accumulation of GBMs in the environment, as has already happened for other synthetic materials such as plastics. 9,10 Thanks to their dimensions from hundreds of nanometers to a few micrometers and to their extreme lightweight, GBM flakes or nanoparticles might be aero-dispersed as particulate matter with a diameter size ≤2.5 μm (PM 2.5 ). Once aero-dispersed, PM 2.5 can be transported for long distances, as reported for carbon black nano-particles, 11 and eventually settle over the vegetation as wet or dry depositions. These premises raise important concerns about possible negative impacts of GBM release on the biota.

So far, widely varying effects (including positive and negative) of GBMs on seed plants have been reported, possibly owing to different experimental conditions (materials, concentrations, exposure time, protocols etc.) and/or species tested. [12][13][14][15] In general, however, there is agreement on the greater toxicity of graphene oxide (GO) with respect to other GBMs. 16 This is likely due to the high number of oxygen functional groups, such as hydroxyl, carboxyl and epoxy, spread over the graphene lattice which render GO more hydrophilic and reactive towards the surrounding environment. 17 Expectedly, the effects of GBMs on seed plants were always evaluated on the vegetative body at different developmental stages, from seed germination to the flowering adult plants. However, nothing is known yet about the effects of GBMs on the complex, delicate reproductive process, which strongly relies on the transfer and interaction of the male pollen grains with the female reproductive structures. The few-celled pollen grain plays the crucial role of delivering the sperm cells, through the formation of a pollen tube, to the egg cell(s) inside the ovule. If this task is successful, fertilization can occur, with the formation of a new embryo, and later, of a seed. This process is fundamental for the reproduction of almost all seed plants but it is also important for humankind since the yield of crop species, largely consisting of seeds and fruits, relies in toto on it. Pollen is more sensitive to environmental pollutants compared to most other structures of the vegetative body of seed plants. 18,19 For this reason, pollen is considered as an experimental model for avellana using a Zeiss Axiocam MRm camera connected to a Zeiss Axiophot microscope (Zeiss, Oberkochen, D). The pollen performance, i.e. germination percentage and pollen tube length, was measured on the control and treated samples. The germination percentage was calculated by scoring at least 150 randomly selected grains which were considered germinated if the pollen tube was longer than the average diameter of the pollen grains. The pollen tube length was measured for at least 50 randomly selected pollen tubes. Further measurements (n > 60) of the angle made by the pollen tube were taken from the control samples and from the samples treated with GO at 100 μg mL -1 .

Viability assay

The viability of N. tabacum pollen grains was evaluated in the control and in the treated samples (100 μg mL -1 GO only) using the fluorochromatic reaction test. 30 A few drops of fluorescein diacetate solution (2 mg mL -1 FDA in acetone) were added to glass slides, waiting 30-40 seconds to allow acetone evaporation. 23 Thereafter, 30 μL of the pollen culture was added to each slide and then covered with coverslips. FDA fluorescence was visualized after 30 min, 1 h and 2 h using a Zeiss Axiophot fluorescence microscope equipped with a Zeiss Axiocam MRn camera; evaluation of viability percentage was done on the captured images by counting at least 150 pollen grains for each sample.

Analysis of pollen tube intracellular pH

The intracellular pH variation along the N. tabacum pollen tubes was assessed in the control samples and in the samples treated with GO at 100 μg mL -1 using the fluorescent dye BCECF-AM (acetoxymethyl). A BCECF-AM (2,7-bis-(2carboxyethyl)-5-(e-6)-carboxyfluorescein) ester probe was used to visualize proton levels (i.e. pH) in the pollen tubes of N. tabacum. 31 A final concentration of 5 μM was obtained from a 1 mM stock solution in dimethyl sulfoxide (DMSO); the required volume was directly added to the germination medium containing resuspended pollen grains. The cytosolic pH was immediately measured after addition of the probe to prevent probe absorption of the organelles (vacuoles). The samples were observed using a Zeiss AxioImager fluorescence microscope equipped with structured illumination in the FITC filter (Ex. λ = 488 nm; Em. λ = 515 nm).

Pollen performances at different pH in relation to GO treatments

The pH of BK suc and BK PEG with or without GO at 100 μg mL -1 was measured using a Crison Basic 20 (Crison, Hach Lange, S) pH-meter in continuously stirred media until the variation measured by the probe was less than ±1 mV per minute. The germination performance of the C. avellana pollen was evaluated in the pH range between 3 and 7 by incubating the pollen for 3 h in the BK medium adjusted with 0.1 M NaOH or H 2 SO 4 , as appropriate. Subsequently, the pollen germination performance was further investigated in the excitation lasers at 633 (red laser, He/Ne) and 532 nm (green laser, diode-pumped solid-state laser). Five spectra were averaged for each sample and the intensity ratio between the D and the G lines (I D /I G ) was measured after baseline correction.

Two additional GOs were purchased from Graphenea (San Sebastián, E) and Sigma-Aldrich (St. Louis, MO, USA); they were provided as water suspensions at concentrations of 4 and 2 mg GO mL -1 , respectively. Data regarding these materials can be found in the safety datasheets available on the websites of the producers.

Pollen material

The pollen of Nicotiana tabacum was collected in the summer of 2016 from dehiscent anthers of plants grown in the experimental fields of the Botanical Garden, University of Siena. Mature undehiscent anthers were sampled, transported to the laboratory and kept at 20 °C until dehiscence. The pollen was gently separated from the anthers and collected in Eppendorf tubes. The pollen of Corylus avellana was collected in January 2017 from 15 native trees in the Classical Karst (Trieste, NE Italy) far from pollution sources. Twigs bearing unripe catkins were sampled, immediately transported to the laboratory, and their bases cut under water and kept immersed under laboratory conditions (20 °C, dim light) until flower ripening and stamen dehiscence (c. 48 h). The harvested pollen of both species was sieved through 100 and then 60 μm mesh sizes to remove debris, dehydrated over silica-gel (RH ∼5%) for 48 h and stored at -20 °C until use.

Germination of N. tabacum and C. avellana pollen

Prior to use, the pollen of N. tabacum and C. avellana was thawed and then rehydrated in a moist chamber for 2 and 3 h, respectively. Germination experiments were performed by suspending the pollen at a concentration of 1 mg mL -1 in Brewbaker & Kwack's culture medium 29 To assess the effects of GO, the pollen of both species was germinated in BK without GO (control samples), and in BK supplemented with GO suspensions to final concentrations of 25, 50 or 100 μg mL -1 (treated samples). Three to four replicates have been prepared for both the control and treated samples. Germination experiments were performed at a controlled temperature of 25 °C. An additional experiment was performed on the N. tabacum pollen replacing sucrose with polyethylene glycol (PEG) to a final concentration of 13% (BK PEG ).

Assessment of pollen performances

Images of the germinating pollen were taken after 1 and 2 h of pollen tube growth for N. tabacum and 3 and 5.5 h for C.

(1 and 2 h for N. tabacum; 3 and 5.5. h for C. avellana) and GO concentrations (0, 25, 50 and 100 μg mL -1 ) as the categorical factors. The same statistical analysis was performed for the pH experiment using germination time and pH as the categorical factors. The effect of GO on the pollen tube angle and on the pollen viability was assessed by running a Mann-Whitney U-test for independent samples. All calculations were performed using Microsoft Office Excel 2003 SP3 (Microsoft corporation, WA, USA) and STATISTICA 8.0 (StatSoft Inc., OK, USA).

Results and discussion

Characterization of GO

In the case of the GO sample provided by Grupo Antolin, the elemental analysis data (Table S1 †) obtained either by classical chemical analysis (data as provided by Grupo Antolin) or measured experimentally by XPS at the CIRIMAT facility revealed some interesting information. In particular, the C/O ratio, which is often reported in publications, is very sensitive to the input data (weight concentrations or atomic concentrations). If XPS analysis is focused only on the C and O elements, neglecting other elements such as hydrogen which is essentially impossible to detect by this technique, the C/O following media: i) BK suc (control), ii) BK suc adjusted with H 2 SO 4 to pH 4.2, iii) BK suc with GO at 100 μg mL -1 and iv) the latter adjusted with NaOH to pH 6.3. For comparison, the pH of water suspensions (100 μg mL -1 ) of commercially available GOs, i.e. GO from Graphenea (San Sebastián, E), and Sigma-Aldrich (St. Louis, MO, USA), was also measured.

Data analysis

Measurements of pollen performance, tube angles and fluorescence signals from FDA and BCECF-AM were performed using the ImageJ software ver. 1.51J8 (Wayne Rasband, NIH, USA) after calibration of the images with the scalebar command of the Axiovision software (Zeiss, D).

In order to exclude bias due to slow growing pollen tubes or to pollen which stopped growing, only the 30% longest tubes of both the control and treated samples were selected for statistical analyses. The length of selected pollen tubes was expressed as the percentage ratio with respect to the mean tube length measured in the control samples. Factorial ANOVA followed by Tukey's HSD post-hoc test were used to evaluate the statistical differences among the control and treated samples in terms of both pollen germination percentage and pollen tube length using optimal germination times the same conditions (Fig. 2D). In any case, GO flakes were found adhering to the pollen grains and the pollen tube surfaces in both species (Fig. 2D andE). The pollen grains of N. tabacum are generally coated by a lipoid "pollen coat" or "pollenkitt" which can account for 10-15% of the total pollen mass. 32,33 This coat is typical of entomophilous species: it has adhesive properties allowing pollen grains to stick to insects and to the stigmatic surfaces of flowers. 34 Moreover, the germinating pollen grains of N. tabacum release lectins and polylectins, which may also cause grain agglutination. 33 Both may explain the pollen-GO aggregation and why this did not occur in the anemophilous C. avellana. As reported by Chen et al., 35 PEG has a good affinity to lipoid substances, and therefore it might have partially removed the pollen coat reducing the aggregation between the pollen and GO.

Effect of GO on pollen performances and viability

The pollen of N. tabacum and C. avellana germinated in BK suc (control samples) had normal tube extrusion and elongation (Fig. 2A) with final germinations of 75% and 81% (Fig. 3A and4A), respectively, suggesting that the pollen materials were healthy at the beginning of the experiments. The pollen performance of both species was significantly affected by GO (see Tables S2 and S3 †). The pollen germination and tube length of N. tabacum decreased to 55% and 81% with respect to the control values (Fig. 3B andD; P < 0.05), respectively, after 2 h of germination at 100 μg GO mL -1 . The decrease in pollen performance was greater in C. avellana; the germination remained steady ratio is largely over-evaluated (2.3) while the use of the weight concentrations obtained from the classical chemical elemental analysis leads to a C/O ratio close to 1. This is important to know because in most publications the chemical elemental analysis data are absent; so, it is impossible to know which elements are also present in addition to C and O. In this case, it is obvious that the presence of a high atomic % of hydrogen leads to important differences in terms of C/O ratio depending on considering it or not.

The GO nanoparticles (Fig. 1C) had a thickness of 3.9 nm, calculated using the Scherrer equation applied to the (002) reflection (Fig. S1 †) and assuming an interlayer spacing of 0.8 nm obtained from the position of the same reflection (peak at 11.05°o f Fig. S1 †). These data suggest that the estimated number of layers was between 4 and 5 per particle. The intensity ratios between the D and the G bands (Fig. S2 a and b †) measured by Raman spectroscopy were 1.03 and 1.17 at 532 and 633 nm, respectively, revealing in both cases a rather disordered material.

Physical interaction between GO and the pollen

The Nicotiana tabacum pollen germinated in BK suc with GO formed clusters intermingled with GO flakes especially at the highest concentration (100 μg mL -1 ; Fig. 2A andB), which made the measurements of pollen performance unreliable (n. a. data in Fig. 3A andC). Interestingly, the replacement of sucrose with PEG reduced the aggregation effect allowing the measurements (Fig. 2B and C and 3B and D). Such pollen-GO clusters were not observed in C. avellana germinated under until 50 μg GO mL -1 , then it decreased to 12% and 14% at 100 μg mL -1 after 3 and 5.5 h, respectively (Fig. 4A). Differently, the pollen tube length decreased progressively with the increase in GO concentration, reaching 62% and 43% of the control value at 100 μg GO mL -1 after 3 and 5.5 h, respectively (Fig. 4B). Furthermore, although the viability of the N. tabacum pollen was not affected at 100 μg GO mL -1 (Fig. 5), the percentage of pollen tubes showing a bent growth doubled with respect to that of the control samples (Fig. 6; P < 0.05). A similar trend was observed in the pollen tubes of C. avellana, but the difference between the control and treated samples was not statistically significant (Fig. 6). The fact that GO had an effect on the pollen performance of both species, although with a different intensity, agrees well with the literature reporting the toxic effects of GO on phylogenetically unrelated model organisms. 16 However, mechanisms of toxicity are still unclear: some authors reported that GO can mechanically damage cells or tissues with its sharp edges 36 and reduce metabolic activity and viability when it is internalized into cells. 37,38 The great majority, though, agrees that GO causes oxidative stress inducing the formation of reactive oxygen species. [39][40][41] In our case, the observed effects are mostly related to the GO composition in terms of functional groups and their interaction with the surrounding media, as it will be discussed below.

Effect of GO on the pH of germination media

Graphene oxide is characterized by a high number of functional groups containing oxygen, with hydroxyl and epoxy groups spread over the graphene lattice, and carboxyl groups at the flake edges. 17 Carboxyl groups are weak acids and together with hydroxyl groups deprotonate in water solutions lowering the pH. Furthermore, GO could still have some residual sulfuric acid (as suggested by the elemental composition in Table S1 †), which is used in the modified Hummers' method for the oxidation process of the graphite flakes. 42 For these reasons, the pH of the BK media and water suspensions of GO of different origins was investigated (Table 1). At 100 μg GO mL -1 , the pH of BK suc and BK PEG decreased by 2.00 and 0.97 points, respectively (Table 1). Water suspensions of GO, GO washed to remove potential impurities from the Hummer's method, GO from Graphenea and GO from Sigma-Aldrich had pH values below 4.25 (Table 1), confirming (a) the strong acidic property of all tested GOs, and (b) the negligible sulphuric acid contamination, because the washing of GO did not change the final pH value.

Effect of pH on pollen performances

Considering the strong acidification effect of GO on the BK suc media, we further investigated the effect of pH on the pollen performance of N. tabacum and C. avellana. As reported by Tupí and Říhová, 43 the N. tabacum pollen tube growth has its optimum in the pH range of 5-6.5 and decreases at lower pH (Fig. 7A). The Corylus avellana pollen had the best performance at pH 7 with 84% of germinated pollen and an aver-age tube length of 100.4 μm after 3 h of germination (Fig. 7B). Below pH 7, pollen germination slightly decreased to 78% at pH 5-6, fell to 40% at pH 4 and zeroed at pH 3 (Fig. 7B). Differently, the pollen tube length progressively decreased to 73.9 μm at pH 6 and then to 57.9 μm at pH 4 (Fig. 7B). These results help to clarify the contribution of different "toxicity" mechanisms, direct or not, on the performance of N. tabacum and C. avellana pollen germinated in GO-enriched media (Fig. 3). In the first case, PEG decreased the GO-dependent acidification stabilizing the pH to 5.9, which falls in the pH range where the N. tabacum pollen performance is purportedly the best (see Table 1 and Fig. 7A). For this reason, the GO-induced pH change does not explain the decrease in pollen performance detected in N. tabacum, Table 1 pH of the BK media (BK suc , BK PEG ) without or with GO (Grupo Antolìn) and pH of dH 2 O with GO, GO washed in dH 2 O to remove residues of the production method (GO w ), GO commercialized by Graphenea (GO § ) and Sigma-Aldrich (GO # ) at a concentration of 100 μg mL elles 46 and vesicle trafficking needed to deposit new pectin at the tube tip for its growth. 47,48 For this reason, localized pH decreases could also explain the higher percentage of bent pollen tubes observed in N. tabacum (Fig. 6).

In the case of C. avellana instead, GO lowered the pH of BK suc to 4.35, i.e. the pH value at which the germination performances of C. avellana are the lowest (see Table 1 and Fig. 7B). To verify if the GO-dependent acidification was the main reason for the observed effects, a further experiment was performed, reversing the pH conditions of the media: those of the controls were acidified to pH 4.2, whereas those with GO were adjusted to pH 6.3. The acidified controls decreased similarly to GO-treated samples after both 3 and 5.5 h (Fig. 9) and, conversely, the adjusted GO-treated samples had a pollen performance similar to that of the pristine controls (Fig. 9). These results confirm that GO-induced acidification of the media is the main toxicity mechanism (for P = 0.000, see Table S4 †), causing the impairment of pollen performance in C. avellana.

The pH-decrease induced by high concentrations of GO might impair the pollen performance of a broad number of seed plants. For instance, Cox 49 tested the pollen performances of 13 forest species over a pH range from 2.6 to 5.6, suggesting a different type of interaction between pollen grains and tubes. Inhibition of pollen germination induced by Ag-nanoparticles has already been documented by Speranza et al., 25 who hypothesized a possible physical blockage of the transmembrane ion channels at the site of tube protrusion by the nanoparticles, preventing the formation of the inward ion current preceding pollen germination. 44 Alternatively, GO internalization at the site of the tube protrusion could have impaired the metabolic pathways triggering pollen germination. Further experiments on GBM internalization and GO-induced ROS production are currently ongoing to clarify the mechanisms on the basis of this germination inhibition. The reduction of the pollen tube length by GO can instead be explained, thanks to microscopy observations of GO treated samples stained with BCECF-AM, a fluorescent dye used to detect intracellular changes in pH. 31 The treated pollen suffered a pH variation along their pollen tubes. Fig. 8C andD show a pollen tube covered by GO flakes with regions emitting fluorescence at an intensity of up to two fold higher than that detected in the control samples (Fig. 8B vs. E). This certainly affects tube growth and tube directionality because low intracellular pH impairs the cytoskeleton 45 which in pollen tubes coordinates the movements of organsimulating the effect of acid rains on the reproductive process of seed plants. He found that in half of the species, the pollen performance was strongly reduced at pH 3.6.

Conclusions and environmental relevance

Aero-dispersed GBM flakes and their environmental fate is still an under-studied field. However, graphene flakes were detected in the air of small-scale graphene manufacturing workplaces (a laboratory and a research institution) in which careful measures to protect operators were adopted. 50 This demonstrates the high volatility of GBMs and raises concerns on the possible release of GBMs that will occur inside and around large scale-production facilities or when these materials are used in degradable objects. For these reasons, in this work, we started addressing the potential effects of GO, one of the most studied and potentially toxic GBMs, on the reproductive process of seed plants by in vitro experiments of pollen germination hypothesizing important future releases of GBMs from both point and widespread sources. We showed that GO impaired the pollen performance in N. tabacum and C. avellana at the highest concentration tested and most of this is derived from the acidic properties of the material. This could be relevant from an environmental point of view because strong depositions of GO flakes on the stigmatic surfaces of angiosperms could decrease pH to a level overcoming the buffering capacities of the organ. Pollen landing on an acidified stigma could be affected and consequently have reduced chances to successfully conclude the reproductive process, as already showed by Cox 49 in connection with the devastating phenomena of acid rains. 51 Furthermore, our experiments brought evidence of another effect of GO, potentially not related to the acidic properties but still important since it reduced the germination performance of N. tabacum by 20%, and doubled the frequency of bent tubes. Still unknown in its nature, this might be related to intracellular ROS development or to a perturbation of Ca 2+ availability, since this ion plays a fundamental role in pollen tube growth. 44 There are no doubts that more experiments involving other GBMs and seed plants, especially of economic value such as cereals, need to be conducted to clarify the mechanisms inducing this effect and its selectivity towards diverse species. It has to be noted though, that the scenario depicted here could be far more complex. Hypothetically, GO flakes deposited on the receptive surface could also affect the reproductive process in other ways, such as reducing the capacity of the stigma to intercept pollen in anemophilous species or hindering physically or chemically the signalling occurring between the stigma and the pollen grain. Also, sexual reproduction of entomophilous seed plants relies on insects for pollination which could collect GO flakes from flower organs or even eat them if the flakes land over nectar or water drops or over the pollen in the open anthers, entering de facto in other trophic chains. Furthermore, the effect of GO and other airborne metallic nanoparticles could be amplified by synergistic effects such as those reported by Tang et al. 52 in the cyanobacterium Microcystis aeruginosa.

Fig. 1

 1 Fig. 1 XPS analysis of GO showing details of the C1s (A) and O1s (B) bands. A typical TEM image of the GO material is shown in (C).

Fig. 2

 2 Fig. 2 Pollen of N. tabacum in BK PEG after 1 h (A), BK suc (B) or BK PEG (C) at 100 μg GO mL -1 after 2 h. Pollen of C. avellana in BK suc at 50 μg GO mL -1 after 5.5 h (D). Straight (E) or bent (F and G) grown pollen tubes of N. tabacum at 100 μg GO mL -1 . Bars: A, G = 20 μm; D = 50 μm; B, C, E, F = 100 μm.

Fig. 3

 3 Fig. 3 Germination percentage (A and B) and pollen tube growth (C and D) measured in Nicotiana tabacum. Pollen germinated in BK suc (A and C) or BK PEG (B and D) without (CTRL) or with 25, 50 or 100 μg GO mL -1 after 1 h (grey bars) and 2 h (dark grey bars). Values are reported as means ±1 st. dev. (n = 150 for A and C; n ≥ 50 for C and D). Statistically different groups are marked with different letters (two-way ANOVA, Tukey's HSD post-hoc test). n.a. = not available data.

Fig. 4

 4 Fig. 4 Germination percentage (A) and pollen tube length (B) measured in Corylus avellana. Pollen germinated in BK suc without (CTRL) or with 25, 50 or 100 μg GO mL -1 after 3 h (grey bars) and 5.5 h (dark grey bars). Values are reported as means ± 1 st. dev (n = 150 for A; n ≥ 50 for B). Statistically different groups are marked with different letters (two way ANOVA, Tukey's HSD post-hoc test).

Fig. 5

 5 Fig. 5 Percentage of N. tabacum viable pollen grains (n = 150) in the control (CTRL) and treated samples with GO at 100 μg mL -1 (GO) after 2 h of germination.

Fig. 6

 6 Fig. 6 Percentage of bent pollen tubes in the control (CTRL) and treated samples (100 μg mL -1 GO) of N. tabacum and C. avellana germinated for 2 h and 5.5 h, respectively. Values are means ± 1 st. dev. (n > 60); statistical differences between CTRL and GO are marked with an asterisk for P < 0.05 (Mann Whitney U-test).

18 Fig. 7

 187 Fig. 7 Percentage of pollen tube mass development in N. tabacum [redrawn from Tupí and Říhová (1984)] (A) and germination percentage (-○-) and pollen tube length (-▲-) of Corylus avellana (B) as a function of pH. Values in B are means ± 1 st. dev (n = 450 in A; n = 90 in B).

Fig. 8

 8 Fig. 8 Variation of the pollen tube pH in germinated Nicotiana tabacum L. pollen without or with GO at 100 μg mL -1 . The pollen tube from a CTRL sample stained with BCECF-AM and observed by fluorescence microscopy after 2 h of germination (A) and the fluorescence profile of the internalized dye measured along the same pollen tube (B). A pollen tube stained with BCECF-AM from a sample treated with GO and observed by fluorescence (C) and light (D) microscopy after 2 h of germination, and the fluorescence profile of BCECF measured along the same pollen tube of figure C (E) and in the pollen tube of samples treated with GO for 1 h (F). The arrows in C and D indicate GO flakes adhering to the pollen tube, whereas the arrowheads in C and E indicate the maximum of BCECF fluorescence in reference to the distance from the pollen tube apex, respectively. Bars = 10 μm.

Fig. 9

 9 Fig. 9 Germination percentage (A) and pollen tube growth (B) in Corylus avellana in BK suc (CTRL) or in BK suc acidified to pH 4.2 (CTRL A ) and in BK suc with GO at 100 μg mL -1 (GO100) or in GO100 neutralized to pH 6.3 (GO100 N ). Measurements were taken after 3 (grey bars) and 5.5 h (dark grey bars). Values are reported as means ± 1 st. dev. (n = 450 in A; n = 90 in B). Statistically different groups are marked with different letters for P < 0.05 (two-way ANOVA, Tukey's HSD post-hoc test).

  (BK) containing 1.62 mM H 3 BO 3 , 1.25 mM CaIJN0 3 ) 2 •4H 2 O (or CaCl 2 •2H 2 O), 2.97 mM KNO 3 (or KCl) and 1.65 mM MgSO 4 •7H 2 O. Germination of N. tabacum and C. avellana pollen was carried out in BK supplemented with sucrose (BK suc ) to a final concentration of 12% and 15% (w/v), respectively.

Acknowledgements

We gratefully acknowledge financial support from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement No. 696656 Graphene Flagship CORE 1. The authors thank Prof. Johanna Wagner (University of Innsbruck) for training FCC in germinating pollen in vitro, Dr. Alice Montagner for her valuable help in graphene handling and Drs. Susanna Bosi and Zois Syrgiannis (University of Trieste) for their technical advices. We are indebted to Ilaria Corsi (University of Siena) for the discussion and suggestions about planning of the experimental procedures. The personnel of the Botanical Garden of the University of Siena are also acknowledged for cultivating the plants of Nicotiana tabacum. Prof. Ester Vázquez Fernández-Pacheco is thanked for providing washed and non-washed GO materials. P. Lonchambon and G. Chimowa are acknowledged for their help with the physico-chemical characterization of the Antolin GO material.

Conflicts of interest

The authors declare no competing financial interests.


Jean-Frédéric Gerbeau

Damiano Lombardi

Eliott Tixier

How to choose biomarkers in view of parameter estimation

Keywords: feature selection, sparse optimization, inverse problems, electrophysiology, hemodynamics

In numerous applications in biophysics, physiology and medicine, the system of interest is studied by monitoring quantities, called biomarkers, extracted from measurements. These biomarkers convey some information about relevant hidden quantities, which can be seen as parameters of an underlying model. In this paper we propose a strategy to automatically design biomarkers to estimate a given parameter. Such biomarkers are chosen as the solution of a sparse optimization problem given a user-supplied dictionary of candidate features. The method is in particular illustrated with two realistic applications, one in electrophysiology and the other in hemodynamics. In both cases, our algorithm provides composite biomarkers which improve the parameter estimation problem.

Introduction

In numerous applications in biophysics, physiology and medicine, the system of interest is studied by monitoring a number of quantities, called biomarkers. For example in electrophysiology, an action potential can be viewed through simple quantities like the action potential duration, the amplitude, the rate of depolarization, etc. In hemodynamics, the systolic and diastolic pressures, or the pulse wave velocity, are typical biomarkers extracted from pressure measurements.

The biomarkers are obtained by applying a nonlinear map to the signal measured during experiments or clinical observations. They convey some information on hidden quantities, that are not directly measured. For example in hemodynamics, the pulse wave velocity can be linked to the arterial stiffness. When performing parameter estimation for a biophysical model, it is often much more convenient to work with biomarkers than with the whole signals. A natural question is therefore: which biomarker should be chosen to estimate a given parameter? Biomarkers are usually proposed by the community, based on physical intuition and experimental observations. They are often relevant in qualitatively describing the hidden quantities. However, in most practical applications, although the biomarkers exhibit a good correlation with respect to the hidden quantity they are designed to monitor, they have a non-negligible correlation with respect to others, making them less robust or of difficult interpretation.

In the present work, we propose a strategy to automatically design biomarkers. The basic ideas of our approach are: [START_REF] Gul | Parametric uncertainty and global sensitivity analysis in a model of the carotid bifurcation: Identification and ranking of most sensitive model parameters[END_REF] to design composite biomarkers that are maximally correlated with the hidden quantities they have to reveal, and minimally correlated with respect to all the others; (2) to provide a set of quantities making the parameter estimation better conditioned.

The biomarker design problem may be interpreted as a feature selection problem. Most of the literature considers the problem of selecting features in the input space in order to predict a given output (that may be the output of a computational model). Even though the aim of this work is reverse we will momentarily consider, for the sake of comparison, the biomarkers to be inputs and the parameters of interests to be the outputs. A common strategy to select a subset of the available features is by ranking or eliminating them according to a given criterion or score. This score may be based on a sensitivity analysis (e.g. first-order sensitivity indices [START_REF] Gul | Parametric uncertainty and global sensitivity analysis in a model of the carotid bifurcation: Identification and ranking of most sensitive model parameters[END_REF]), based on information theory (e.g. Fisher information matrix [START_REF] Cintrón-Arias | A sensitivity matrix based methodology for inverse problem formulation[END_REF] or mutual information [START_REF] Krier | Feature scoring by mutual information for classification of mass spectra[END_REF]) or on the input covariance matrix [START_REF] Ottesen | Structural correlation method for model reduction and practical estimation of patient specific parameters illustrated on heart rate regulation[END_REF]. For other feature selection techniques and an overview of the field, the interested reader is referred to [START_REF] Guyon | An introduction to variable and feature selection[END_REF]. Other methods consist in selecting directions in the parameter space. In other words, instead of selecting a subset of features, linear combinations of the features are sought. In Principle Component Analysis directions of maximum variance in the parameter space are sought [START_REF] Arquès | Study of a perturbation in the coding periodicity[END_REF]. The same principle holds for functional-PCA [START_REF] Gokulakrishnan | A functionalpca approach for analyzing and reducing complex chemical mechanisms[END_REF], its counterpart applied to the case where the input space is a function space. Neither of these approaches take into account the relationship between inputs and outputs. In Active Subspaces [START_REF] Constantine | Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies[END_REF], directions in the input space are sought so that the gradient of the output with respect to these directions is maximum. In Partial Least Squares [START_REF] Wold | The collinearity problem in linear regression. the partial least squares (pls) approach to generalized inverses[END_REF], directions are sought in the parameter space so that their covariance with the output is maximum. In this regard, this approach bears some similarities with the present work. Indeed, we look for biomarkers that are maximally correlated with their respective parameters. The main difference is that we also add the constraint that they are minimally correlated with all remaining parameters. Another way to reduce the input space is to perform a sparse linear regression using the Lasso algorithm [START_REF] Young | Kernel methods and haplotypes used in selection of sparse dna markers for protein yield in dairy cattle[END_REF]. The use of the 1 norm penalization makes this approach similar to the present work even though the cost function to be minimized is different.

Another aspect which makes the present approach different from previous works is that the feature, or biomarker, selection is performed in order to simplify future inverse problems. This issue has been addressed in [START_REF] Davidson | The dimensional reduction method for identification of parameters that trade-off due to similar model roles[END_REF] but is rarely the focus of feature selection studies. For a comprehensive review of inverse problem techniques, the reader is referred to [START_REF] Kaipio | Statistical and computational inverse problems[END_REF].

Our method is based on a semi-empirical approach. A mathematical model of the system of interest is considered and a database of simulations is built, by taking meaningful scenarios into account. Then, a dictionary of linear and nonlinear forms of the observable is considered. Such a dictionary is user-defined and therefore the efficiency of the proposed method relies on a relevant choice of its entries. Typically, such a dictionary is built by considering first common or "classical" biomarkers defined in the literature and then by incorporating agnostic features to enrich the space spanned by the dictionary entries. The composite biomarker is defined as a linear combination of the elements of the dictionary. The linear combination is sought such that the resulting biomarker is maximally correlated to the hidden quantity it refers to, and minimally correlated with respect to all the others.

From a practical point of view, at the computational cost of one single offline database computation (done once for all), the expansion coefficients of the biomarkers on a dictionary of observable forms are computed, for a given experimental setup or physical system. Then this result can be exploited for an unlimited number of experiments. As a by-product, when doing parameter identification, the 2 -distance in the space of the biomarkers defines a metric which is the 2 -distance in the space of the hidden quantities, up to a controlled perturbation (see Proposition 1). This makes the inverse problem less ill-conditioned and, in general, easier to solve.

The structure of the work is as follows. In section 2, the notation is introduced and the method is described. In section 3, a property of the proposed composite biomarkers is analyzed. The numerical methods used to solve the problem in practice are described in Section 4. In section 5, three numerical experiments are shown. In the first one, a synthetic case is considered to illustrate the approach and highlight its features and performances. Then, two realistic examples are considered, one in the context of electrophysiology and the other in hemodynamics. In both cases, our algorithm provides composite biomarkers which lead to more efficient inverse problems.

Problem setting

Notation

Let (Θ, A, P) be a complete probability space, Θ being the set of outcomes, A a σ-algebra and P a probability measure. Let E [•] denote the expectation operator. In the following, random variables are denoted by boldface letters.

We assume we have an ODE or PDE model, depending on p parameters θ 1 , . . . , θ p , which, for the sake of simplicity, are assumed to be zero-mean unit-variance and mutually uncorrelated (see remark 1).

The observable of the model is defined as a function v from R p to R M which corresponds to the output of the parameterized equations. The observable can be viewed as a random vector v of the following form:

v : (Θ, A, P) -→ R M (θ 1 , . . . , θ p ) -→ v (θ 1 , . . . , θ p ) . (1)
Let us now define a biomarker g as a function from R M to R. A biomarker can also be viewed as a random variable g defined as follows:

g : (Θ, A, P) -→ R (θ 1 , . . . , θ p ) -→ g (v (θ 1 , . . . , θ p )) . (2)
We introduce the biomarkers dictionary, a set of N g biomarkers G = g 1 , . . . , g Ng . The biomarker g k is referred to as the k-th entry of the dictionary. The purpose of this work is to derive p composite biomarkers corresponding to the p parameters. The composite biomarker y (j) associated with the parameter θ j is sought as a linear combination of the dictionary entries:

y (j) = Ng k=1 w (j) k g k , ∀j ∈ {1, . . . , p}, (3)
where w

(j)
k , 1 ≤ k ≤ N g , 1 ≤ j ≤ p, are the weights of the dictionary decomposition and the unknowns of the problem. Analogously, a random composite biomarker is defined as follows:

y (j) = Ng k=1 w (j) k gk , ∀j ∈ {1, . . . , p}, (4)
where gk = g k -E [g k ], g k is the k-th random dictionary entry as defined in [START_REF] Cintrón-Arias | A sensitivity matrix based methodology for inverse problem formulation[END_REF]. The vector of the unknowns (w

(j) 1 , . . . , w (j)
Ng ) is denoted by w (j) in what follows. Example: To illustrate the notation, consider the following parametric ODE model:

u(t; θ 1 , θ 2 ) = θ 1 u(t; θ 1 , θ 2 ), t ∈ [0, T ], u(0; θ 1 , θ 2 ) = θ 2 .
(

) 5
The output of the model is assumed to be made of the values of the solution u at M time instants t 1 , . . . , t M . Thus the observable v = (v 1 , . . . , v M ) is defined as:

v : R 2 -→ R M (θ 1 , θ 2 ) -→ (u(t 1 ; θ 1 , θ 2 )), . . . , u(t M ; θ 1 , θ 2 ))) . ( 6
)
The dictionary is assumed to be made of N g = 3 biomarkers G = {g 1 , g 2 , g 3 }, g i : R M -→ R, defined for example by:

g 1 : (v 1 , . . . , v M ) -→ max (v 1 , . . . , v M ) , g 2 : (v 1 , . . . , v M ) -→ min (v 1 , . . . , v M ) , g 3 : (v 1 , . . . , v M ) -→ 1 M M i=1 v i .
We will not discuss the choice of the entries of the dictionary: we just assume that they are given and can be applied to the measurements of the phenomenon of interest.

Our wish is to build a quantity y (1) (resp. y (2) ) that is well suited to the identification of the parameter θ 1 (resp. θ 2 ). More precisely, the algorithm proposed in this work will find two vectors w (1) and w (2) of R 3 which define the two composite biomarkers as:

y (1) = w (1)
1 g 1 + w

(1) 2 g 2 + w

(1)

3 g 3 , y (2) = w (2) 1 g 1 + w (2) 2 g 2 + w (2) 3 g 3 .
Once the problem is solved, Proposition 1 below can give an a posteriori hint about the quality of the dictionary with respect to the parameter estimation problem.

The following quantities will be useful to the definition of the algorithm:

C k,l := E [g k θ l ] , for 1 ≤ k ≤ N g , 1 ≤ l ≤ p, V k,l := E [g k gl ] , for 1 ≤ k ≤ N g , 1 ≤ l ≤ N g , ν(w (j) ) = w (j)T V w (j) , for 1 ≤ j ≤ p, e (j) = (e (j)
l ) l=1,...,p , with e (j)

l := δ lj , for 1 ≤ j ≤ p, 1 ≤ l ≤ p, (7)
where δ lj is 1 if l = j and 0 if l = j.

Problem formulation

For j ∈ {1, . . . , p}, we define the composite biomarker y (j) by Ng k=1 w (j) k g k , where (w (j) k ) k=1,...,Ng is solution to the minimization problem:

w (j) * = arg min w (j) ∈R Ng J λ (w (j) ) (8)
with

J λ (w (j) ) = J (w (j) ) + λ (j) N g w (j) 1 (R Ng ) , (9)
and J (w

(j) ) = 1 2 C T w (j) -e (j) 2 2 (R p ) + ξ 2 ν w (j) -1 2 , ( 10
)
where the parameters λ (j) and ξ are assumed to be given for now. Note that the first part of J (w (j) ) is convex while the second one is not.

It follows from the definition of y (j) and the definitions in (7) that

C T w (j) = E θy (j) , ν(w (j) ) = E y (j) 2 . ( 11
)
Let us explain the rationale of problem [START_REF] Constantine | Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies[END_REF]. We remind the reader that the parameters θ 1 , . . . , θ p are supposed to be zero-mean unit-variance random variables. Thus, the minimization of C T w (j)e (j) 2

2 (R p ) implies that y (j) , defined by

Ng k=1 w (j)
k gk , has a covariance with parameter θ j (i.e. E y (j) θ j ) close to 1, and a covariance with the other parameters (i.e. E y (j) θ l , ∀l = j) close to 0. The penalty term ξ 2 ν w (j) -1 2 will force the variance of y (j) to be close to one. The penalty term λ (j) Ng w (j) 1 (R Ng ) tends to promote sparsity. The motivation behind the sparsity promotion is twofold. First, the obtained biomarkers, if most values of w (j) are zero, can be more easily interpreted. Second, 1 -penalized models tend to be more predictive and less prone to overfitting than 2 -penalized ones, as proved for example in the case of logistic regression in the comparative study [START_REF] Ng | Feature selection, l 1 vs. l 2 regularization, and rotational invariance[END_REF].

Remark 1. The parameters θ i have been assumed to be zero-mean unit-variance and mutually uncorrelated. In the general case, the mean and the covariance of the parameters read:

θ i := E[θ i ], S ϑ ij := E[(θ j -θ j )(θ i -θ i )].
When the matrix S ϑ is positive definite, there exists an affine change of coordinates such that the transformed parameters have the properties assumed in the derivation of the method.

Remark 2. Problem (8) corresponds to a penalization of a minimization problem which has a unique solution (minimization of the strictly convex function J λ over the compact set defined by ν w (j) = 1).

Analysis of a property of the composite biomarkers

In this Section, we propose an estimation which illustrates the benefits of the composite biomarkers as defined above. For the sake of simplicity, the analysis is carried out in absence of 1 -penalization, i.e. it is assumed that λ (j) = 0.

From a parameter identification perspective, an ideal configuration would be to have the j-th composite biomarker y (j) proportional to the parameter θ j . So we propose here to assess how far we are from this configuration and to relate this distance to the value of the functional J (j) . The result can also be viewed as an indication about the quality of the biomarkers dictionary G.

The composite biomarker y (j) , after the coefficients w (j) have been found, can be seen as a function of the random variables θ 1 , . . . , θ p . We remind that the θ j are zero-mean unit-variance and mutually uncorrelated random variables, i.e E[θ j ] = 0, E[θ j θ l ] = δ jl , ∀j, l. The expression of the biomarker y (j) can be decomposed into a linear part with respect to the parameter vector θ = (θ 1 , . . . , θ p ) and a non-linear part as follows:

y (j) = θ T α (j) + q (j) , (12)
where α (j) ∈ R p and q (j) : (Θ, A, P) → R are such that E[q (j) θ l ] = 0, ∀l.

Proposition 1. Decomposing the biomarker y (j) as in (12), we have:

q (j) 2 ≤ K J (j) 1/4 ,
where K is a constant, J (j) = J (w (j) ) is defined as the minimum of (10), and

v 2 2 = E[v 2 ].
Proof. For convenience, the vector α (j) is decomposed in e (j) + γ (j) , where e j) . Inserting the formulas in [START_REF] Davidson | The dimensional reduction method for identification of parameters that trade-off due to similar model roles[END_REF] into Eq.( 10) leads to:

(j) l = δ lj . Let ∆ (j) = θ T γ (
J (j) = 1 2 E θy (j) -e (j) 2 2 (R p ) + ξ 2 E y (j) 2 -1 2 , ( 13
) = 1 2 E [θθ j ] + E θθ T γ (j) -e (j) 2 2 (R p ) + ξ 2 E θ 2 j + ∆ (j)2 + q (j)2 + 2∆ (j) θ j -1 2 , (14)
= 1 2 γ (j) 2 2 (R p ) + ξ 2 E ∆ (j)2 + E q (j)2 + 2E ∆ (j) θ j 2 , (15)
= 1 2 ∆ (j) 2 2 + ξ 2 q (j)2 2 2 + ∆ (j) 2 2 + 2E ∆ (j) θ j 2 . ( 16
)
First, we notice that ∆ (j) 2 ≤ √ 2 J (j) 1/2 . Then, applying the Cauchy-Schwarz inequality and observing that θ j 2 = 1, ∀j, we have:

q (j) 2 2 ≤ 1 √ ξ 2J (j) -∆ (j) 2 2 1/2 -∆ (j) 2 2 -2E[∆ (j) θ j ] ≤ 2 ξ J (j) 1/2 -∆ (j) 2 2 + 2 ∆ (j) 2 ≤ 2 ξ J (j) 1/2 + 2 ∆ (j) 2 , ≤ K 2 J (j) 1/2 (17) with K 2 = 2/ξ + 2 √ 2. ♦
The result of the proposition shows that q (j) , which is the nonlinear part (with respect to the parameters) of the biomarker y (j) , has a norm • 2 controlled by the value of the functional J (j) , which is minimized. If the value of the functional is close to zero, then, the norm of q (j) is close to zero. Thus, the biomarker, which is a nonlinear map of the observable, is close to a linear map of the parameters. In addition, this linear map is a perturbation of the identity. The amplitude of this perturbation is measured by the quantity ∆ (j) 2 introduced in the proof of Proposition 1, and it was noticed that

∆ (j) 2 ≤ √ 2 J (j) 1/2
. This shows that if the value of J (j) is small, then the 2 -distance of the biomarkers is close to the 2distance of the parameters. This property is particularly appealing when considering parameter identification problems, as will be shown below.

In conclusion, it is desirable that the value of J (j) be as small as possible. Thus, the value reached by J (j) once the minimization has been performed can be seen as an a posteriori indicator of the quality of the dictionary. Remark 3. If the smallness of J (j) is not deemed satisfactory, the biomarkers dictionary may be augmented with additional entries. Knowing a priori the typical number of entries in the dictionary required to obtain a small enough cost function is a challenging task. This important question goes beyond the scope of the present paper.

Numerical methods

Discretization

The expectation operators appearing in the definition [START_REF] Gokulakrishnan | A functionalpca approach for analyzing and reducing complex chemical mechanisms[END_REF] are approximated by a Monte-Carlo method. In the offline phase, N samples of the random parameters (θ 1 , . . . , θ p ) are drawn and stored in the matrix P ∈ R N ×p . For each sample (θ 1,i , . . . , θ p,i ), i ∈ {i, . . . , N }, the computational model is evaluated and v (θ 1,i , . . . , θ p,i ) is computed. The entries of the matrix G ∈ R N ×Ng are then defined by:

G i,k = g k (v (θ 1,i , . . . , θ p,i )) - 1 N N i=1 g k (v (θ 1,i , . . . , θ p,i )). ( 18
)

Solution of the optimization problem

The 1 penalty term makes the objective function J λ not differentiable. A standard way to deal with this difficulty is to introduce the change of variables w (j) = w

(j) + -w (j)
-, where w

(j) + (resp. w (j)
-) is the positive (resp. negative) part of w (j) [START_REF] Bach | Convex optimization with sparsity-inducing norms[END_REF]. The minimization problem can then be recast as follows:

       min w (j) + ,w (j) - J (w (j) + -w (j) -) + λ (j) Ng 1 T w (j) + + 1 T w (j) - s.t. w (j) + ≥ 0 w (j) - ≥ 0 , ( 19
)
where 1 denotes the vector of ones in R Ng . Note that while the non-differentiability issue of the 1 -norm is avoided, the number of unknowns is doubled and 2N g inequality constraints are added.

Many algorithms can be used to solve the constrained optimization problem in [START_REF] Davies | An in silico canine cardiac midmyocardial action potential duration model as a tool for early drug safety assessment[END_REF]. In the present work, we choose the Nesterov accelerated gradient descent method [START_REF] Nesterov | A method for unconstrained convex minimization problem with the rate of convergence o (1/k2)[END_REF][START_REF] Donoghue | Adaptive restart for accelerated gradient schemes[END_REF] in which, at each step, the unknowns are projected onto the constraint set if they do not satisfy the inequality constraints. Given the type of constraints, this projection is straightforward. The gradient to be used in the Nesterov gradient descent iterations reads:

∂ w (j) + J λ = 1 N G T P ( 1 N P T Gw (j) -e (j) ) + 2ξ N ( 1 N w (j)T G T Gw (j) -1)G T Gw (j) + λ (j) N g 1, ∂ w (j) - J λ = - 1 N G T P ( 1 N P T Gw (j) -e (j) ) - 2ξ N ( 1 N w (j)T G T Gw (j) -1)G T Gw (j) + λ (j) N g 1.
(20)

Choice of the penalty parameter ξ

A very simple automatic tuning of the penalty parameter ξ is performed along the resolution of the minimization problem [START_REF] Davies | An in silico canine cardiac midmyocardial action potential duration model as a tool for early drug safety assessment[END_REF]: starting from an initial guess, it is increased until the constraint is verified up to a certain tolerance parameter c . The strategy is summarized in Algorithm 1 in the Annex A1.

Choice of the penalty parameter λ (j)

Two strategies are proposed to set the regularization parameters λ (1) , . . . , λ (p) : "L-curve" criterion and training error threshold.

"L-curve" criterion. A common way to set the regularization parameters involves computing the so-called "L-curve" which represents the balance between the sparsity of the solution and the corresponding training error. Here, the sparsity of the solution is measured by w (j)

1 (R Ng ) . The training error corresponding to the biomarker y (j) is defined as 1 N P T Gw (j)e (j) 2 (R p ) and quantifies how well the discrete biomarker y (j) = Gw (j) fulfills the objectives for a given training set {P, G}. The problem in [START_REF] Davies | An in silico canine cardiac midmyocardial action potential duration model as a tool for early drug safety assessment[END_REF] is solved for different values of λ (j) and therefore one L-curve is obtained for each biomarker y (j) (i.e. for each parameter θ j ) for j = 1, . . . , p. In Figure 1 an example of L-curves is provided when considering the first two parameters of the model presented in Section 5.3. A natural way to determine λ (j) from the 0.0 0.1 0.2 0.3 0.4 0.5 Here, λ (j) varies between 10 -6 and 20.

1 N P T Gw (j) -e (j) ℓ 2 (R p ) 0.0 0.1 0.2 0.3 0.4 0.5 1 N g w (j) ℓ 1 (R Ng ) λ =1e-06 λ =0.001 λ =0.02 λ =2.0 λ =20.0 j = 4 j = 5
L-curve is to choose the λ (j) that maximizes its curvature. It is, in practice, a good compromise between an over-fitted (λ (j) is too small) and an over-regularized (λ (j) is too large) solution, as detailed in [START_REF] Johnston | Selecting the corner in the l-curve approach to tikhonov regularization[END_REF] where several approaches are compared. It also has the advantage of not requiring other hyperparameters.

Training error threshold. Another way to set the regularization parameters is to choose the largest λ (j) so that a user-defined training error is respected. While easy to implement, this requires an additional tuning parameter.

Numerical Experiments

Parameter identification algorithm

One of the outcomes of the proposed methodology is to make the parameter estimation easier. Before illustrating this feature in numerical experiments, we briefly present in this section the strategy used to identify parameters.

The parameter estimation is formulated as an optimization problem, in which the misfit between measurements and the simulated observables is minimized. The cost function has the following form:

J 1 (θ) = 1 2 v(θ) -v * 2 2 (R M ) , (21)
where v * ∈ R M denotes the vector of measurements from which the parameters true values θ * are to be estimated. In the following numerical experiments, v * is simply a noisy evaluation of the model:

v * = v(θ * ) + (η 1 , . . . , η M ) ( 22
)
where the η i are independent realizations of η ∼ N (0, τ 2 ).

When considering the composite biomarkers, the cost function reads:

J 2 (θ) = 1 2 y(v(θ)) -y(v * ) 2 2 (R p ) . ( 23
)
The minimization is performed with a simplification of Powell's BOBYQA algorithm [START_REF] Powell | The bobyqa algorithm for bound constrained optimization without derivatives[END_REF] which relies on successive quadratic approximations of the cost function.

Starting from an initial guess, the cost function is locally evaluated and approximated by a quadratic function from which the minimum is easily computed. This iterative process does not require any gradient evaluation and converges, in practice, faster than the classical gradient descent. The algorithm, later referred to as gradient-free minimization, is presented in Algorithm 2 in the Annex A2.

A synthetic model

To start with, we illustrate our algorithm with a simple algebraic model defined as follows:

u(x; θ 1 , θ 2 ) = x 2 exp - (θ 2 1 + θ 2 2 -R 2 ) 2 2σ 2 0 + θ 1 exp - (x -(x 1 + αθ 1 )) 2 2σ 2 1 + θ 2 exp - (x -(x 2 + αθ 2 )) 2 2σ 2 2 , x ∈ [-1, 1] (24)
with the following fixed parameters: R = 0.75, x 1 = -0.5, x 2 = 0.5, α = 0.1, σ 0 = 0.2 and σ 1 = σ 2 = 0.1. Here, x is the state variable and the uncertain parameters are θ 1 and θ 2 . The segment [-1, 1] is discretized into a uniform grid of M = 201 nodes x 1 , . . . , x M . The observable v is defined as follows:

v : R 2 -→ R M (θ 1 , θ 2 ) -→ (u(x 1 ; θ 1 , θ 2 )), . . . , u(x M ; θ 1 , θ 2 ))) (25)
The biomarkers dictionary is composed of N g = M entries defined as follows:

g k : R M -→ R (v(θ 1 , θ 2 ) 1 , . . . , v(θ 1 , θ 2 ) M ) -→ v(θ 1 , θ 2 ) k , 1 ≤ k ≤ N g . (26)
Note that, in this particular example, (g 1 (v), . . . g Ng (v)) and v coincide but this is not always the case. More complex transformations of the observable will be presented in the next numerical examples.

The biomarkers dictionary entries are the values of the model at each node of the discretized domain u(x 1 ), . . . , u(x Ng ). The training set is generated with N = 10 3 samples of (θ 1 , θ 2 ) drawn from the normal distribution N (µ, Σ), where µ = (0, 0) and Σ = 2.5 • 10 -1 × I 2 .

Biomarkers computation

The procedure described in Section 4 is applied to the present test case. The thresholding method is used to choose the regularization parameters λ (1) and λ (2) . The method in Algorithm 1 is used with a constraint tolerance c = 10 -3 . The threshold on the training error is set to 10 -2 which yields λ (1) = 0.11 and λ (2) = 0.09. The obtained biomarkers weights w (1) and w (2) are represented in Fig. 3 as well as the composite biomarkers y (1) and y (2) with respect to the parameters. Note that most coefficients are zero, as expected with the sparsity promoting 1 penalty.

Parameter estimation

Here, we highlight how using the previously obtained composite biomarkers leads to easier inverse problems. In the context of a parameter estimation problem, one seeks the true values of the uncertain parameters θ * from a measurement v * as defined in [START_REF] Kaur | Fitting membrane resistance along with action potential shape in cardiac myocytes improves convergence: application of a multi-objective parallel genetic algorithm[END_REF] with θ * = (0.3, 0.3) and τ = 2 • 10 -2 . We propose to perform a parameter estimation using J 1 and J 2 with a gradient-free minimization. The cost functionals J 1 and J 2 as function of θ 1 , θ 2 are plotted in Fig. 4. While J 1 features local minima and is non convex, the functional J 2 has only one global minimum, corresponding to the true solution. The results of the parameter estimations are summarized in Table 1. As expected, the parameter estimation fails when using J 1 in the sense that the distance to the true solution at the end of the minimization procedure is large. The minimization method falls, for J 1 , into a local minimum as seen in Fig. 4. On the contrary, the minimization succeeds in finding the true solution when using J 2 (using the composite biomarkers).

Influence of hyperparameters

Using the analytical model, we now investigate the influence of the regularization parameters λ (1) , λ (2) and the training set parameters samples (denoted in matrix form by P ). For the sake of simplicity and because θ 1 and θ 2 have essentially symmetrical effects, the study will be restricted to the cases where λ (1) = λ (2) . In Figure 5, the cost function J 2 (θ 1 , θ 2 ) associated with the inverse problem in ( 23) is represented for a regularization parameter increasing from 10 -5 to 10 1 . For high values of λ (j) , the cost function loses its convexity since the penalty term in (8) becomes predominant. In Figure 6 is represented J 2 (θ 1 , θ 2 ) for different training set parameters samples P . Using the same number of samples (N = 10 3 ) and same standard deviation, the samples are drawn from a normal distribution with different means. The method is robust with respect to the choice of training set in the sense that the cost function remains convex even when the training set center is far from the true solution, except in the last case (µ = (-2, -2)).

A model in electrophysiology

In this section, the proposed methodology is applied to a cardiomyocyte (heart muscle cell) electrophysiology model. Cardiomyocytes are contractible cells that react to an electrical activation. Detailed electrophysiology models have been developed to reproduce the cardiomyocytes complex electrical activity. Such models consist in general in a set of nonlinear coupled ODEs of the following form:

           dV dt = m i=1 Φ i (V, t), Φ i (V, t) = κ i f i (V, z i ) for 1 ≤ i ≤ m, dz i dt = ζ i (V, z i ) for 1 ≤ i ≤ m. ( 27
)
θ (n iter ) -θ * 2 (R p ) n iter J 1 1.31 × 10 0 7 J 2 6.14 × 10 -3 18
Table 1: Parameter estimation using Powell's method. Both approaches (J 1 and J 2 ) are compared in terms of number of iterations to reach convergence and distance from the converged solution θ (niter) to the true solution θ *

14

-1

0 1 θ 1 -1 0 1 θ 2 λ = 10 -5 -1 0 1 θ 1 -1 0 1 θ 2 λ = 10 -3 -1 0 1 θ 1 -1 0 1 θ 2 λ = 10 -1 -1 0 1 θ 1 -1 0
J 2 (θ 1 , θ 2 )
Figure 5: Cost function J 2 (θ 1 , θ 2 ) with different regularization parameters λ = λ (1) = λ (2) and parameters samples centered on µ = (0, 0). The red dot materializes the position of the true solution θ * = (0.3, 0.3).

-1 Despite a wide variety of models, most of them share this common structure. The integer m corresponds to the number of different kinds of ionic currents and ranges from 1 to more than 30 depending on the model's precision level. In this paper, the model by Davies et al. [START_REF] Davies | An in silico canine cardiac midmyocardial action potential duration model as a tool for early drug safety assessment[END_REF] was chosen. It counts m = 17 ionic currents, 29 ODEs and is designed to model cardiomyocytes belonging to the ventricular region of the canine heart. The observable quantity is the action potential (AP), denoted by V in [START_REF] Reymond | Validation of a one-dimensional model of the systemic arterial tree[END_REF]. It corresponds to the electrical potential difference between the inside and the outside of the cell. It drives the heart contraction and its alteration is at the origin of many cardiac diseases. In Figure 8 is plotted a typical AP along with some biomarkers that are commonly extracted from it. In our case, the parameters of interest are the conductances, κ i in [START_REF] Reymond | Validation of a one-dimensional model of the systemic arterial tree[END_REF]. They are the parameters commonly calibrated to fit different kind of cells or to model diseases and drug-induced AP alterations. We chose to retain 5 of these conductances as the parameters of interest, namely D Kr , D K1 , C to1 , D Ca,L and D N aCa . Each parameter θ j , 1 ≤ j ≤ is a multiplicative factor of the corresponding parameter reference value, which is equivalent to making the following modifications in the original model:

0 1 θ 1 -1 0 1 θ 2 µ = (0, 0) -1 0 1 θ 1 -1 0 1 θ 2 µ = (-0.5, -0.5) -1 0 1 θ 1 -1 0 1 θ 2 µ = (-1, -1) -1 0 1 θ 1 -1 0 1 θ 2 µ = (-2, -2) 0.0 0.6 1.2 1.8 2.4 3.0 3.6 J 2 (θ 1 , θ 2 )
D Kr = θ 1 DKr , D K1 = θ 2 DK1 , C to1 = θ 3 Cto1 , D Ca,L = θ 4 DCa,L , D N aCa = θ 5 DNaCa ,
where the bar symbol refers to the published values. Estimating this kind of parameters has many applications in electrophysiology, ranging from fitting models to experiments [START_REF] Syed | Atrial cell action potential parameter fitting using genetic algorithms[END_REF][START_REF] Sarkar | Regression analysis for constraining free parameters in electrophysiological models of cardiac cells[END_REF][START_REF] Kaur | Fitting membrane resistance along with action potential shape in cardiac myocytes improves convergence: application of a multi-objective parallel genetic algorithm[END_REF] to obtaining insights into disease-induced cell modifications [START_REF] Sánchez | Inter-subject variability in human atrial action potential in sinus rhythm versus chronic atrial fibrillation[END_REF]. Most often, this inverse problem is solved by using biomarkers.

Figure 9 shows how correlated these 5 parameters of interest are with their classically associated biomarkers. The latter are, respectively, APD 90 , RMP, V notch , APD 50 /APD 90 , APD 50 (see Figure 8). This choice of biomarkers reflects the common associations made in the literature [START_REF] Davies | An in silico canine cardiac midmyocardial action potential duration model as a tool for early drug safety assessment[END_REF][START_REF] Courtemanche | Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model[END_REF][START_REF] Koivumäki | In silico screening of the key cellular remodeling targets in chronic atrial fibrillation[END_REF]. Without surprise, this result shows that the classically adopted biomarkers are indeed correlated with the parameters of interest. Nevertheless, it also shows that each feature is not maximally correlated with its corresponding parameter and that it is also correlated with the other parameters. There is therefore room for improvement, and we propose to test if our strategy to define composite biomarkers can lead to a better parameter estimation.

We apply our method to a dictionary of 104 entries computed from V (t) + η, where η ∼ N (0, τ ) with τ = 0.5. More details about the dictionary entries may be found in the Annex A3. This corresponds to a SNR of 43 dB. The dictionary entries are the AP values at different times, the products of AP values at different times, the time derivatives, the integral over time, the PCA coefficients and the "classical" biomarkers (see Figure 8). The training set is generated with N = 5 × 10 4 samples drawn from a uniform distribution over the hypercube [0.5, 3] 5 . This means the conductances κ i are allowed to vary from half to three times their reference values. The penalization parameter ξ is fixed and set to 1. A posteriori verifications show that the error on the variance constraint does not exceed 10 -2 in practice.

The regularization parameters are set using the L-curve criterion.

As an example, Figure 7 shows the resulting biomarker for parameter θ 3 . Interestingly, the algorithm finds that the composite biomarker is a correction of V notch , which is indeed the biomarkers classically associated with parameter θ 3 . The obtained composite biomarkers are now used to solve parameter estimation problems. Figure 10 shows an example of such an inverse problem. Given a certain true value θ * and initial guess θ (0) , both in [0.5, 3] 5 , the parameter estimation is carried out by minimizing J 1 (θ) and J 2 (θ) with the gradient-free method described in Algorithm 2. In this specific example, minimizing J 1 fails to reach the true solution because it falls into a local minimum, whereas minimizing J 2 leads to the true solution.

V notch

We now propose to illustrate the advantage of using composite biomarkers instead of the time series to identify parameters. Since the efficiency and the accuracy of the parameter estimation strongly depends on the initial guess and the sought values, we carry out a statistical survey: 25 samples are randomly drawn from [0.5, 3] 5 and an inverse problem is executed for each two-point combination. This amounts to 625 inverse problems. Each inverse problem is carried out twice by minimizing J 1 and J 2 and the estimation error is stored. Here, the estimation error is simply defined as the 2-norm between the true solution and the solution of the minimization problem. Figure 11 shows the histogram of estimation errors in both cases. When using the composite biomarkers instead of the whole time series, the estimation error is in average smaller and is less spread (lower standard deviation). In the first case, the estimation error has mean 0.12 and standard deviation 0.16. In the second case, the estimation error has mean 0.44 and standard deviation 0.43.

A model of systemic blood flow circulation

A last an application is provided concerning the hemodynamics in a networks of arteries. The motivation is for example the monitoring of hypertensive patients, for which it is important to assess arteries stiffening from non-invasive measurements. The biomarker which is typically used in practice is the Pulse Wave Velocity (PWV), which is the average speed of the pressure waves travelling from carotid to femoral arteries. But arterial stiffness is not the only property affecting the PWV. Daily fluctuations of PWV can be explained as the interplay of different phenomena occurring in the cardiovascular system, such as the different metabolic need of the peripheral organs or the fluctuation in the heart rate due to physical activity. In what follows, the possibility of correcting PWV is studied, in order to have a biomarker which is more intrinsically related to arteries stiffness.

The main 55 arteries of the human body are considered, and the fluid-solid interaction occurring in blood vessels is described by means of a simplified 1-D model (see e.g. [START_REF] Formaggia | Cardiovascular Mathematics: Modeling and simulation of the circulatory system[END_REF][START_REF] Reymond | Validation of a one-dimensional model of the systemic arterial tree[END_REF][START_REF] Matthys | Pulse wave propagation in a model human arterial network: assessment of 1-d numerical simulations against in vitro measurements[END_REF]). The resulting system is a set of coupled hyperbolic partial differential equations in the unknowns (A i , q i ), where A i (x, t) is the crosssectional area of the i-th vessel and q i (x, t) is the sectional averaged blood velocity, Figure 10: Example of a parameter estimation from a synthetic AP using the gradient-free minimization algorithm. v * is the synthetic measurement, θ (0) the initial guess and θ (J1) (resp. θ (J2) ) the solution of the minimization of J 1 (θ) (resp. J 2 (θ)). x ∈ [0, i ] being the axial coordinate of the i-th vessel, t ∈ R + being the time. The system takes the following form:

∂ t A i + ∂ x (A i q i ) = 0, (28)
∂ t q i + ∂ x q 2 i 2 + p i = -κ q i A i , (29)
where is the blood density, κ a friction coefficient, p i (x, t) is the pressure in the vessels which is related to the cross-sectional area through the following algebraic relation:

p i = p 0 + β A 1/2 i -A (0) 1/2 i
, where β is the elastic coefficient [START_REF] Formaggia | Cardiovascular Mathematics: Modeling and simulation of the circulatory system[END_REF] and A (0) i is the cross sectional area of the i-th vessel, at rest. The boundary conditions for the model are of three natures. First, the flow Q(t) is imposed at the inlet of the ascending aorta as follows:

A 1 (0, t)q 1 (0, t) = Q(t). (30)
Q(t) is user-defined and is supposed to mimic the heart activity. Second, the effects of the organs and micro-circulation are modeled by standard three-element Windkessel models at the outlets [START_REF] Matthys | Pulse wave propagation in a model human arterial network: assessment of 1-d numerical simulations against in vitro measurements[END_REF]. Third, boundary conditions at the bifurcations are set by enforcing mass conservation and continuity of total pressure. The model is discretized and solved following [START_REF] Audebert | Kinetic scheme for arterial and venous blood flow, and application to partial hepatectomy modeling[END_REF] where a kinetic scheme is considered. The model, having more than 100 free parameters, suffers from a severe identifiability issue when using realistic data. Hence, the model has been reparametrised by introducing 8 hyper-parameters, accounting for the main source of variability in the observed data. In particular, the main goal is to be able to represent (in a simple way) the daily fluctuations of the blood flow in the network as well as the variability in the material properties of the arteries. The stiffness of the different segments has been parametrized as follows:

β i = θ 1 + θ 2 β (ref ) i , (31)
where

β (ref ) i
is a reference elastic coefficient for the i-th vessel, available in the literature [START_REF] Reymond | Validation of a one-dimensional model of the systemic arterial tree[END_REF]. The three-element Windkessel model parameters have been grouped into three sets: the first one gathers the terminal models of the upper part of the body, the second one gathers those of the central part and the last one the lumped models of the lower body. For each, the distal resistance is rescaled, and the capacitance is set in such a way that the characteristic time of the circuit remains unchanged:

R u,j = θ 3 R (ref ) u,j , C u,j = C (ref ) u,j θ 3 , (32)
R c,j = θ 4 R (ref ) c,j , C c,j = C (ref ) c,j θ 4 , (33)
R l,j = θ 5 R (ref ) l,j , C l,j = C (ref ) l,j θ 5 . (34)
Three additional parameters account for variations in the cardiac output and cardiac rhythm:

Q(t) = θ 6 sin π θ 7 Ts t 0 ≤ t ≤ θ 7 T s 0 θ 7 < t ≤ θ 8 T , (35)
where T s = 0.25s is the systolic period and T = 0.8s is a reference cardiac period.

The training set is generated by drawing N = 1024 samples of θ = (θ 1 , . . . , θ 8 ) from the normal distribution N (µ, Σ), where µ = (0, 1, 1, 1, 1, 1, 1, 1) and Σ = 0.05 2 × I 8 . The flow and pressure signals are computed by integrating the model for these samples. With the variability considered, the fluctuation in terms of pressure and flow are about 20%, which is enough to cover meaningful scenarios of daily fluctuations.

The dictionary is built as follows. In three distinct locations -the carotid, the subclavian artery and the femoral artery -average in time, maximum and minimum value of pressure, flow, and mechanical power P m = pAq are taken. The last three elements of the dictionary are PWV taken from carotid to femoral, from carotid to subclavian and from subclavian to femoral. The dictionary is ordered as follows: G = {(Aq) 0 , . . . , (Aq) 8 , p 0 , . . . , p 8 , P m,0 , . . . , P m,8 , P W V 0 , P W

V 1 , P W V 2 } .
The method is used with λ (j) = 10 -3 , ∀j ∈ {1, . . . , 8}. The penalization parameter ξ is fixed and set to 1. A posteriori verifications show that the error on the variance constraint does not exceed 10 -2 in practice.

Figure 12: Correlation of PWV with respect to all parameters (black) and correlation of y (2) , the composite biomarker associated with θ 2 , with respect to all parameters (red).

In Fig. 12, the red squares (resp. the black dots) show the correlation coefficients between the biomarker y 1 provided by our algorithm (resp. the PWV) and each of the 9 parameters. The PWV is significantly correlated to θ 1 , which is not a surprise since the arterial stiffness is known to be linked to the PWV. But the PWV also has a significant correlation (of about 0.3 -0.4) with the peripheral resistances of the central body and with the cardiac output. This provides a possible explanation for the observed daily fluctuations of the PWV. On the contrary, the biomarker y 1 provided by our method has a higher correlation with respect to θ 1 and a negligible correlation with respect to all the other parameters.

Interestingly, the decomposition of the biomarker onto the dictionary sheds some light onto the physical interpretation of the correction. In Fig. 13 the expression of w (2) is shown. We can observe that it is sparse. As expected, the PWV from carotid to subclavian is selected. Then, the correction is provided by a combination of the mechanical power. It is remarkable that the relevant quantity proves to be the product of the flow and the pressure, but not the flow and the pressure separately.

PWV

Figure 13: Biomarker expression on the dictionary, for y (2) In conclusion, if we are interested in estimating the arterial stiffness, the algorithm suggests that it could be more relevant to consider a combination of PWV and mechanical power rather than just the PWV. This result is preliminary and would require a deeper investigation. Nevertheless, it shows an example of possible applications of the proposed method to correct the expression of standard biomarkers used in hemodynamics.

Conclusions and Perspectives

In this work a method is proposed to define composite biomarkers. A semiempirical approach is used, that consists in building a database of simulations and exploiting the results in order to define biomarkers which are maximally correlated to the parameters they are meant to monitor and minimally correlated with respect to all the others. Such an objective is translated into an optimization problem and the biomarkers are sought as the solution of its discretized version. The biomarkers are, by construction, linear combinations of dictionary entries that comprise relevant features found in the literature and additional quantities computed from the observable signal. Furthermore, a sparsity-promoting penalization ensures that the obtained biomarker representation on the dictionary entries are sparse and, often, easily interpretable. As a by-product, a metric in the space of the observable is obtained which is equivalent, up to a small perturbation, to the 2 -distance in the space of parameters. The final cost function of the optimization problem can be used as an a posteriori estimate of the biomarkers dictionary quality. In case this cost is deemed to high, additional entries may be added to the dictionary using agnostic or model-driven quantities. The method can be exploited in order to enhance the robustness of the biomarkers which are used in biology and medicine by using numerical simulations. Furthermore, it can be used to regularize inverse problems and make them less ill-conditioned.

The proposed approach is particularly suited to biomedical applications where biomarkers are often preferred to raw measurements. However, the authors believe the methodology presented in this work could find applications beyond this scope. For all the systems in which several signals are measured and predictions, identification, statistical learning tasks have to be performed, the proposed approach could be of some help. For instance, applications could be considered in environmental sciences, ecology, climate monitoring or geophysics.

Acknowledgments

Eliott Tixier was funded by a doctoral grant from the French Ministry of Research and Higher Education.

A2. Algorithm to solve the parameter identification problem

Input:

• Cost function: J(θ) (see Eqs. ( 21),( 23))

• Initial guess: θ (0)

• Bracket size: γ

• Stopping criteria: absTol and relTol i ← 0; r (0) := J(θ (0) ); while r (i) < absTol or r (i)r (i-1) < relTol × r (i) do for j = 1 to p do Fit quadratic model Q(s) to f (s) = J(. . . , θ

j-1 , s, θ end Update residual: r (i+1) := J(θ (i+1) ); i ← i + 1 ; end Algorithm 2: Alternate directions descent quadratic approximations for the parameter identification problem (simplification of Powell's BOBYQA algorithm [START_REF] Powell | The bobyqa algorithm for bound constrained optimization without derivatives[END_REF]).

A3. Action potential biomarkers

The 104 dictionary entries for the action potential model used in Section 5.3 are hereby described.

Action potential durations. The first entries correspond to action potential durations (APDx) at x% repolarization. They are defined as the amount of time during which the action potential remains above its peak value minus x% of its amplitude. There are 8 APD values in the dictionary: APD90, APD80, ..., APD20.

Other classical biomarkers. Some classical biomarkers are also added, respectively APD50/APD90, APA (AP amplitude), RMP (resting membrane potential).

Action potential values.

The following entries are action potential values at 61 different time steps between 0ms and 600 ms. V notch is the 6 th of these entries.

Derivatives and integrals. The following entries are dV/dt max (maximum derivative during the upstroke), dV/dt min (minimum derivative during the repolarization) and AUC (area under curve). Auto-correlation. 21 entries correspond to products of the action potential values at different time steps of the form: V (t i ) * V (t j ).

PCA coefficients. Finally, 8 entries correspond to the first 8 coefficients of a truncated PCA expansion of V (t).

Figure 1 :

 1 Figure1: The "L-curve" is obtained for the first two parameters of the test case presented in Section 5.3. The 1 -norm of w (j) is plotted against the training error for different values of λ(j) . Here, λ (j) varies between 10 -6 and 20.

5 Figure 2 :

 52 Figure 2: Model outputs for different values of θ 1 , θ 2 .

Figure 3 :

 3 Figure 3: Analytical model: biomarkers weights (top) and corresponding biomarkers (bottom).

Figure 4 :

 4 Figure 4: Contours of the inverse problem cost function with true solution θ * = (0.3, 0.3) represented by a red dot. J 1 using the raw model outputs (left) and J 2 using the composite biomarkers (right).

Figure 6 :

 6 Figure 6: Cost function J 2 (θ 1 , θ 2 ) with different training sets and λ (1) = λ (2) = 10 -3 . The red dot materializes the position of the true solution θ * = (0.3, 0.3).

Figure 7 :

 7 Figure 7: Biomarker weights w (3) associated with θ 3 . The composite biomarker is actually a correction of the classically associated one V notch .

Figure 8 :

 8 Figure 8: Action potential (V (t) in (27)) computed using the Davies model. Common features associated with the parameters of interest used in the literature are also represented.

Figure 9 :

 9 Figure 9: Covariance matrix between the five parameters of interest of the Davies model and their classically associated biomarkers.

Figure 11 :

 11 Figure11: Histograms of the parameter estimation error for the 625 inverse problems carried out with both J 1 and J 2 as the inverse problem cost function.

  component of θ = (θ 1 , . . . , θ 8 )

Annexes

A1. Algorithm to solve the minimization problem [START_REF] Davies | An in silico canine cardiac midmyocardial action potential duration model as a tool for early drug safety assessment[END_REF] Input:

• Sparsity promoting parameter: λ (j)

• Initial guess: w (j),0

• Stopping criteria: absTol and relTol

// increase constraint penalization parameter end Compute gradient: ∂ w (j) ± J λ (w (j),i );

// using (20) Update weights: w (j),i+1 := f (w (j),i , ∂ w (j) J λ (w (j),i ); // Nesterov formula Update residual: r (i+1) := J λ (w (j),i+1 ); i ← i + 1 ; end Algorithm 1: solve the minimization problem [START_REF] Davies | An in silico canine cardiac midmyocardial action potential duration model as a tool for early drug safety assessment[END_REF] with an automatic tuning of parameter ξ.