Marie Meyer-Monatha

Claudine Chatelliera

Deirdre Cabooterc

Florence Rougetd

Isabelle Morelb

Francois Lestremaua
email: francois.lestremau@ineris.fr

Development of liquid chromatography methods coupled to mass spectrometry for the analysis of substances with a wide variety of polarity in meconium

Keywords: reversed phase, ion exchange, HILIC, LC/MS/MS, meconium, micropollutants Highlights

Meconium is a fetal matrix which integrates a large exposure window to xenobiotics Analysis of 28 target compounds with a large range of polarity Comparison of six columns for analytical method development by LC/MS/MS in meconium Optimization of HILIC conditions for polar compounds in meconium analysis Abstract Meconium is the first fecal excretion of newborns. This complex accumulative matrix allows assessing the exposure of the fetus to xenobiotics during the last six months of pregnancy. To determine the eventual effect of fetal exposure to micropolluants in this matrix, robust and sensitive analytical methods must be developed. This article describes the method development of liquid chromatography methods coupled to triple quadrupole mass spectrometry for relevant pollutants. The 28 selected target compounds had different physico chemical properties from very polar (glyphosate) to non-polar molecules (pyrethroids). Tests were performed with six different columns: reversed phase, ion exchange and HILIC. As a unique method could not be determined for the simultaneous analysis of all compounds, three columns were selected and suitable chromatographic methods were optimized. Similar results were noticed for the separation of the target compounds dissolved in either meconium extract or solvent for reversed phase and ion exchange columns. However, for HILIC, the matrix had a significant influence on the peak shape and robustness of the method. Finally, the analytical methods were applied to "real" meconium samples.

Introduction

Like the general population, pregnant women are commonly exposed to many environmental pollutants (pesticides, solvents…) which have been determined harmful to humans due to their carcinogenic properties and/or their effects on reproduction. To determine the exposure to micropollutants of women during pregnancy, the PELAGIE project (Perturbateurs Endocriniens: Étude Longitudinale sur les Anomalies de la Grossesse, l'Infertilité et l'Enfance), wherein maternal urine and cord blood samples were analyzed, was carried out in Brittany, France [1]. The results of this study identified some important compounds of interest, such as organophosphorus pesticides and dialkylphosphates (metabolites of organophosphorus) [2]. The constant exposure of pregnant women to these compounds could affect the development of the fetus and especially induce the production of congenital malformations. To determine if an association exists between congenital malformations and fetal exposure to organic pollutants, a new project, the PENEW project (Pregnancy Environment and NEWborn malformations), has recently been launched by the Registry of congenital malformations in Brittany, France. This project includes the study of compounds quantified in the PELAGIE project (organophosphorus pesticides) and other relevant molecules widely used in Brittany and suspected to be toxic for the fetus such as other pesticides (triazoles, pyrethroids, glyphosate…) and VOCs; with some associated metabolites.

To evaluate direct fetal exposure to xenobiotics throughout pregnancy, fetal matrix was considered. Commonly used fetal matrices, such as cord blood, urine and newborn blood, only reflect the last days of exposure. Therefore, the meconium matrix was selected to represent a larger exposure window. Meconium is the first stool of a newborn. It starts forming during the 12th-13th week of gestation in intestinal compartments and accumulates until birth [3,4,5]. It is a very complex matrix [6] composed of water (~70%) and lipophilic compounds (~30%): lipids, proteins, bile acids, enzymes, lanugo… Usually, it is expelled by the newborn within 24 hours after birth. The collection of meconium is non-intrusive for the newborn and simple, unlike some other fetal matrices (urine, blood).

Meconium was already used to quantify fetal exposure to drugs and metabolites [7,8], alcohol metabolites [9,10] or pesticides [11,12] with the analysis performed by liquid chromatography coupled with tandem mass spectrometry LC/MS/MS. Usually these analytical methods were focused only on restricted classes of compounds and eventually their metabolites [8,9,13]. To determine a possible link with malformations, the target compounds of the PENEW project regrouped a large variety of pollutants within several families of pesticides, their metabolites, VOCs [14] and metabolites. Compared to a previous study carried out by our group where the analysis of certain micropollutants was performed on a C18 column [12], the current list of compounds included additional pesticides (epoxiconazole, tebuconazole, 2,4-D…) and also, notably, substances with stronger polarity (glyphosate) or different polar functional groups (acid compounds with the metabolites of VOCs). As the considered molecules encompassed a much larger variety in physicochemical properties, ranging from very polar to non-polar, the previously developed method was inadequate for the extended list of substances.

Analytical method development was therefore performed for the 28 substances in meconium.

Extraction methods for these compounds have been investigated and described elsewhere [15]. The aim of the method development by LC/MS/MS was to obtain one or several methods suited for all target substances with acceptable retention times (not eluted in the void volume), peak shapes and sensitivity with the fewest constraints. Analysis time was also considered due to the large number of samples to be analyzed (around 235 meconium samples). Tests on several columns with different properties, buffer optimization and matrix effects are described in this article. Chromatographic tests with meconium matrix were carried out after the sample preparation optimized in this matrix [15]. The final analytical method was applied to real PENEW "cases".

Experimental

Chemicals and materials

O-cresol (o-C, purity: 99.9%), hippuric acid (HA, 98.0%), 2-methylhippuric acid (2-MHA, 98.0%), 3-methylhippuric acid (3-MHA, 98.0%), 4-methylhippuric acid (4-MHA, 98.0%), phenylglyoxylic acid (PGA, 98.0%), S-phenylmercapturic acid (S-PMA, 99.0%), diethylthiophosphate (DETP, 98.0%), 2,4-dichlorophenoxyacetic acid (2,4-D, ≥ 98.0%), diazinon (purity: 98.3%), cypermethrin (98.5%), cyfluthrin (99.8%), deltamethrin (99.7%), clopyralid (99.3%), glyphosate (99.2%), tebuconazole (99.5%), propoxur (99.8%), dichlorvos (99.9%) and benzoic acid-d5 (BA-d5, 99.0%) were purchased from Sigma-Aldrich Chemie GmbH (Steinheim, Germany). S-benzylmercapturic acid (S-BMA, purity: 98.0%) was obtained from Sigma-Aldrich (Milwaukee, WI, USA). Chlorpyrifos (purity: 99.9%), malathion (97.2%), epoxiconazole (99.2%) were purchased from Riedel-de Haën (Seelze, Germany). Dimethylthiophosphate (DMTP, 99.2%), dimethyldithiophosphate (DMDTP, 99.4%), 3-(2,2-dichlorovinyl)-2,2-dimethylcyclopropanecarboxylic acid (DCCA, purity ≥ 95.0%, at 100 µg.mL-1 in acetonitrile), diethylthiophosphate-d10 (DETP-d10, purity: 98.0%, at 100 µg.mL-1 in methanol) were obtained from Cambridge Isotope Laboratories Inc.

(Andover, MA, USA). Dimethylphosphate (DMP, 98.0%) was obtained from Acros Organics (Geel, Belgium). Permethrin (99.8%) was purchased from Ultra Scientific (Rhode Island, USA). Mandelic acid (MA, 99.0%) was obtained from ChromaDex (Irvine, CA, USA). Diethylphosphate (DEP, 99.5%) was obtained from Chem Service Inc. (West Chester, PA, USA). 3-(2,2-dibromovinyl)-2,2-dimethylcyclopropane-1-carboxylic acid (Br2CA, 99.0%, at 10 µg.mL-1 in methanol), diazinon-d10 (96.5%, at 100 µg.mL-1 in acetonitrile), propoxur-d3 (99.5%, at 100 µg.mL-1 in acetonitrile), dichlorvos-d6 (96.0%, at 100 µg.mL-1 in cyclohexane), trans-cypermethrin-d6 (98.5%, at 100 µg.mL-1 in acetonitrile), chlorpyrifos-d10 (98.0%, at 100 µg.mL-1 in acetonitrile), tebuconazole-d6 (95.0%, at 100 µg.mL-1 in acetonitrile), 3-(2,2-dichlorovinyl)-2,2-dimethylcyclopropanecarboxylic acid-d6 (DCCA-d6, 96.0%, at 100 µg.mL-1 in acetonitrile), 2,4-dichlorophenoxyacetic acid-d3 (2,4-D-d3, 97.0%, at 100 µg.mL-1 in acetonitrile) were purchased from Dr. Ehrenstorfer (Augsburg, Germany).

Hippuric acid-d5 (HA-d5, 98.0%) was obtained from Toronto Research Chemicals Inc.

(North York, Canada).

Sodium acetate (NaOAc), ammonium hydroxide solution (NH4OH) (99.9%), and LC-MS grade methanol (MeOH) and acetonitrile (ACN) (≥ 99.9%) were obtained from Sigma-Aldrich (Steinheim, Germany). Ammonium acetate (> 98%) was obtained by Merck (Darmstadt, Germany). Ultrapure water was supplied by a Mili-Q water purifier system from Millipore (Bedford, MA, USA). Acetic acid (100%) and sodium sulfate anhydrous (NaSO4) were purchased from VWR BDH Prolabo (Leuven, Belgium).

Meconium

For method development, a pool of meconium was obtained from meconium samples collected from autopsied fetus and autopsied newborns by the anatomical pathology service of the University Hospital of Rennes, France. When tested, some of these meconium did not contain any of the analytes. These "blank" meconium were pooled to form a representative pool of "blank" meconium for use in all subsequent experiments. Meconium samples were kept at -80°C until analysis to prevent bacteriological developments.

Sample preparation

The extraction procedure of meconium specimens for quantification using the LC/MS/MS was the same as described previously [15].

Briefly, 1 g of meconium sample was diluted with 8 mL of water. Ammonium hydroxide and the internal standards were added to 6 mL of the diluted sample. After vortex and centrifugation, 4 mL of supernatant were filtered. To the remaining meconium phase, ACN, ammonium sulfate and ammonium acetate were added. The supernatant was evaporated at 500 µL. The two supernatants were pooled together before performing purification steps with a Bond Elut SAX cartridge (Agilent). The elutions of loading and washing steps of the Bond Elut SAX cartridge were purified with a Strata-X cartridge (Phenomenex). The elutions of the two SPE cartridges were pooled together before LC/MS/MS analysis.

Pyrethroids are not soluble in water. To avoid precipitation of pyrethroids or glyphosate in a unique solution, standard solutions were prepared in 100 µL of water/ACN (50/50 v/v) for analysis with Ascentis Express RP-Amide and Acclaim Trinity P1 columns (one injection of 10 µL for each column). For analysis with LUNA HILIC column, the 80 µL of remaining solvent (mixture of ACN and water) in vial were evaporated to dryness and reconstituted in 80 µL of ACN.

Liquid chromatography columns tests

The meconium and calibration samples were analyzed on an Acquity UPLC H-Class from Waters (Milford, MA, USA). Chromatographic development was accomplished with six columns: Ascentis Express C18, 150 × 2.1 mm, 2.7 µm from Supelco (Bellefonte, USA); Ascentis Express RP-Amide, 150 × 2.1 mm, 2.7 µm from Supelco (Bellefonte, USA); Ascentis Express Phenyl-Hexyl, 150 × 2.1 mm, 2.7 µm from Supelco (Bellefonte, USA); Acclaim Trinity P1, 150 × 2.1 mm, 3 µm from Thermo Fisher Scientific (Asheville, NC, USA); Ascentis Express HILIC, 150 × 2.1 mm, 2.7 µm from Supelco (Bellefonte, USA); and LUNA HILIC, 150 × 2.0 mm, 3 µm from Phenomenex (Torrance, CA, USA).

Optimized conditions for liquid chromatography

Chromatographic separation was carried out on three columns equipped with a pre-column: Acclaim Trinity P1, Ascentis Express RP-Amide and Luna HILIC. The mobile phase consisted of ammonium acetate buffer 20 mM set at pH 3 with acetic acid (solvent A) and acetonitrile (solvent B). The gradient program for the Ascentis Express RP-Amide column was as followed: isocratic hold at 100% A for 2 min; from 100% A to 100% B in 8 min; hold at 100% B for 2 min (total gradient time: 12 min). The gradient program for the Acclaim Trinity P1 column was: isocratic hold at 95% A for 2 min; from 95% A to 50% A in 3 min (total gradient time: 5 min). For the Luna HILIC column, a linear gradient from 100% B to 60% B in 5 min was used. The aqueous mobile phase was ammonium acetate buffered with acetic acid (mobile phase A) and acetonitrile as solvent (mobile phase B). The columns were constantly heated at 30°C for a better repeatability and the injection volume was set to 10 µL.

All the runs were performed at a flow rate of 400 µL min-1.

Using the current instrumentation, among different usage, the so called purge solvent is used to push the sample along the injection pathway. Therefore, the purge solvent comes into contact with the sample so its composition must be set close to the sample solvent composition. For Ascentis Express RP-Amide and Acclaim Trinity P1, the purge solvent consisted of water/MeOH (90/10: v/v). For the experiments performed with the LUNA HILIC column, the purge solvent consisted of ACN alone.

Mass spectrometry

The LC was coupled to a triple quadrupole mass spectrometer Acquity UPLC TQ Detector from Waters (Milford, MA, USA) in ESI positive and negative mode. The electrospray source was used at 150°C with a desolvation temperature of 400°C and a capillary voltage of ± 3.30 kV (Tables 1 and2). Nitrogen was used as nebulizer gas with desolvation gas flow set at 800 L Hr-1, cone gas flow set at 50 L Hr-1 and argon as collision gas with flow rate set at 0.17 mL min-1.

For MS detection, whenever possible, two ions were used for each compound in Multiple Reaction Monitoring (MRM) mode: one target ion for quantification and one qualifier ion for confirmation (Tables 1 and2). MassLynx software (Waters) was used for data acquisition and processing.

Results and discussion

The analysis of meconium is challenging because of the complex nature of this matrix. Each step of the analytical protocol has to be thoroughly optimized in order to reduce matrix effects. The determination of suitable chromatographic parameters is particularly critical to ensure a reliable and sensitive analysis. Separation from the matrix will ensure that interferences are reduced on peak shapes and sensitivity. Many components extracted from the matrix elute in the void volume so all target components must have a retention time of at least 2 times the void volume of the column. Moreover, even if detection by mass spectrometry is selective, two transitions are used to confirm the identity of the detected components by MS/MS, this significantly decreases the possibility of obtaining a false positive by wrongly identifying a component. For compounds with adequate retention time, the selectivity provided by mass spectrometry implies that compounds that co-elute can still be identified and quantified. However, the complexity of the matrix, despite intensive sample preparation, still has a large effect on the spectrometric response as demonstrated in [15].

Therefore, the developed method must represent the best compromise between the time of analysis and the efficiency of the separation.

For this particular study, the chromatographic analysis was also a significant part of this method development because of the large polarity range of the target molecules: from very polar (log Kow = -4.0 for glyphosate) to non-polar (log Kow = 5.7 for deltamethrin). The aim of this optimization was therefore to develop (an) easy and fast chromatographic method(s) for all target compounds.

The selection of columns was based on two criteria: adequate retention of the molecules on the column and good peak shape to maximize the sensitivity of the analysis.

Evaluation of the chromatographic separation with one column

To reduce the total analysis time, quantification of all target compounds using only one column would represent the ideal case. Considering the nature of the investigated compounds, from polar to non-polar, this column should display different types of interaction suitable for these molecules. Mixed-mode columns allow retaining different compounds with several interactions. Acclaim Trinity P1 is a trimodal column which allows different phase/molecule interactions: reversed-phase and ion exchange. These interaction modes depend on the ionization of the molecules. For chromatographic analysis, this state mainly depends on the pH of the buffer. Tests with different pH were therefore performed. These experiments were carried out with ACN and 20 mM buffer at pH 3 and 6. The flow rate was set at 400 µL min-1 and the gradient program was: isocratic hold at 95% A for 2 min; from 95% A to 95% B in 13 min; hold at 95% B for 2 min; then back to 95% A in 2 min and hold for 3.5 min at 95% A (total gradient time: 22.5 min). Most of the target compounds were separated both at acidic and neutral pH (fig. 1a and1b). For analysis at pH 3 (fig. 1a), only two molecules could not be eluted: PGA and clopyralid while DETP was detected at 18.5 min. At pH 6 (fig. 1b), glyphosate and DMP were not or very little retained, and 2,4-D eluted at 18.9 min. For the dialkylphosphates (DAPs), glyphosate, PGA, clopyralid and 2,4-D, the principal mechanism of retention was ion exchange. These compounds were more retained at pH 3, except for 2,4-D which was more retained at pH 6. Indeed, at pH 3, this molecule was partially ionized whereas it is completely ionized at high pH. As shown in fig. 1a and1b, several compounds such as the organophosphorus pesticides, propoxur, triazoles and pyrethroids were retained through a reversed-phase mechanism as their retention was not influenced by the buffer pH.

Overall, even if most of the target compounds eluted with a retention time ≥ 2 times the void volume, some were not retained and/or presented poor peak shapes at each pH. This method could therefore be relevant for studies where high concentrations are expected but for this particular work aiming at traces analysis, this constituted a limitation. Therefore the present analytical method was not selected for the quantification of all target compounds in meconium. Another method using multiple columns was rather envisaged to obtain good peak shapes and retention times for all studied compounds.

Evaluation of the separation with multiple columns

Reversed-phase optimization

As an analysis on only one column was not possible, tests with reversed phase columns were carried out to analyze most of the target compounds. Literature showed that this type of column was suitable for pesticide analysis in different matrices, such as water [16][17][18] and biological matrices like urine [19] or meconium [12]. In a previous work by our group [12], a C18 column was used to analyze pesticides and their metabolites in meconium. However, this column was not suitable for this work due to the larger range of polarity and functional groups for the considered compounds. Three reversed-phase columns with different functional groups were therefore tested: Ascentis Express C18, Ascentis Express RP-Amide, and Ascentis Express Phenyl-Hexyl. Ascentis Express C18 was chosen as a reference C18 reversed phase.

Ascentis Express RP-Amide was tested because it can improve the peak shape of some molecules compared to Ascentis Express C18. Finally, most of the target molecules had a phenyl group which could interact with a stationary phase containing a phenyl group such as Ascentis Express Phenyl-Hexyl.

For the tests of the 3 reversed phase columns, the same experimental conditions were used (gradient, mobile phase with buffer composed of ammonium acetate and acetic acid/ACN).

The flow rate was set at 400 µL min-1 and gradient program was isocratic hold at 95% A for 2 min; from 95% A to 95% B in 13 min; hold at 95% B for 2 min; then back to 95% A in 2 min and hold for 3.5 min at 95% A (total gradient time: 22.5 min).

The mobile phase pH is an important parameter for column selectivity particularly in this study given the high number of substances with ionizable groups. Therefore, buffers were tested at pH 3 and pH 6. Overall, compared to the results obtained with the mixed mode Acclaim Trinity P1 column, better peak shapes were obtained for most compounds using the reversed phase columns.

While most compounds, such as diazinon, remained unaffected by pH changes, some molecules with particular functional groups such as mandelic acid or S-PMA presented distorted peaks depending on the pH or the column (fig. S1).

At pH 6, acidic compounds such as mandelic acid or S-PMA are ionized and therefore elute at or close to the void volume. Moreover, peak shapes were very broad on all columns for S-PMA at pH 6. For this type of compound, the same behavior was noticed on the Acquity BEH C18 used at this pH in our previous work (data not shown) [12]. At pH 3, mandelic acid and S-PMA are mainly neutral. Therefore, interactions with these stationary phases are mainly governed by hydrophobic interactions.

At pH 3, broad peak shapes were observed for mandelic acid with Ascentis Express C18 and for Ascentis Express Phenyl-Hexyl columns but the Ascentis RP-Amide provided an acceptable result in term of peak shapes and retention times. At this pH, the Ascentis RP-Amide appeared more suitable for the analysis of polar compounds.

The peak shape and the retention of diazinon were not influenced by pH because this compound is in neutral form at pH 3 or 6.

Finally, with the 3 tested columns, glyphosate and DAPs (except DMDTP) were also not or poorly retained because of their high polarities. Therefore, the Ascentis Express RP-Amide column at pH 3 was selected (fig. 2) for the analysis of all compounds apart glyphosate and DAPs (except DMDTP).

HILIC optimization

Several published methods use C18 columns to analyze dialkylphosphates (DAPs) [20][21][22]. In our previous work [12], a C18 column was also used to analyze DAPs along with other pesticides. But as demonstrated, the buffer conditions which were used here (pH 6) were not suitable for the analysis of the acidic compounds in this study. Moreover, considering the sample preparation was changed to account for the wider range of polarities of the investigated compounds in this study, interferences during the spectrometric detection were noticed with the analysis of the very polar DAPs under these conditions. Therefore, Hydrophilic Interaction Chromatography (HILIC) was evaluated. As an alternative to normal phase, this interaction mode is used for the liquid chromatographic separation of polar and hydrophilic molecules that are not retained using reversed-phase liquid chromatography. The separation is based principally on electrostatic interactions (ionized analytes, ions contained in buffer (salt) and deprotonated silanols of stationary phase) and hydrophilic partitioning (between the buffer layer at the surface of the stationary phase and the organic layer) [23,24]. Many HILIC stationary phases are available: bare silica, polar functionalities bonded to silica or polymer coating bonded to silica. This type of column is very sensitive to chromatographic conditions and to matrix effects.

Optimization of the separation by HILIC was carried out by testing 2 different types of HILIC columns and by evaluating the nature and concentration of salt and the buffer pH.

The two HILIC columns tested were an Ascentis Express HILIC and a LUNA HILIC. The first column contains a fused core silica phase and the second has a cross-linked diol on the silica phase.

Tests with glyphosate demonstrated that this molecule was too strongly retained on the HILIC column and could hence not be eluted. Therefore, the HILIC method optimization was focused on the analysis of DAPs.

A mixture of DAPs prepared in ACN was analyzed with the Ascentis Express HILIC column with a gradient starting at 100% ACN. Indeed, when experiments with different gradients starting with low percentage of water (95% or 98%) were carried out, DETP was eluted in the void volume of the column. With the gradient starting at 100% ACN, all the injections performed with the pool of meconium were reproducible in term of retention times.

An experiment was performed to check the influence of the purge solvent. The purge solvent was a mixture of water/MeOH (90/10) as used for the analysis with the reversed-phase and Acclaim Trinity P1 columns. Under these conditions, two DAPs, DETP and DMTP, eluted in the void volume of the column. Therefore, the rest of the method optimization was carried out with a mixture of standards prepared in ACN and with a purge solvent consisting of ACN only.

Influence of buffer

Addition of salt (compatible with MS detection) to the mobile phase allows improving retention time, selectivity and peak shapes in HILIC [25,26]. Indeed, charges induced by the salt create electrostatic interactions with ionized compounds [27]. To determine the influence of this parameter on the two HILIC columns, two buffers were tested: 20 mM of ammonium acetate/acetic acid at pH 3 and 20 mM of ammonium formate/formic acid at pH 3. At the same pH, molecule retention time and peak shape depend on buffer type [26]. Better peak shapes and an increase in retention time were obtained for DEP and DMP for the Ascentis Express HILIC column with ammonium acetate as observed in other studies (fig. 3a and3b).

For the LUNA HILIC column, better peak shapes were also observed for DEP and DMP with this buffer (fig. 3b) while for DETP and DMTP, the two buffers provided poor retention (fig. 3b and3d). With these results, subsequent tests were carried out in ammonium acetate solution.

The effect of the salt concentration was investigated using 10, 20, 50, 80 and 100 mM of salt in aqueous solution. On the LUNA HILIC column (fig. S2a), a slight increase of retention was observed for DETP and DMTP when the salt concentration in the buffer was increased.

On the Ascentis Express HILIC column (fig. S2b), the retention times of DETP and DMTP decreased slightly with increasing ammonium acetate level. A decrease in DMP retention was observed from 10 to 50 mM of salt which then stabilized at higher concentrations. The same phenomenon (decrease of retention time) was noticed for DETP but in the range from 10 to 20 mM. An ammonium acetate buffer concentration of 20 mM was therefore selected.

Influence of pH

Four pH-values were tested: 3, 4, 5 and 6.5. The effect of mobile phase pH on the LUNA HILIC column is presented in fig. 4a. With 20 mM of ammonium acetate (full markers), retention times of DETP (pKa 2.7) and DMTP (pKa 2.9) slightly increase from pH 3 to pH 4 but decrease for higher pH-values to elute close to the void volume at pH 6.5. This trend can be correlated with the pKa of these 2 compounds which are mainly in neutral form up to pH 4 and are then ionized for higher pHs. For DMP and DEP, a decrease was observed from pH 3 to pH 6.5. Indeed, DMP (pKa 2.0) and DEP (pKa 1.9) are charged negatively at pH 3. From pH 4, the retention time of all compounds decreased because competition phenomena increased between ionized target molecules and increasing levels of acetate ions to create electrostatic interactions with the cross-linked diol phase. Overall, the same trend was observed with the formate buffer (empty symbols), retention times were stable from pH 3 to pH 5 and then decreased significantly.

On the Ascentis Express HILIC column (fig. 4b) and using ammonium acetate buffer (full symbols), no effect on the retention times of DMP and DEP was noticed when the pH varied between pH 3 and pH 5. At pH values higher than 5, the ionized DEP and DMP are in competition with acetate ions (more electronegative) for electrostatic interactions accounting for the decrease in retention time. For DETP/DMTP, retention times are stable over the entire pH range. As shown by the empty symbols, the retention times of all compounds were independent of the pH-value when using formate buffer. The same trend as observed with the acetate buffer was again seen for DETP and DMTP. However, for DMP and DEP, no decrease in retention was observed. The electronegativity of the formate ions might be less strong than those of DEP and DMP, so their electrostatics interactions are not modified with a pH change.

At pH 6.5, several compounds were poorly retained so this pH was not chosen. In general, variations of retention times between pH 3, 4 and 5 were minimal. Therefore, pH 3 was selected. With this choice, the buffer was the same as for the analysis with Ascentis Express RP-Amide.

The gradient was optimized to improve MS sensitivity, peak shapes and retention times: 100% of ACN to 60% at 5 min, then 100% of ACN at 7 min, hold 3 min; with 20 mM of ammonium acetate buffer at pH 3 with acetic acid used. In these conditions, all compounds were retained on the Ascentis Express HILIC column (fig. 3b) but peak shapes of DMTP and DETP remained unacceptable on the LUNA HILIC column (fig. 3a).

Choice of HILIC column for meconium analysis HILIC columns are very sensitive to matrix effects. To select a suitable column, a meconium extract containing DAPs was analyzed by LC/MS/MS. As shown in fig. 5, the chromatogram obtained for the meconium extract was very different from the chromatogram obtained in solvent (fig. 3). DMP was slightly affected by a matrix effect with good peak shape but the retention time shifted by 0.2 min. The retention time of DEP was not modified but diffusion effects were observed as for DETP and DMTP resulting in severe band broadening.

Interactions between these two last compounds and meconium extract were very important. Parts of the meconium matrix which were not removed during the purification process created interferences and interacted with the stationary phase and/or the target compounds. These interactions prevented the DAPs to link strongly with the bare silica phase.

Under identical analytical conditions, the same meconium extract was analyzed with the LUNA HILIC column (fig. 5a). In contrast to the Ascentis Express HILIC column (fig. 5b), the peak shapes were improved for DEP, DMTP and DETP compared to the analysis in solvent (fig. 3b and3d). It can be assumed that the meconium matrix had a positive effect by neutralizing the groups of the stationary phase responsible for the distorted peak shapes observed in solvent. This column was therefore selected for the analysis of the target molecules in the meconium samples.

To evaluate the importance of matrix effects further, a comparison with the chromatographic method proposed in this article and the one published by Odetokun et al. [28] in human urine was performed with the LUNA HILIC column. The chromatographic method of Odetokun et al. consisted of 93% of ACN and 7% of 100 mM ammonium acetate in isocratic mode with a flow rate at 400 µL min-1. In solvent, the best peak shapes were obtained for DETP and DMTP with the method of Odetokun et al. (fig. 6a) but for DEP and DMP, our optimized method provided better results. For the meconium extract (fig. 6b), our method again lead to the best peak shapes for all DAPs. This experiment demonstrates the particular nature and effect of the meconium matrix on the chromatographic HILIC conditions. Therefore, when dealing with complex matrix such as meconium, all final optimizations should preferentially be performed in the studied matrix.

Acclaim Trinity P1 optimization for glyphosate Usually, glyphosate is derivatized for chromatographical analysis on reversed-phase columns.

For underivatized glyphosate, columns with weak anion-exchange properties can be used to obtain sufficient retention and selectivity [29,30]. Indeed, in our tests, glyphosate was not retained on the reversed-phase columns and could not be eluted with the HILIC columns.

Further tests were therefore performed with the Acclaim Trinity P1 to benefit from the ion exchange interactions as seen in § 2.1. Retention with the Acclaim Trinity P1 depends very much on the buffer composition: salt concentration and pH. These two parameters were evaluated for glyphosate analysis.

Influence of pH

Several articles on underivatizated glyphosate analysis with LC using an anion-exchange column with acidic buffer have been published [29,31]. To verify these methods, two buffer pHs were tested with 20 mM of ammonium acetate: at pH 3 with acetic acid addition and pH 6 (pH limit conditions indicated by the manufacturer). The flow rate was set at 400 µL min-1 and ACN was used as solvent. The gradient program was isocratic hold at 95% A for 2 min; from 95% A to 95% B in 13 min; hold at 95% B for 2 min; then back to 95% A in 2 min and hold for 3.5 min at 95% A (total gradient time: 22.5 min). At pH 6, glyphosate eluted in the void volume of the column. Indeed, at this pH, anion exchange interactions are weak and hydrophobic interactions were insufficient to retain this compound as observed with reversed phase columns. Finally, a buffer set at pH 3 was selected since glyphosate was retained under these conditions.

Variation of salt concentration

The separation optimization on a trimodal column is influenced mainly by the effect of the mobile phase ionic strength (concentration and pH of buffer) [32]. The influence of the buffer salt concentration on the Acclaim Trinity P1 column was tested for different concentrations of ammonium acetate (5, 10, 15, 20 and 30 mM) at pH 3. The chromatographic conditions were the same as described in the previous paragraph. As shown in fig. S3, glyphosate was sensitive to the variation of ammonium acetate concentration. Indeed, the retention time of this molecule decreased with increasing buffer concentration. With 5 mM buffer, it was not eluted because the mobile phase had no sufficient ionic strength: interactions with glyphosate and the stationary phase were stronger than the glyphosate/mobile phase interactions. From 10 mM of ammonium acetate on retention increased slightly (6.42 min for 10 mM; 6.07 min for 15 mM and 4.62 min for 30 mM) and an improvement of the peak shape was observed with a reduction of peak tailing. A buffer containing 20 mM of ammonium acetate was therefore selected for the analysis of glyphosate on the Acclaim Trinity P1 column since it provided an acceptable result and for convenience of the overall analysis as this mobile phase was also used for the two other chromatographic methods in this work (reversed phase and HILIC).

Application to real samples

The full analytical protocol including sample preparation was validated and is described elsewhere [15]. It was applied overall to 235 meconium samples as part of the PENEW project. The meconium matrix is specific for each newborn; therefore the composition of the matrix varies for each analysis. Indeed, for the analysis of the 235 meconium samples with the LUNA HILIC column, an average retention time of 2.46 min was observed for the internal standard DETP-d10. The minimal retention time for this compound was 1.69 min and the maximum 3.35 min. This deviation was not present when the matrix matched calibration was carried out with the same pool of meconium and with the quality controls inserted within the sequence of samples analysis. This phenomenon was noticed only with the analysis performed with the HILIC column demonstrating this column is particularly sensitive to matrix effects.

Therefore, the use of a specific internal standard is required to take into account the variation of retention times produced by the difference in composition of the meconium samples.

Conclusion

Meconium is a complex accumulative matrix that can reflect fetal exposure to xenobiotics. To quantify 28 target compounds presenting a wide range of physico chemical properties in this matrix by LC/MS/MS, a test using only one column (of the most widely used type) did not provide satisfying results. Therefore an approach with three different complementary chromatographic methods was developed to obtain an orthogonal system suitable for all target compounds: one with a reversed-phase column for most of the target molecules (semi-polar to non-polar compounds), one with a HILIC column for polar metabolites (dialkylphosphates) and one with an ion-exchange column for a very polar pesticide (glyphosate). Analysis of meconium extracts was particularly critical using HILIC which demonstrated significant variation when analytical conditions/ matrix were changed. This analytical approach was successfully tested on 235 meconium samples. clopyralid;3,MA;4,5,HA;6,MHA;7,8,9,dichlorvos;10,propoxur;11,2,12,DCCA;13,epoxiconazole;14,malathion;15,DMDTP;16,Br2CA;17,tebuconazole;18,diazinon;19,cyfluthrine;20,chlorpyrifos;21,cypermethrin;22,deltamethrin;23,permethrin;24,glyphosate;25,DMP;26,DEP;27,DMTP;28,DETP Fig. 2 Analysis of meconium spiked with target compounds by Ascentis Express RP-Amide

Fig. 1

 1 Fig. 1 Acclaim Trinity P1 chromatogram in solvent : a) at pH3 and b) at pH6: 1, PGA; 2,

Fig. 3

 3 Fig. 3 Chromatogram of DAPs in solvent analyzed with 20 mM of ammonium acetate at pH 3

Fig. 3 Fig. 4

 34 Fig. 3 Chromatogram of DAPs in solvent analyzed with 20 mM of ammonium acetate at pH 3

Fig. 5 Fig. 6

 56 Fig.5Chromatogram of spiked meconium analyzed with 20 mM of buffer at pH 3, with a)

Table 1

 1 Chromatographic and MS parameters for target compounds

Acknowledgements

This work was funded par the Programme Hospitalier de Recherche Clinique national (PHRC) 2010 assigned to Rennes University Hospital (CHU) for the PENEW project. It was supported by Dr Florence Rouget (Brittany Registry of Congenital Malformations) in collaboration with the INSERM unities 1085 and 991 of Institut de Recherche en Santé, Environnement et Travail (IRSET) (Université Rennes I), and the Centre de Ressources Biologiques (CRB) of Centre d'Investigations Cliniques (CIC) of CHU de Rennes.

The authors thank Dr Loget of anatomical pathology service of Centre Hospitalier Universitaire (CHU) de Rennes for providing meconium samples.


N Zemzemi

R Aboulaich

N Fikal

E El Guarmah

Sensitivity of the electrocardiography inverse solution to the torso conductivity uncertainties

Sensitivity of the Electrocardiography Inverse Solution to the Torso Conductivity Uncertainties

Nejib Zemzemi, Rajae Aboulaich, Najib Fikal, Emahdi El Guarmah

Introduction

The ECGI procedure helps medical doctors to target some triggers of cardiac arrhythmia and consequently plan a much more accurate surgical interventions [START_REF] Shah | Non invasive ecg mapping to guide catheter ablation[END_REF]. The mathematical problem behind the ECGI solution is known to be ill posed [START_REF] Hadamard | Lectures on Cauchy's problem in linear partial differential equations[END_REF]. Therefore, many regularization techniques have been used in order to solve the problem [START_REF] Ghosh | Application of l1-norm regularization to epicardial potential solution of the inverse electrocardiography problem[END_REF][START_REF] Zakharov | Algorithms and numerical analysis of dc fields in a piecewisehomogeneous medium by the boundary integral equation method[END_REF][START_REF] Li | Localization of the site of origin of cardiac activation by means of a heartmodel-based electrocardiographic imaging approach[END_REF][START_REF] Doessel | Localization of the origin of premature beats using an integral method[END_REF][START_REF] Zemzemi | A machine learning regularization of the inverse problem in electrocardiography imaging[END_REF]. Thus, the obtained solution depends on the regularization method and parameters [START_REF] Hansen | The use of the l-curve in the regularization of discrete ill-posed problems[END_REF]. Furthermore, the inverse solution does not depend only on the regularization method but also to the physical parameters and the geometry of the patient. These dependencies have not been taken into account in most (if not all) of the studies. In particular, in most of the studies in the literature the torso is supposed to be a homogenous. Moreover, when the conductive heterogeneities are considered, they are fixed according to data obtained from textbooks. The problem is that these data are different from a paper to another [START_REF] Foster | Dielectric properties of tissues and biological materials: a critical review[END_REF][START_REF] Faes | The electric resistivity of human tissues (100 hz-10 mhz): a meta-analysis of review studies[END_REF][START_REF] Gabriel | The dielectric properties of biological tissues: Ii. measurements in the frequency range 10 hz to 20 ghz[END_REF]. Their variability is mainly due to the difference between the experiment environments and other factors related to the measurement tools. Only few works have been performed in order to study the effect of conductivities uncertainty in the propagation of the electrical potential in the torso [START_REF] Van Oosterom | The effect of torso inhomogeneities on body surface potentials quantified using tailored geometry[END_REF][START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ecg forward modeling to organ conductivity[END_REF][START_REF] Weber | Predicting tissue conductivity influences on body surface potentialsan efficient approach based on principal component analysis[END_REF]. In [START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ecg forward modeling to organ conductivity[END_REF], authors use the stochastic finite elements method (SFEM) to evaluate the effect of lungs muscles and fat conductivities on the forward problem. In [START_REF] Weber | Predicting tissue conductivity influences on body surface potentialsan efficient approach based on principal component analysis[END_REF], authors used a principal component approach to predict the effect of conductivities variation on the body surface potential.

However, to the best of our knowledge, no work in the literature has treated the influence of conductivity uncertainties of the ECGI inverse solution. In this work, we propose to use an optimal control approach to solve the inverse problem and to compute the potential value on the heart, we use an energy functional [START_REF] Aboulaϊch | Missing boundary data reconstruction via an approximate optimal control[END_REF]. with SPDE constraints as proposed in [START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ecg forward modeling to organ conductivity[END_REF]. We are interested in assessing the effect of the variability of tissue conductivity on the solution of the ECGI problem. The resolution of optimal control will be made using an iterative procedure based on the conjugate gradient method and the numerical approximation will be performed using the SFEM.

Methods

Solving stochastic forward problem of electrocardiography

Following [START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ecg forward modeling to organ conductivity[END_REF], we represent the stochastic characteristics of the forward solution of the Laplace equation by the generalized chaos polynomial. For the space domain we use simplified analytical 2D model representing a cross-section of the torso (see Figure 1) in which the conductivities vary stochastically. Let us first denote by D the spatial domain and Ω the probability space. By supposing that the conductivity parameter depends on the space and on the stochastic variable σ(x, ξ), the solution of the Laplace equation does also depend on space and the stochastic variable u(x, ξ)). The stochastic forward problem of electrocardiography can be written as follows:

   .(σ(x, ξ) u(x, ξ)) = 0 in D × Ω, u(x, ξ) = u 0 on Γ int × Ω, σ(x, ξ) ∂u(x,ξ) ∂n = 0 on Γ ext × Ω, (1)
where, Γ int and Γ ext are the epicardial and torso boundaries respectively, ξ ∈ Ω is the stochastic variable (it could also be a vector) and u 0 is the potential at the epicardial boundary.

Numerical descretization of the stochastic forward problem

We use the stochastic Galerkin (SG) method in order to solve equation [START_REF] Shah | Non invasive ecg mapping to guide catheter ablation[END_REF]. A stochastic process X(ξ) of a parameter or a variable X is represented by weighted sum of orthogonal polynomials known as probability density functions PDFs {Ψ i (ξ)}. More details about the different choices of PDFs could be found in [START_REF] Le Maître | A stochastic projection method for fluid flow: Ii. random process[END_REF]. In our case we use the Legendre polynomials which are more suitable for uniform probability density .

X(ξ) = p i=0 Xi Ψ i (ξ)
where Xi are the projections of the random process on the stochastic basis {Ψ i (ξ)} p i=1 . The mean value and the standard deviation of X over Ω are then computed as follows

E(X) = Ω p i=0 Xi Ψ i (ξ) = X0 , stdev[X] = p i=1 X2 i Ω Ψ i (ξ) 2 1
Since in our study we would like to evaluate the effect of the conductivity randomness of the different torso organs on the electrical potential. Both of σ and u are now expressed in the Galerkin space as follows:

σ(x, ξ) = p i=0 σi (x)Ψ i (ξ), u(x, ξ) = p j=0 ûj (x)Ψ i (ξ)
By substituting (2),(3) into the elliptic equation [START_REF] Shah | Non invasive ecg mapping to guide catheter ablation[END_REF] and by projecting the result on the polynomial basis {Ψ k (ξ)} p k=1 we obtain the following system: For k = 0, ..., p :

p i=0 p j=0 T ijk ∇.(σ i (x)∇)û j (x)) =0 in D, û0 (x) =u 0 (x) on Γ int , ûj (x) =0 on Γ int ∀j = 1, ...p, σi (x) ∂ ûj (x) ∂n =0 on Γ ext ∀ i, j = 0, ...p, (2)
where

T ijk = E[Ψ i (ξ), Ψ j (ξ), Ψ k (ξ)]
. By applying the standard finite elements variational formlation and Galerkin projections we obtain a linear system of size (p × dof ), where dof is the number of the degrees of freedom for the Laplace equation in the deterministic framework.

Fig. 1: 2D computational mesh of the torso geometry showing the different regions of the torso considered in this study (fat, muscle, lungs, torso cavity).

Stochastic inverse problem problem

In this section we propose to build an optimal control problem taking into account the variability of the tissue conductivities in the torso. Let us first propose a gold standard solution representing the true electrical potential in the torso. Supposing that the true conductivity distribution in the torso is σ T , at a given time the true electrical potential in the torso u T is governed by the Laplace equation, it depends on the extracellular potential at the epicardium u ex as follows:

   . σ T u T = 0 in D, u T = u ex on Γ int , σ T ∂u T ∂n = 0 on Γ ext . (3)
From the deterministic forward problem solution u T , we extract the electrical potential at the external boundary and we denote it by f = u T /Γext . For the inverse problem formulation, we write a cost function that takes into account the uncertainties in the torso conductivities. We then use an energy cost function as used in [START_REF] Andrieux | Solving cauchy problems by minimizing an energylike functional[END_REF][START_REF] Aboulaϊch | Missing boundary data reconstruction via an approximate optimal control[END_REF] constrained by the stochastic conductivity Laplace formulation as presented in the previous paragraph. We look for the current density and the value of the potential on the epicardial boundary

(η, τ ) ∈ L -1 2 (Γ int ) × L 1 2 (Γ int ) minimizing the following cost function                J(η, τ ) = 1 2 E v(x, ξ) -f 2 L 2 (Γext) + 1 2 σ(x, ξ) ∂u(x,ξ) ∂n -η 2 L 2 (Γint) with v(x, ξ) solution of : .(σ(x, ξ) v(x, ξ)) = 0 in D × Ω, v(x, ξ) = τ on Γ int × Ω, σ(x, ξ) ∂v(x,ξ) ∂n = 0 on Γ ext × Ω. (4)
In order to solve this minimization problem, we use a conjugate gradient method as used in [START_REF] Aboulaϊch | Missing boundary data reconstruction via an approximate optimal control[END_REF] where the components of the gradient of the cost function are computed using an adjoint method. The gradient of the functional J is given by:

                 < ∂J(η,τ ) ∂η .φ >= -E[ Γint (σ ∂v ∂n -η)φdΓ int ] ∀ φ ∈ L 2 (Γ int ), < ∂J(η,τ ) ∂τ .h >= E[ Γint σ ∂λ ∂n hdΓ int ] ∀ h ∈ L 2 (Γ int ), with λ solution of : ∇.(σ(x, ξ)∇λ(x, ξ)) = 0 on D × Ω, λ(x, ξ) = σ(x, ξ) ∂v(x,ξ) ∂n -η on Γ int × Ω, σ(x, ξ) ∂λ(x,ξ) ∂n = -(v -f ) on Γ ext × Ω.
(5)

The conjugate gradient algorithm.

The stochastic Cauchy problem has been reformulated as a minimization one in the previous section. For the numerical solution, we use a conjugate gradient optimization procedure. The algorithm is then as follows,

Step 1. Given f ∈ L 2 (Γ ext ), choose an arbitrary initial guess

(ϕ p , t p ) ∈ L 2 (Γ int ) × L 2 (Γ int )
Step 1.1. Solve the well-posed stochastic forward problem :

   ∇.(σ(x, ξ)∇v p (x, ξ)) = 0 in D × Ω, σ(x, ξ) ∂v p (x,ξ) ∂n = 0 on Γ ext × Ω, v p (x, ξ) = t p on Γ int × Ω. ( 6
)
Step 1.2. Solve the stochastic adjoint problem:

   ∇.(σ(x, ξ)∇λ p (x, ξ)) = 0 in D × Ω, λ p (x, ξ) = σ(x, ξ) ∂v(x,ξ) ∂n -ϕ p on Γ int × Ω, σ(x, ξ) ∂λ p (x,ξ) ∂n = -(v(x, ξ) -f ) on Γ ext × Ω. (7)
step 1.3. : We evaluate the gradient of the cost function:

∇J(ϕ p , t p ) = E ϕ p -σ(x, ξ) ∂v p (x, ξ) ∂n , E σ(x, ξ) ∂λ p (x, ξ) ∂n (8)
Step 1.4. Determine the descent direction d p as follows::

γ p-1 = ∇J(ϕp,tp) 2 ∇J(ϕp-1,tp-1) 2 d p := (d p 1 , d p 2 ) = -∇J(ϕ p , t p ) + γ p-1 d p-1 (9)
We compute: (ϕ p+1 , t p+1 ) = (ϕ p , t p ) + α p d p where the scalar α p is obtained through a linear search by :

α p = - E[ Γext z p (v p -f )dΓ ext ] + E[ Γint (σ ∂z p ∂n -d p 1 )(σ ∂v p ∂n -ϕ p )dΓ int ] E[ Γext (z p ) 2 dΓ ext ] + E[ Γint (σ ∂z p ∂n -d p 1 ) 2 dΓ int ] (10)
where z p is the solution of the following stochastic problem:

   ∇.(σ(x, ξ)∇z p (x, ξ)) = 0 in D × Ω z p (x, ξ) = d p 2 on Γ int × Ω σ(x, ξ) ∂z p (x,ξ) ∂n = 0 on Γ ext × Ω (11)
Step 2. Having obtained (ϕ p , t p ) for p ≥ 0, set p = p + 1 and repeat from step 1.1 until the prescribed stopping criterion is satisfied. As a stopping criterion we choose

E[v p (x, ξ)] -f L 2 (Γext) ≤ or ∇J(ϕ p , t p ) ≤ .
In practice, we take = 0.001.

Results

In this section we present the numerical results of the stochastic inverse problem. In order to assess the effect of the conductivity uncertainties of each of lungs and fat conductivities on the reconstructed electrical potential on the heart boundary, we start by generating our ground truth solution by solving the equation (3). For the sake of simplicity we take a harmonic function on the heart boundary: where the exact extracellular potential u ex (x, y) = exp(x) sin(y).

Since we assume that the uncertainty of the conductivity value follows a uniform probability density, as probability density functions {Ψ i } we use the Legendre polynomials defined on the interval Ω = [-1, 1]. We also suppose that the true conductivity uncertainty interval is centered by σ T , the true conductivity. We propose to study three cases. In the first case we suppose that there is no uncertainties. In the second (respectively, third) case we suppose that all the conductivities values are known but the conductivity of fat (respectively, lungs) where we gradually increase the uncertainty from zero to ±50% of the true conductivity value σ T . We solve the stochastic inverse problem following the algorithm described in the previous section. We measure the effect of the uncertainties using relative error (RE) and the correlation coefficient (CC). In table1, we show the performance of the stochastic method. We used different level of uncertainties: 0%, ±3%, ±10%, ±20%, ±30% and ±50%. We find that the relative error of the inverse solution has been barely affected by the uncertainties of the fat conductivity even for high uncertainty levels. In fact the RE (respectively, CC ) is 0.1245 (respectively, 0.993) when there is no uncertainties and for ±50% of uncertainties the RE (respectively, CC) is 0.1276 (respectively, 0.991). On the contrary the effect of the lung conductivity uncertainties is high: The RE increases from 0.1245 when we don't consider the uncertainties to 0.4852 when we have ±50% of uncertainties on the lung conductivity. The effect of the uncertainty on the correlation coefficient could also be qualitatively seen in Figure 2, where the pattern of the mean value of the stochastic inverse solution looks the same in Figure 2 compare these last figures to Figure 3 (right) we see that the deviation from the exect solution has been substantially increased. The same behavior is seen on the traces of the electrical potential on the int boundary: In Figure 4, we see that he inverse solution on the heart boundary does not change too much when considering the fat conductivity uncertainties and considerably changes when considering ±50% lungs conductivity uncertainties. compare these last figures to Figure 3 (right) we see that the deviation from the exect solution has been substantially increased. The same behavior is seen on the traces of the electrical potential on the int boundary: In Figure 4, we see that he inverse solution on the heart boundary does not change too much when considering the fat conductivity uncertainties and considerably changes when considering ±50% lungs conductivity uncertainties. compare these last figures to Figure 3 (right) we see that the deviation from the solution has been substantially increased. The same behavior is seen on the tra the electrical potential on the int boundary: In Figure 4, we see that he inverse so on the heart boundary does not change too much when considering the fat cond ity uncertainties and considerably changes when considering ±50% lungs conduc uncertainties. In Figure 2 we show the mean value of the stochastic inverse problem solution for the three cases. Qualitatively, the mean of the solutions didn't change too much. In Figure 3 we show the deviation of the inverse solution from the exact solution. Comparing Figure 3 (left) to Figure 3 (middle), we remark that the deviation from the exact solution did not change while considering the fat conductivity uncertainties. Whereas if we compare these last figures to Figure 3 (right) we see that the deviation from the exect solution has been substantially increased. The same behavior is seen on the traces of the electrical potential on the int boundary: In Figure 4, we see that he inverse solution on the heart boundary does not change too much when considering the fat conductivity uncertainties and considerably changes when considering ±50% lungs conductivity uncertainties. 3. As shown in Table 1, the error does not change too much from no uncertainties (Figure 3 (a)) to ±50% of fat conductivity uncertainty (Figure 3 (b)). Whereas the error is high for ±50% of lungs conductivity uncertainty (Figure 3 (c)). The propagation of uncertainties from the conductivities to the the inverse problem solution is reflected in the deviation of the stochastic inverse solution from the ground truth presented in Figure 4. We remark that the error is concentrated in the heart boundary Γ int , it reaches 0.8 for ±50% of lungs conductivity uncertainty and 0.25 for ±50% of fat conductivity uncertainty.

Discussion

Solving the inverse problem in electrocardiography imaging based on an optimal control problem allowed us to quantify the effect of the torso organs conductivity uncertainties using a stochastic finite element method. This work is the first mathematical study of the effect of the conductivity uncertainties on the inverse problem solution. Our results show that increasing the level of the fat conductivity uncertainty from zero to ±50% of its original value does not alter the quality of the reconstructed potential. This is in line with the results presented in [START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ecg forward modeling to organ conductivity[END_REF] for the forward problem when introducing ±50% uncertainties in the fat conductivity. On the contrary, the results that we obtained for the uncertainties on the lungs conductivity show an important effect on the ECGI solution. In fact the relative error is about 50% when introducing ±50% of uncertainty and the CC is significantly altered. This result is different from he results presented in [START_REF] Geneser | Application of stochastic finite element methods to study the sensitivity of ecg forward modeling to organ conductivity[END_REF] for the forward solution with ±50% uncertainties in the lungs conductivity where the standard deviation does not exceed 3% of the mean value.

conclusions

In this work we presented a novel approach to study sensitivity to parameters values in data completion inverse problem and that could have application in a wide range of bioelectric and biomedical inverse problems resolution. We used a stochastic finite element method in order to take into account the variability of the conductivity values in the ECGI inverse problem formulated in an optimal control manner. We used a conjugate gradient method to solve this problem where the gradient of the cost function was computed using an adjoint method. We have described the different steps of the algorithm used to solve this stochastic inverse problem. The numerical simulation that we conducted in a simple 2D case showed that there is an important sensitivity of the solution to the lungs conductivity, whereas the uncertainties on the fat conductivity did not affect too much the inverse solution. It is still not clear if the difference in the sensitivity to the uncertainties is due to volume of the organs and/or their positions and proximity to the heart and/or the magnitude of their mean conductivities. This question would be a subject of a mathematical and numerical investigation as a continuation of this work.

  (a) (no uncertainties) and (b) (±50% of uncertainties on the fat conductivity) and different in Figure 2 (c) (±50% of uncertainties on lungs conductivity). Similarly the effect of uncertainties on the relative error could be (a) (b) (c)

Fig. 2 :Fig. 3 :

 23 Fig. 2: Left (respectively, middle, right):Mean of the SFE solution with respect to null uncertainty in all organs model (respectively,±50% from the reference fat and lungs conductivity)

Fig. 2 :Fig. 3 :

 23 Fig. 2: Left (respectively, middle, right):Mean of the SFE solution with respect to null uncertainty in all organs model (respectively,±50% from the reference fat and lungs conductivity)

Fig. 2 :Fig. 3 :

 23 Fig. 2: Left (respectively, middle, right):Mean of the SFE solution with respect t uncertainty in all organs model (respectively,±50% from the reference fat and conductivity)

Fig. 2 :Fig. 3 :

 23 Fig. 2: Left (respectively, middle, right):Mean of the SFE solution with respect to null uncertainty in all organs model (respectively,±50% from the reference fat and lungs conductivity)

Fig. 2 :

 2 Fig. 2: Mean value of the SFEM inverse solution : No uncertainties(a), ±50% uncertainty in the fat conductivity(b) and ±50% uncertainty in lungs conductivity (c)

Fig. 3 :Fig. 4 :

 34 Fig. 3: Comparison of the SFEM solution to the ground truth on the heart boundary Γ int : (a) no uncertainties, (b) ±50% uncertainties in the fat conductivity, (c) ±50% uncertainties in lungs conductivity. Exact solution (blue continuous line). Stochastic inverse solution mean value (red dashed line). X-axis polar coordinate angle from -π to π. Y-axis value of the electrical potential on the boundary Γ int . (a) (b) (c)

Table 1 :

 1 Relative error and correlation coefficient of the stochastic inverse solution for different levels of uncertainty on the fat and lungs conductivities .

		conductivity uncertainties 0% ±3% ±10% ±20% ±30% ±50%
	relative error	fat	0.1245 0.1245 0.1248 0.1248 0.1251 0.1276
		lungs	0.1245 0.1263 0.1439 0.2208 0.3333 0.4852
	Corr coeff	fat	0.9930 0.9930 0.9945 0.9943 0.9980 0.991
		lungs	0.9930 0.9933 0.9899 0.9767 0.9660 0.885

Acknowledgments

We would like to thank the LIRIMA Laboratory which financially supported the teams ANO and EPICARD to perform this work.


Carlos Eduardo Chavez

Nejib Zemzemi

Yves Coudiere

Felipe Alonso-Atienza

Diego Álvarez

Inverse Problem of Electrocardiography: estimating the location of cardiac isquemia in a 3D geometry

The inverse problem in cardiology (IPC) has been formulated in different ways in order to non invasively obtain valuable informations about the heart condition. Most of the formulations solve the IPC under a quasistatic assumption neglecting the dynamic behavior of the electrical wave propagation in the heart. In this work we take into account this dynamic behavior by constraining the cost function with the monodomain model. We use an iterative algorithm combined with a level set formulation allowing us to localize an ischemic region in the heart. The method has been presented by Alvarez et al in [1] and [4], in which the authors developed a method for localize ischemic regions using a simple phenomenological model in a 2D cardiac tissue. In this work, we analyze the performance of this method in different 3D geometries. The inverse procedure exploits the spatiotemporal correlations contained in the observed data, which is formulated as a parametric adjust of a mathematical model that minimizes the misfit between the simulated and the observed data. We start by testing this method on two concentric spheres and then analyze the performance in a 3D real anatomical geometry. Both for analytical and real life geometries, numerical results show that using this algorithm we are capable of identifying the position and, in most of the cases, approximate the size of the ischemic regions.

Introduction

Cardiac arrhythmias are disturbances in the normal rhythm of the heart produced by an abnormal electrical activity. They are one of the major causes of death worldwide [START_REF] Shah | Body surface electrocardiographic mapping for non-invasive identification of arrhythmic sources[END_REF]. Cardiac arrhythmias can be analyzed by means of noninvasive procedures aiming to characterize the cardiac electric sources (membrane potential, epicardial or endocardical potentials, activation times) and/or the cardiac substrate (ischemia regions, post-infarct scars) from voltages recorded by electrodes systems which are not in direct contact with the cardiac tissue [START_REF] Rudy | Noninvasive electrocardiographic imaging of arrhythmogenic substrates in humans[END_REF][START_REF] Berger | Single-beat noninvasive imaging of cardiac electrophysiology of ventricular pre-excitation[END_REF]. From a mathematical point of view, the problem of recovering the heart electrical activity from remotes voltages can be formulated as an inverse problem, known as the inverse problem of electrocardiography (IPE) [START_REF] Gulrajani | The forward and inverse problems of electrocardiography[END_REF]. The IPE is a hard technological challenge since in its general formulation is an ill-posed problem so a number of regularization approaches have been developed over the years to obtain stable and realistic solutions [START_REF] Huiskamp | The depolarization sequence of the human heart surface computed from measured body surface potentials[END_REF][START_REF] Greensite | An improved method for estimating epicardial potentials from the body surface[END_REF][START_REF] Brooks | Inverse electrocardiography by simultaneous imposition of multiple constraints[END_REF].

In the present study we analyze the IPE in terms of localizing cardiac ischemic regions from remote voltage measurements. Cardiac ischemia is a pathology produced by the lack of blood supplied to the heart muscle that generates electrophysiologically-abnormal substrate which may lead to life threatening arrhythmias, and ultimately to a heart attack [START_REF] Lazzara | Ventricular arrhythmias and electrophysiological consequences of myocardial ischemia and infarction[END_REF]. In the clinical practice, determining the size and location of cardiac ischemic regions has shown limited accuracy, specially in severe damaged hearts [START_REF] Tokuda | Direct comparison of adjacent endocardial and epicardial electrograms: Implications for substrate mapping[END_REF]. Ischemia alters the propagation properties of the electrical impulse through the cardiac muscle and this results in alterations in the associated voltage measurement recordings. These alteration patterns are temporally stable, making it feasible to inversely reconstruct ischemic regions from voltage measurements recordings [START_REF] Wang | Inverse electrocardiographic source localization of ischemia: An optimization framework and finite element solution[END_REF].

Previous works have analyzed the localization of cardiac ischemia from remote voltage measurements [START_REF] Wang | Inverse electrocardiographic source localization of ischemia: An optimization framework and finite element solution[END_REF][START_REF] Li | Non-invasive estimation of myocardial infarction by means of a heart-model-based imaging approach: a simulation study[END_REF][START_REF] Maclachlan | Computing the size and location of myocardial ischemia using measurements of st-segment shift[END_REF][START_REF] Farina | Model-based approach to the localization of infarction[END_REF][START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF]. In [START_REF] Li | Non-invasive estimation of myocardial infarction by means of a heart-model-based imaging approach: a simulation study[END_REF][START_REF] Farina | Model-based approach to the localization of infarction[END_REF] the size and position of myocardial infarction is estimated by minimizing the difference between real voltage measurements and model-simulated ones. In [START_REF] Maclachlan | Computing the size and location of myocardial ischemia using measurements of st-segment shift[END_REF][START_REF] Ruud | A computationally efficient method for determining the size and location of myocardial ischemia[END_REF][START_REF] Nielsen | Computing ischemic regions in the heart with the bidomain model; first steps towards validation[END_REF] ischemic regions are assessed by reconstructing the epicardial potentials at a single time-instant during the plateau phase of the cardiac membrane potential by using a level-set formulation. Wang et al. [START_REF] Wang | Inverse electrocardiographic source localization of ischemia: An optimization framework and finite element solution[END_REF] formulated the IPE as a constraint minimization problem, in which the the size and position of ischemic regions were estimated by inversely computing the membrane potential at a single time instant during the plateau phase. While showing promising results towards clinical validation [START_REF] Wang | Inverse electrocardiographic source localization of ischemia: An optimization framework and finite element solution[END_REF][START_REF] Nielsen | Computing ischemic regions in the heart with the bidomain model; first steps towards validation[END_REF], in these studies time instants are treated independently from the others, thus ignoring the spatiotemporal correlation information contained in the voltage measurements. On the other hand, Álvarez et al. [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF] presented an algorithm being able to reconstruct disconnected cardiac ischemic regions with a limited number of recording sites by exploiting the spatiotemporal correlation contained in the measured data.

Here, we extend the work presented in [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF] to 3D geometries. We analyze the performance of the method in a model of spheres and then we study its use on 3D realistic anatomical model. Our results show that the algorithm is capable of identifying the position and, in most of the cases, approximates the size of the ischemic regions.

Forward Calculation

The forward problem refers to the mathematical models that describes the ionic processes involved in the generation of the action potential, its propagation through the cardiac tissue and remote recording model. In this section, we present the models used for the forward simulations.

Action potential model

Cellular electrical activity was simulated using a modified version of the Two-Current model (TC), proposed by Mitchell and Schaeffer [START_REF] Mitchell | A two-current model for the dynamics of cardiac membrane[END_REF]. The TC consists of just two ordinary differential equations for two variables: the transmembrane potential or voltage v(t) and the inactivation gate variable h(t). The voltage, which is dimensionless and scaled so that is varies between zero and one is defined as follows

dv dt = J T C = J stim (t) + J in (v, h) + J out (v) , (1)
where J stim represents the initial stimulus, J in and J out denotes the sum of all inward and outward currents, respectively. In order to simulate ischemia, we incorporated a parameter (v rest ) associated to the resting potential formulated by Alvarez et al [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF],

J in (v, h) = h(1 -v)(v -v rest ) 2 τ in , (2)
J out (v, h) = - v -v rest τ out . (3)
The gating variable h(t) is dimensionless and varies between zero and one. This variable regulates inward current flows and obeys the following equation

dh dt = G(v, h) = (1 -h)/τ open , v < v crit -h/τ close , v ≥ v crit (4)
This model contains four time constants (τ in , τ out , τ open and τ close ) which correspond to the four phases of the cardiac action potential: initiation, plateau, decay and recovery. The parameter v crit is the change-over voltage. We simulate the effects of ischemia modifying the values of τ in and v rest [START_REF] Chávez | Avoiding the inverse crime in the inverse problem of electrocardiography: estimating the shape and location of cardiac ischemia[END_REF]. For ischemic cells, we set the parameter τ in and v rest equal to 0.8 and 0.1, respectively.

Cardiac Tissue Model

The monodomain equation was use to simulate the electrical activity of the heart. Despite the discrete nature of cardiac cells structure, at the macroscopic scale cardiac tissue behaves as a functional syncytium. This permits to consider cardiac tissue (Ω H ) as an excitable medium in which the membrane potential propagation can be mathematically described according to the following reaction-diffussion equation

∂v ∂t = ∇ • (D∇v(r H , t)) + J T C ( 5
)
where D is the intracellular conductivity tensor, which is set to a constant value of D = 1.4mm 2 /ms. Imposing v(r H , 0) = v rest and h(r H , 0) = 1 as initial conditions and no-flux boundary conditions. J T C is taken as the ion current provided by TC model.

Remote Recording Model

The resulting potential distribution at position r T = (x T , y T , z T ) ∈ Ω T outside the cardiac tissue Ω H , is found by the addition (volume conductor theory) of the contribution of all source elements [START_REF] Malmivuo | Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields[END_REF][START_REF] Van Oosterom | Genesis of the p wave: Atrial signals as generated by the equivalent double layer source model[END_REF]. In this work, we consider remote recording measurements modeled as point electrodes at a distance R(r

H -r T ) = ||r T -r H || from source location r H = (x H , y H , z H ) ∈ Ω H to observation point r T i given by ϕ i (r T i , t) = 1 4πσ 0 Ω H ∇ • (D∇v(r H , t)) R(r H , r T i ) dΩ H (6)
where σ 0 represents the medium conductivity (assumed homogeneous), and v(r H , t) is solution of [START_REF] Farina | Model-based approach to the localization of infarction[END_REF]. Expressions ( 1)-( 6) comprise a complete description of the forward problem.

Inverse Procedure

The aim of the inverse procedure is to estimate the shape and locations of the ischemic areas adjusting the parameter distribution of τ in and v rest from the knowledge of the remote recordings measurements ϕ i (r T , t). The functional to be minimized is the L 2 norm of the misfit between the observed data and the simulated torso measurement, ϕ i R (r T i , t) and ϕ i S (r T i , t) respectively:

J (τ in (r H ), v rest (r H )) = 1 2 T 0 N i=1 ϕ i R (r T i , t) -ϕ i S (r T i , t; τ in , v rest ) 2 dt (7) = 1 2 T 0 N i=1 ϕ i R (r T i , t) - 1 4πσ 0 Ω H ∇ • (D∇v(r H , t)) R(r H , r T i ) dΩ H 2 dt
where T represents the recording time and N the number of electrodes.

In this approach, we follow an iterative scheme in which the cost functional is reduced at each step (ξ) of the reconstruction process. Therefore, we find a direction in the space parameter (τ in (r H ), v rest (r H )) such that the cost functional decreases. Since both parameters, τ in (r H ) and v rest (r H ), define the same region, we only consider the variation of τ in (r H ) for the gradient computation. The ischemic region is defined by a level set function as in [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF]. Subsequently, as we change the parameter distribution of τ in (r H ), distribution of v rest (r H ) is modified as well. Following [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF], the gradient of the function J over the parameter τ in is given by

grad τin J (τ in (r H )) = T 0 w(r H , t) h(r H , t) [v(r H , t) -v rest (r H )] 2 [v(r H , t) -1] τ 2 in (r H ) dt (8)
where adjoint state w is the solution of the following problem

       ∂w ∂t + ∇D • ∇w(r H , t) - ∂J T C ∂τ in w(r H , t) = - 1 4πσ 0 N i=1 R i Ω H ∇ 2 1 R(r H , r T i ) dΩ H dt w(r H , t = T ) = 0 (9) where R i = ||ϕ i R (r T i , t) -ϕ i S (r T i , t)||, and
∂J T C ∂τ in = h(r H , t) -v rest (r H ) τ in (r H ) (2 -3v(r H , t) + v rest (r H )) - 1 τ out (r H ) . ( 10
)
In this section, we present some numerical experiments in order of analyze the behavior of the inverse algorithm on 3D geometries. We conducted two scenarios:

(1) initialize the method with different initial guesses using a two concentric spheres geometry, (2) study the performance of the algorithm in a real anatomical geometry. For all experiments, the reconstruction algorithm is applied to a single cardiac cycle of length T = 240 ms in the steady state. The inverse procedure stops when the functional cost becomes stationary.

First Experiment

Simulation Setup. For the first scenario, we consider a spherical cardiac tissue of radius 50.0mm which might contain one or several ischemic regions. The cardiac tissue domain, Ω H , was discretized by a triangular finite element mesh generated with Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with builtin pre-and post-processing facilities[END_REF]. For the remote measurement points Ω T , we set a concentric sphere of radius 65.0mm.

Scenario. We suppose a circular ischemic region of radius 14 mm as is shown in Figure 1 panel (0). The aim of this experiment is to verify if the inverse procedure is able to reconstruct a single ischemia independently of the initial guess. We simulated three different cases: (a) supposing a healthy cardiac tissue i.e., there is not an ischemic region, (b) initialize the method with a ischemic region close to the real position and (c) considering an ischemic region elsewhere in the cardiac tissue (see Figure 1 panel (a1), (b1) and (c1), respectively). Figure 1 panel (a2), (b2) and (c2) shows reconstructed τ in for each case, respectively. A qualitative comparison between the obtained reconstruction and the real ischemic region shows that position and size of the ischemia was reconstructed successfully for all cases. We have computed the correlation coefficient between the true ischemia distribution and the distribution obtained from the inverse problem. For all cases, the correlation coefficient (CC) excess a value of 0.95 (CC1 = 0.9571, CC2 = 0.9543 and CC3 = 0.9532) which is consistent with the reconstructions observed in Figure 1.

Second Experiment

Anatomical model As a geometry, we use a 3D mesh of torso and heart obtained from a CT scan of a 43 years old woman (Figure 2 panel (a)). We used the medical imaging software Osirix to segment the heart and the torso from the CT scan DICOM files. We then construct the meshes using the CardioViz3D software.

Scenario. We suppose a circular ischemic region of radius 13 mm as is shown in Figure 2 panel (a). Contrary to the methodology used in Alvarez et al [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF], for this experiment we do not assume an ischemic initial guess. We start the iterative method assuming healthy conditions for the entire tissue. At each iteration, the functional cost gradient is computed and both parameters, τ in and v rest , are updated. After 28 iterations, stop criteria is accomplished obtaining the final reconstruction. For this scenario, we obtain a correlation coefficient between the reconstructed ischemic region and the ground truth of CC = 0.8241. In literature, we can find different methodologies to address the problem of locatize ischemic regions using potential measurements taken at the body surface. Many of them use voltage measurements in an instant of time, without taking into account the spatio-temporal correlation presented in EGM measures. In this paper, we used the method proposed by [START_REF] Álvarez | Shape reconstruction of cardiac ischemia from non-contact intracardiac recordings: A model study[END_REF] and tested on 2D geometries in order to reconstruct ischemic areas. This algorithm allows to solve the inverse problem including, into the formulation of level set, a regularization with the spatio-temporal correlation presented in the recordings.

The methodology was tested on different 3D geometries. First, the correct performance of the method was tested in a two concentric spheres geometry. Where it was found that the solution to this problem is not dependent on the initial guess for the inverse procedure. Subsequently, we analyze the behavior to a more realistic geometry. We solved the inverse problem considering a single ischemia in a real anatomical geometry. The quality of the results obtained for each of the proposed scenarios shows that the ischemic region has been located satisfactorily, obtaining a correlation coefficient exceeds the value of 0.95 for the spheres scenario and for the anatomical geometry is around 0.85. We think that this difference in the CC is due to the fact that the torso surface in the spherical case is much closer to the heart surface than it is for the anatomical geometry.

In future work, we aim to incorporate more realistic settings using an electrophysiology detailed model, with a physiological representation of the ischemia model. We will also investigate the performance of this method in localizing ischemia from clinical measurements. This assumes that the electrophysiological model that would be used to constrain the minimization problem is sufficiently accurate.

Fig. 1 .

 1 Fig. 1. First row: The image represents the spatial distribution of the ischemia: red ischemia, green border zone and blue health tissue. Second row represents initial guesses (τin configuration) for different simulation cases. Case (a): a healthy cardiac tissue, Case (b): initializing with an small ischemic region close to the real one and Case (c): considering an ischemic region elsewhere in the cardiac tissue. Third row shows the final reconstruction of the ischemic region

Fig. 2 .

 2 Fig. 2. Second Experiment: (a) A snapshot of the torso potential at time 170 ms. (b): Representations of the true ischemic region (τin configuration). (c): The reconstructed ischemia region.

Acknowledgments

This work was supported in part by Spanish MICINN grants TEC-2013-46067-R, FIS2013-41802-R and by Carlos III of Madrid University PIF Grant to Carlos E. Chavez.


Julien Bouyssier

Nejib Zemzemi

Jacques Henry

INVERSE PROBLEM IN ELECTROCARDOGRAPHY VIA THE FACTORIZATION METHOD OF BOUNDARY VALUE PROBLEMS

Keywords: Inverse problem, electrocardiography, Riccati equations, factorization method, boundary value problems, electrocardiagraphic imaging (ECGI)

We present a new mathematical approach for solving the inverse problem in electrocardiography. This approach is based on the factorization of boundary value problems method. In this paper we derive the mathematical equations and test this method on synthetical data generated on realistic heart and torso geometries using the state-of-the-art bidomain model in the heart coupled to the Laplace equation in the torso. We measure the accuracy of the inverse solution using spatial Relative Error (RE) and Correlation Coefficient (CC).

INTRODUCTION

ECGI is a new imaging technique that noninvasively images cardiac electrical activity on the heart surface. In ECGI, a multi-electrode vest records body-surface potential maps (BSPMs). Then, using BSPMs, geometrical information from CT-scans, a mathematical modelling and a numerical algorithm, we reconstruct electrical potentials, electrograms and activation maps on the heart surface [START_REF] Bokeriya | Hardwaresoftware system for noninvasive electrocardiographic heart examination based on inverse problem of electrocardiography[END_REF][START_REF] Ramanathan | Noninvasive electrocardiographic imaging for cardiac electrophysiology and arrhythmia[END_REF][START_REF] Liu | Noninvasive reconstruction of three-dimensional ventricular activation sequence from the inverse solution of distributed equivalent current density[END_REF]. In this work we develop a new method called the factorization of boundary value problems method [START_REF] Henry | Factorization of secondorder elliptic boundary value problems by dynamic programming[END_REF][START_REF] Jday | Complétion de données frontières : la méthode de plongement invariant[END_REF][START_REF] Ben Abda | Boundary data completion: the method of boundary value problem factorization[END_REF] based on the invariant embedding approach. The electrical potential in the ECGI inverse problem is governed by the Laplace equation with a homogeneous Neumann boundary condition and an extra Dirichlet boundary condition on the body surface. We decompose this problem into two boundary value problems and treat it as an optimal control problem. In [START_REF] Andrieux | Solving cauchy problems by minimizing an energy-like functional[END_REF], this optimal control problem is solved using a minimization procedure of a cost function, allowing to compute both electrical potential and its normal derivative over the heart surface. In our approach, the two sub-problems are embedded into a family of similar problems on subdomains bounded by a boundary moving from the torso skin to the heart surface. We can This work is sponsored from the french government managed by the National Research Agency under the program "Investissements d'Avenir" by the grant number ANR-10-IAHU-04.

define Neumman-Dirichlet and Dirichlet-Neumann operators on this moving boundary that satisfy Riccati equations and which need to be solved from the torso surface to the heart surface. Then we can define a matrix system whose resolution is equivalent to minimise the cost function of the optimal control problem. As this is an ill-posed Cauchy problem, a regularization term must be applied on the matrix system to get physiologically acceptable solution. In this paper, we first introduce principles of the factorization by deriving equations in the simpler 3D case of a cylinder. Then, we present equations derived in the most 3D generalized case. This is applied to a human torso geometry case where numerical simulations are performed by using synthetical data based on the ECG solver [START_REF] Zemzemi | Étude théorique et numérique de lactivité électrique du coeur: Applications aux électrocardiogrammes[END_REF]. We analyse the accuracy of the method by comparing reconstructed and synthetical solutions with RE and CC computations and illustrations of reconstructed potential fields.

METHODS

In this section we present the mathematical formulation of the factorization method. For the sake of clarity we first derive the equation in the 3D cylindrical case. Then, we present the mathematical formulation of the factorization method in the 3D general case.

Initial problem

For a given BSPs data T , the electrical potential u in the torso domain Ω satisfies the following system (P 0 ) :

       ∆u = 0 in Ω , Ω : cylinder u = 0 on Σ , Σ : lateral surface u = T on Γ T , Γ T : torso surface ∇u • n = 0 on Γ T , n : outward normal to Γ T (1)
The aim of solving the inverse problem is to find the trace of u on the on the heart surface Γ H .

Optimal control problem

Following [START_REF] Andrieux | Solving cauchy problems by minimizing an energy-like functional[END_REF], we decompose (P 0 ) into two auxiliary subproblems :

       ∆u 1 = 0 in Ω u 1 = 0 on Σ u 1 = T on Γ T ∇u 1 • n = η on Γ H and        ∆u 2 = 0 in Ω u 2 = 0 on Σ u 2 = τ on Γ H ∇u 2 • n = 0 on Γ T (2)
We then define the cost function as follows :

E(η, τ ) = Ω (∇u 1 (η) -∇u 2 (τ )) 2 (3)
When E reaches zero, we have u 1 = u 2 , and both of them satisfy (P 0 ). Our goal is then to minimise E. In the following paragraph, we present the factorization approach allowing to solve (P 0 ). This problem is ill-posed [START_REF] Ben Abda | Boundary data completion: the method of boundary value problem factorization[END_REF], we then add a regularization term to avoid the ill-posedness. The cost function becomes :

E(η, τ ) = Ω (∇u 1 (η) -∇u 2 (τ )) 2 + ǫ ΓH (η 2 + τ 2 ) (4)
where ǫ is the regularization parameter.

Principle of invariant embbeding

The main idea of this approach is to "transport" potential data from torso surface to the heart surface as illustrated in Figure 1. For this, boundary value problems ( 2) are embedded into a family of similar problems on subdomains Ω s . They are bounded by a moving boundary Γ s defined at x = s for x = 0 -→ x = a, where a is the length of the cylinder, x = 0 corresponding to the torso surface, x = a corresponding to the heart surface and s is the variable that describes the axis of evolution. At each position x = s, we impose a Neumann boundary condition

       ∆u 1 s = 0 in Ω s u 1 s = 0 on Σ s u 1 s = T on Γ T ∇u 1 s • n = α on Γ s (5)        ∆u 2 s = 0 in Ω s u 2 s = 0 on Σ s u 2 s = β on Γ s ∇u 2 s • n = 0 on Γ T . (6)
At each position x = s, we define on Γ s an affine mapping :

Neumann to Dirichlet α -→ u 1 s : u 1 s = Q(s)α + w 1 (s) (7)
and a linear mapping :

Dirichlet to Neumann β -→ ∂u 2 s ∂x : ∂u 2 s ∂x = P (s)β ( 8
)
where w is the residual associated to the operator Q. Folloowing [START_REF] Henry | Factorization of secondorder elliptic boundary value problems by dynamic programming[END_REF], we substitute equations ( 7) in ( 5) and ( 8) in ( 6) to obtain the following system :

       dP dx + P 2 = -∆ y , P (0) = 0 dQ dx -Q∆ y Q = I, Q(0) = 0 dw1 dx -Q∆ y w = 0, w(0) = T (9)
where I is the identity function, ∆ y is the laplacian operator defined on the section Γ s . Operators P and Q satisfy Riccati equations, whereas the residual term w satisfies an ordinary differential equation depending on the operator Q. All the initial conditions are known and are defined on the torso surface.

According to [START_REF] Ben Abda | Boundary data completion: the method of boundary value problem factorization[END_REF], the energy function E(η, τ ) could be rewritten using P , Q and w as follows :

E(η, τ ) = C + [η, τ ]A[η, τ ] T -2 < P (a)w(a), τ > + ǫ < τ, τ > +ǫ < η, η > (10
) where :

A = Q(a) + ǫI -Q(a)P (a) -P (a)Q(a) P (a) + ǫI (11)
and C = ΓT T ∂ n w is a constant independent of τ and η.

We compute the gradient of E over τ and η and find that the optimum of E is reached for τ and η satisfying :

(1 + ǫ)I -Q(a)P (a) τ = w(a) (1 + ǫ)I -P (a)Q(a) η = P (a)w(a) (12)
Since we are only interested in computing the electrical potential on the heart surface we only solve the first equation of (12) using the Matlab Cholesky solver.

Equations in the 3D generalized case

We name Γ 0 the torso skin, Γ 1 the heart skin and Ω the torso domain. We recall T the electrical potential on the torso surface. For 0 ≤ s ≤ 1, we define a familly surfaces Γ s included between the torso skin and the heart skin. This familly will represent the moving boundary. We also introduce Φ s ′ s , which is the application that maps the surface Γ s to the surface Γ s ′ , and V (x) = d ds Φ s 0 (x) which could be seen as a "velocity" term. We decompose V into a normal term α n and a tangential term V Γs of the surface Γ s . Following the same derivation process from the equation ( 5) to the equation ( 9), we can establish the variational formulation of P and Q operators and the residuals term w. We use the first order finite element method to discretize theses equations on Γ s . We denote (e j ) N 1 the finite element basis allowing to approximate P , Q and w by P h , Q h and w h . Following [START_REF] Bouyssier | 3d generalisation of the factorization of boudary values problem method : application to the electrocardiography inverse problem[END_REF] we obtain :

d ds (M P h ) + P T h M α P h +BP h + (BP h ) T = K α P h (0) = 0 (13) M dQ h ds - dM ds Q h + Q T h K α Q h -B T Q h -(B T Q h ) T = M α , Q h (0) = 0 (14) M dw h ds + Q T h K α w h -B T w h = 0, w h (0) = T ( 15
)
where M and M α are mass matrices, K α a stiffness and B is the matrix defined by B i,j =< V Γs .∇e i , e j >. The size of all the defined matrix is N × N where N is the number of nodes on the heart surface. Then, we solve ( 13)-( 15) using a first order explicit Euler scheme to obtain P h (1), Q h (1) and w h (1). Finally we obtain the heart potential by solving the linear system (12).

RESULTS

A computational geometry of a human torso volume was produced from a 43 years old women anatomical model. We used the bidomain model to simulate one second of heart beat on the heart and body surface. We extract the electrical potential on the torso boundary and use it as an input of the inverse problem. At each time step the torso potential represents the data T as in equation (15). Figure 2 shows a snapshot of the torso potential (left) and the corresponding heart potential (right).

In order to generate the sequence of surfaces morphing from Fig. 2. Example of torso and heart potential field generated by the bidomain model.

the torso surface to the heart surface, we solve the Laplace equation with a Dirichlet boundary condition equal to 1 on the torso surface and 0 on the heart surface as shown in Figure 3 (left). We then extract the surfaces as the iso-values of the solution as shown in Figure 3 (right). These iso-surfaces are Fig. 3. Illustration of the surfaces Γ s generated. Left : isovalues computed in the torso volume. Right : surfaces created from iso-values.

used to compute the operators P h (1), Q h (1) and w h (1) and construct the electrical potential on the heart surface as explained in section 2.3 for various ǫ values (10 -1 , 10 -2 , 10 -3 and 10 -4 ). In Figure 4, we show a comparison between the measured (left column) and the inverse solutions (right column) at times t = 100 ms (top) and t = 250 ms (bottom) for ǫ = 10 -4 . We remark that time t = 100 ms, the inverse solution is diffusive compared to the exact solution which has a sharp wave front. Whereas at time t = 250 ms where the exact solution is smoother and the reconstructed potential is more accurate. In Figure 5 we show the time course of the

DISCUSSION

This work represents the first use of the factorization of boundary value problems method in solving the inverse problem in electrocardiography. The method was firstly introduced for general Laplace equation in 2D [START_REF] Ben Abda | Boundary data completion: the method of boundary value problem factorization[END_REF] and it was technically challenging to implement this method in a 3D framework and adapt it to the ECGI problem. The most challenging parts are : 1) generating the family of moving surfaces between the heart and the torso, 2) solving the non linear Riccati equations on the moving surfaces. We used an explicit Euler scheme to linearize and solve the Riccati equations in order to build the operators P and Q at the heart surface. This scheme allows to conserve the positivity and the symmetry of the operator P , but does not necessarily conserve the symmetry of the operator Q. More adapted numerical schemes should be carried out in future works in order to accurately solve the operator Q.

CONCLUSION

In this paper, we presented a new approach for solving the inverse problem in electrocardiography called the factoriza-tion of boundary value problems method. We derived equations of the method in a cylindrical 3D case to illustrate this approach then in the general 3D case. We have tested this method on synthetical data generated on a realistic human geometry. Since the CC is relatively high, the method could be considered as a good candidate for applications like capturing the wave front position on the heart surface. We think that RE could be increased by improving the numerical discretization of the equation solving the operator Q. We will also take into account the time sequences by adding a temporal regularization term into the cost function. This would be subject of future works.

∂u 1 s

 1 ∂x |Γs = α for the u 1 problem and a Dirichlet boundary condition (u 2 s ) |Γs = β for the u 2 problem

Fig. 1 .

 1 Fig. 1. Illustration of the moving boundary for the cylinder case.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Measured (left column) and reconstructed (right column) potential on the heart surface at t = 100 ms (top) and 250 ms (bottom).relative error (RE) in space and in Figure6we show the time course of the correlation coefficient (CC), for the different values of the regularization parameter ǫ. Although the RE is high, the average CC is close to 0.7 allowing to a relatively

Fig. 6 .

 6 Fig. 6. Correlation coefficient computed for several time steps.