Martina Teichert

Charlotte Van Noord

Andre Uitterlinden

Albert Hofman

Peter N Buhre

Peter Agm De Smet

Sabine Straus

Bruno Hch Stricker
email: stricker@erasmusmc.nl

Loes E Visser

Andr÷ G Uitterlinden

B H Ch

MMed Stricker

Proton pump inhibitors and the risk of overanticoagulation during acenocoumarol maintenance treatment

Keywords: acenocoumarol, proton pump inhibitors, omeprazole, rabeprazole, pantoprazole, lansoprazole, esomeprazole

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

Abstract

Background

In the Netherlands, several reports described a potentiation of acenocoumarol-induced anticoagulation by co-medication of (es)omeprazole.

Objective

To investigate the effects of co-medication with PPIs on overanticoagulation during acenocoumarol maintenance treatment.

Methods

All 2,755 subjects from the Rotterdam Study who received acenocoumarol maintenance treatment between April 1 st, 1991 and September 9 th , 2009 were followed for events of an international normalized ratio (INR) ≥6, until death, end of treatment, or end of the study period. With the Andersen-Gill extension of the Cox proportional hazards model, risks for repeated events of overanticoagulation in relation to concomitant PPI use were calculated.

Results

The risk for overanticoagulation was most pronounced for esomeprazole (HR 1.99, 95% CI 1.55 -2.55) and lansoprazole (HR 1.49, 95% CI 1.05 -2.10). There was also a lower and non-significant risk increase for the other PPIs. We did not detect a modification of these results by CYP2C19*2 genotype.

Conclusion

Caution should be paid to co-medication with esomeprazole and lansoprazole during acenocoumarol treatment and possibly also with the other PPIs.

Introduction

Coumarin anticoagulants are first choice in treatment and prevention of arterial or venous thrombosis. [START_REF] Ansell | The pharmacology and management of the vitamin K antagonists: the Seventh ACCP Conference on Antithrombotic and Thrombolytic Therapy[END_REF] Warfarin, acenocoumarol and phenprocoumon act as vitamin-K antagonists by inhibiting the synthesis of coagulation factor II, VII, IX and X. Due to a particularly narrow therapeutic range, patients treated with these drugs have to be closely monitored by regular assessments of the international normalized ratio (INR) to warrant anticoagulation without serious bleedings. Individual dosing schemes differ widely between patients, mainly due to genetic variation in the vitamin K epoxide reductase complex subunit 1 (VKORC1) and cytochrome P450 2C9 (CYP2C9) genes. In addition, age, sex and body mass index (BMI) are important determining factors. [START_REF] Aquilante | Influence of coagulation factor, vitamin K epoxide reductase complex subunit 1, and cytochrome P450 2C9 gene polymorphisms on warfarin dose requirements[END_REF][START_REF] Carlquist | Genotypes of the cytochrome p450 isoform, CYP2C9, and the vitamin K epoxide reductase complex subunit 1 conjointly determine stable warfarin dose: a prospective study[END_REF][START_REF] D'andrea | A polymorphism in the VKORC1 gene is associated with an interindividual variability in the dose-anticoagulant effect of warfarin[END_REF][START_REF] Geisen | VKORC1 haplotypes and their impact on the inter-individual and inter-ethnical variability of oral anticoagulation[END_REF][START_REF] Oldenburg | Current pharmacogenetic developments in oral anticoagulation therapy: The influence of variant VKORC1 and CYP2C9 alleles[END_REF][START_REF] Qazim | Dependency of phenprocoumon dosage on polymorphisms in the VKORC1 and CYP2C9 genes[END_REF][START_REF] Sconce | The impact of CYP2C9 and VKORC1 genetic polymorphism and patient characteristics upon warfarin dose requirements: proposal for a new dosing regimen[END_REF][START_REF] Teichert | A genome-wide association study of acenocoumarol maintenance dosage[END_REF] Therefore, coumarin dosage is carefully titrated during an initiation period by reference to the INR to achieve a stable individual maintenance dosage. However, also within the maintenance period, effects of coumarins can vary due to interactions with co-medication, development of co-morbidity or changes in lifestyle. In the Netherlands, acenocoumarol and phenprocoumon are most frequently used whereas warfarin is only given on rare occasions. Although management of coumarin use in the Netherlands is constantly monitored in anticoagulation clinics by regular INR-measurements to obtain optimal INR-levels, bleeding associated with coumarins is among the leading causes of druginduced hospital admissions. [START_REF] Leendertse | Frequency of and risk factors for preventable medication-related hospital admissions in the Netherlands[END_REF][START_REF] Van Der Hooft | Adverse drug reaction-related hospitalisations: a population-based cohort study[END_REF] Over the years, different drugs were reported to increase the anticoagulation effect of coumarins. [START_REF] Blume | Pharmacokinetic drug interaction profiles of proton pump inhibitors[END_REF][START_REF] Wells | Interactions of warfarin with drugs and food[END_REF] four are increasingly associated with overanticoagulation, and above six the chance of serious bleedings strongly increases. [START_REF] Cannegieter | Optimal oral anticoagulant therapy in patients with mechanical heart valves[END_REF] Competitive inhibition of CYP2C19 has been suggested as a possible mechanism for interactions of PPIs with other drugs. [START_REF] Blume | Pharmacokinetic drug interaction profiles of proton pump inhibitors[END_REF] Although earlier studies did not show changes in acenocoumarol pharmacokinetics or pharmacodynamics [START_REF] Hoorn De | No effect of short-term omeprazole intake on acenocoumarol pharmacokinetics and pharmacodynamics[END_REF] or a need of acenocoumarol dosage adjustment [START_REF] Vreeburg | Lack of effect of omeprazole in oral acenocoumarol anticoagulant therapy[END_REF]

Results

Of the 14,926 subjects in the Rotterdam cohorts, 2,755 had used acenocoumarol during the study period for longer than 42 days continuously as maintenance therapy. In 887 subjects an (HWE,). Subjects homozygous with the variant A-allele of CYP2C19*2 had a decreased risk on overanticoagulation during acenocoumarol treatment (HR 0.54; 95% CI 0.36-0.83). As this group was quite small (2.4%), we pooled subjects homozygous and heterozygous for the variant A-allele. Any variant A-allele was then no longer associated with the risk of overanticoagulation (HR 0.99; 95% CI 0.89-1.11).

In our cohort in total 457 (16.6%) subjects used omeprazole, 160 (5.8%) pantoprazole, 102 (3.7%) lansoprazole, 286 (10.4%) rabeprazole and 149 (5.4%) esomeprazole (table 2). In 223 events of an INR ≥6 a PPI (10.4% of all events) was used on the index date for at least three days. Co-medication during acenocoumarol maintenance treatment with esomeprazole was associated with a doubled risk of an INR ≥6 (HR 1.99; 95% CI 1.55 -2.55) and comedication with lansoprazole increased the risk by 49 percent (HR 1.49; 95% CI 1.05 -2.10, table 2). The other PPIs also increased the risk of overanticoagulation less pronounced between 1.12 and 1.23 times, but associations were just below the significance threshold after adjustment for co-medication with NSAIDs. Co-medication with ranitidine, used by 229 subjects (10.5%), showed no association with overanticoagulation (HR 1.06, 95% CI 0.83 -1.35).

For esomeprazole and lansoprazole we investigated whether the CYP2C19 genotype modified the interaction with acenocoumarol (table 3). We did not find a multiplicative

Discussion

In our study population, all PPIs tended to increase the risk of overanticoagulation during acenocoumarol maintenance treatment. This was most pronounced for esomeprazole, doubling the risk of overanticoagulation and for lansoprazole with a risk increase of 49%.

Risk increases for omeprazole, pantoprazole and rabeprazole were between 12 and 23% and just failed to reach significance after adjustment for NSAID co-medication. PPIs have been mentioned occasionally as potential risk factors for overanticoagulation during treatment with vitamin K antagonist in case reports, however case control studies and trials concluded no effects of clinical relevance. [START_REF] Ehrlich | Lack of pharmacodynamic and pharmacokinetic interaction between pantoprazole and phenprocoumon in man[END_REF][START_REF] Hoorn De | No effect of short-term omeprazole intake on acenocoumarol pharmacokinetics and pharmacodynamics[END_REF][START_REF] Sutfin | Stereoselective interaction of omeprazole with warfarin in healthy men[END_REF][START_REF] Unge | A study of the interaction of omeprazole and warfarin in anticoagulated patients[END_REF][START_REF] Vreeburg | Lack of effect of omeprazole in oral acenocoumarol anticoagulant therapy[END_REF] As a possible mechanism a competitive inhibition of CYP2C19 by PPIs on coumarin clearance has been mentioned as both drug groups are metabolized by this enzyme and the INR increases in the case reports were reported within 11 days after start of the PPI. [START_REF] Blume | Pharmacokinetic drug interaction profiles of proton pump inhibitors[END_REF][START_REF] Gerson | Proton pump inhibitors and their drug interactions: an evidence-based approach[END_REF][START_REF] Mccarthy | Impact of cotherapy with some proton pump inhibitors on medical claims among HMO patients already using other common drugs also cleared by cytochrome P450[END_REF] (S-)acenocoumarol, the enantiomeric form being mainly responsible for the drug effect, is for about 20% metabolized via CYP2C19. For the five PPIs, in vitro studies showed that all five PPIs were able to competitively inhibit CYP2C19, but to a different extent. Lansoprazole had the highest inhibitory potency of CYP2C19 in vitro. [START_REF] Li | Comparison of inhibitory effects of the proton pump-inhibiting drugs omeprazole, esomeprazole, lansoprazole, pantoprazole, and rabeprazole on human cytochrome P450 activities[END_REF] Rabeprazole had a relatively lower in vitro inhibition of CYP2C19 than the other PPIs, possibly as it is metabolized mainly via a non-enzymatic reaction to a thioether compound with only minor CYP2C19 involvement. [START_REF] Gerson | Proton pump inhibitors and their drug interactions: an evidence-based approach[END_REF][START_REF] Hata | Proton pump inhibitors may increase the risk of delayed bleeding complications after open heart surgery if used concomitantly with warfarin[END_REF][START_REF] Li | Comparison of inhibitory effects of the proton pump-inhibiting drugs omeprazole, esomeprazole, lansoprazole, pantoprazole, and rabeprazole on human cytochrome P450 activities[END_REF] Omeprazole and esomeprazole in vitro also reduced CYP2C19 activity [START_REF] Li | Comparison of inhibitory effects of the proton pump-inhibiting drugs omeprazole, esomeprazole, lansoprazole, pantoprazole, and rabeprazole on human cytochrome P450 activities[END_REF] and were more potent than pantoprazole. [START_REF] Simon | Faster in vitro biotransformation of S-omeprazole by the cytochrome P450 isoenzyme system compared to pantoprazole[END_REF] Theoretically a competition for scarce CYP2C19 enzyme might prolong the effectiveness of the coumarins which would result in an increased INR. In subjects with variant allele genotypes and decreased enzyme availability this effect may be even more pronounced. However, from our analysis, no multiplicative or additive effect modification of CYP2C19*2 variant alleles on the association between esomeprazole or lansoprazole and acenocoumarol was detected. In our study population there were no subjects with the CYP2C19*2 AA genotype using esomeprazole or lansoprazole and this limited our possibilities for stratification. Possibly the effect of the variant allele on acenocoumarol elimination was too small and more cases are needed to warrant enough power for a genotype-stratified analysis.

Our study is the first observational cohort study of an increase of acenocoumarol effectiveness by co-medication of all PPIs separately with a substantial number of subjects during the whole period of acenocoumarol maintenance treatment within a subject. We consider the chance of bias and confounding due to study design as negligible. First, the Rotterdam Study is a prospective cohort study and the regional anticoagulation clinic covered a complete area of more than one million inhabitants in the Rotterdam area. Consequently, everyone who is treated with a coumarin anticoagulant as an outpatient will be registered as such and selection bias is unlikely. Second, all medication use of all subjects was almost completely covered by the pharmacy data we retrieved and during the study period PPIs were only available on prescription via pharmacies. Any lack of compliance would move our results into the direction of the null hypothesis, and would tend to make our results conservative. During our study period, an association between PPI and acenocoumarol overanticoagulation was not yet common knowledge. Third we followed the cohort members during their whole period of acenocoumarol maintenance therapy and compared cases to all other cohort members available at the index dates for PPI exposure. The analysis made use of all data available and adjusted for time varying effects. However, in observational studies there is always a risk of confounding by indication. We therefore performed the same analysis for rantidine, a H2-receptor antagonist used for the same indication but not known for inhibitory potential of CYP2C19. We did not find an association for ranitidine with an increased risk on achieving INR ≥6 during acenocoumarol treatment. As PPIs might have been started to treat the symptoms of an increased INR, in our analysis use of PPIs -and also ranitidine -had to be started for at least 3 days prior to the index date. We did not take bleedings as an outcome because PPI s are likely to be started with complaints preceding gastrointestinal bleeding. Associations found might be due to protopathic bias.

In conclusion, in this population-based cohort study among outpatients of an anticoagulation clinic using acenocoumarol for maintenance treatment, esomeprazole doubled the risk on an INR ≥6 and lansoprazole tended to increase this risk by approximately 50%.

For omeprazole, pantoprazole and rabeprazole, a risk increase was less pronounced and nonsignificant. Our results suggest that extra monitoring during acenocoumarol treatment may be warranted in patients on esomeprazole or lansoprazole and this perhaps applies to other PPIs as well.

Methods

Setting

We selected all subjects from the three cohorts of the Rotterdam Study (RS-I , RS-II and RS-III). The rationale and design of the Rotterdam Study have been described elsewhere. [START_REF] Hofman | The Rotterdam Study: 2010 objectives and design update[END_REF][START_REF] Hofman | The Rotterdam Study: objectives and design update[END_REF][START_REF] Hofman | Determinants of disease and disability in the elderly: the Rotterdam Elderly Study[END_REF] Coumarin doses are adjusted on the basis of computerized dose calculations. More than 99% of participants fill their drug prescriptions at seven regional pharmacies, which are fully computerised. Complete data on drug use from these pharmacies were available as of January 1st , 1991. The pharmacy data include the Anatomical Therapeutical Chemical (ATC)-

Cohort and outcome definition

Our study population consisted of all patients of the three RS cohorts who started with acenocoumarol in the study period from April 1 st , 1991 through September 9 th , 2009 and used it consecutively for at least 42 days. The start date of April 1 st was chosen to ensure that at least 3 months of medication history from the pharmacy was available for each cohort member. We regarded the period starting 42 days after initiation with acenocoumarol as maintenance period. In general, steady state of a drug is usually achieved within 5-7 half-lives of drug elimination. For (R-)acenocoumarol, the enantiomeric form on which treatment effects mainly depend, with a half-life of 8 hours, a period of 6 weeks taken as initiation period was considered extensive enough to reach a steady state. [START_REF] Flockhart | Pharmacogenetic testing of CYP2C9 and VKORC1 alleles for warfarin[END_REF] In patients from the Rotterdam study using acenocoumarol during the maintenance period, we took an event of an INR of 6 and greater after baseline study enrolment as an outcome. INR levels ≥6 are associated with an extensively increased risk of bleedings. [START_REF] Cannegieter | Optimal oral anticoagulant therapy in patients with mechanical heart valves[END_REF] We excluded INR events of six and greater that happened within 21 days of an earlier event.

Cofactors

The following baseline patient characteristics were considered as potential confounders or effect modifiers: sex, age, BMI, and target INR level. BMI was defined as (kg/m 2 ) and missing values were imputed with a linear regression model consisting of INR≥6, age, sex and target INR as variables. We further adjusted for oral non-steroidal anti-inflammatory drugs (NSAIDs) as co-medication (ATC code M01A), since PPIs might have been started because of treatment with an NSAID and NSAIDs have been shown earlier to increase acenocoumarol effectiveness. [START_REF] Visser | Allelic Variants of Cytochrome P450 2C9 Modify the Interaction Between Nonsteroidal Anti-inflammatory Drugs and Coumarin Anticoagulants[END_REF] In a subanalysis we studied the association of acenocoumarol maintenance treatment with use of ranitidine (A02BA02). This drug is used for the same indication and not known to (Hancock, et al).

Suggestion to skip the following details and just to keep the reference of Hancock et al. 449) were genotyped with 550k (V.3) single and duo chips, while RS-II participants (n=2,516) were genotyped with 550k (V.3) duo and 610k Quad chips. RS-III (n=2.420) participants were genotyped with the Human 610 Quad Arrays of Illumina.

Genotype calling was performed in RS-I using BeadStudio software (version 0.3.10.14), GenomesStudio in in RSIII. Participants with call rates <97.5%, excess autosomal heterozygosity, sex mismatch or outlying identity-by-state clustering estimates were excluded. After quality control, 5,974 RS-I participants, 2,157 RS-II and 2,078 RS-III participants remained with complete data on genotyping. (Hancock, et al) For imputation, 512,349 autosomal SNPs in RS-I and 466,389 autosomal SNPs in RS-II and RS-III were used after exclusions for call rate < 98%, HWE P<10 -6 , and MAF <1%, in MACH (version 1.00.15 for RS-I, 1.00.16 for RS-II and RS-III) with reference to the 2,543,886 SNPs of the HapMap CEU (release 22, build 36). (Hancock, et al).

Skip until here.

Statistical analysis

Within a subject, an INR ≥6 could occur more than once and co-medication with PPIs could change during acenocoumarol treatment period. In order to include all information available for the whole study period, we used the Andersen-Gill extension of the Cox proportional b For the CY2C19*2 genotype within total cohort, allelic frequency of the G-allele was 84.3%, of the A-allele 15.7%, and HWE was 0.93 Abbreviations: INR, international normalized ratio; HR, hazard ratio; CI, confidence interval; HWE, Hardy-Weinberg equilibrium.

F o r P e e r R e v i e w

  Until April 2009, the Netherlands' national pharmacovigilance center (LAREB) received nine reports about an increase of the anticoagulant effect of acenocoumarol within eleven days after start of treatment with omeprazole and seven reports in combination with esomeprazole.(Lareb Dutch National ) Cases of coumarin potentiation by other PPIs have not been reported. In four of these reports, the INR rose above six. INR measurements rising above

  , the Dutch Federation of Anticoagulation Clinics added esomeprazole and omeprazole to the list of drugs potentially interacting with coumarins as of January 2010.(Federatie van Nederlandse Trombosediensten visited 28th September 2010) Before that date, Dutch computerized medication surveillance systems of GPs and pharmacists did not flag interactions between acenocoumarol and PPIs. In a large prospective population-based cohort study, we investigated whether comedication with the PPIs omeprazole, pantoprazole, lansoprazole, rabeprazole or esomeprazole was associated with an increased risk of overanticoagulation during acenocoumarol maintenance treatment and whether an effect was modified by CYP2C19 variant alleles.

  drugs separately or grouped for use of either of the drugs and at least one CYP2C19*2 variant allele. To investigate presence of an additive effect modification we formed 4 groups for the combinations of drug use yes and no with absence and presence of CYP2C19*2 variant alleles. No additive effect modification was detected with this dominant model by genetic variation in the CYP2C19 gene as the risk for overanticoagulation for esomeprazole-or lansoprazole did not differ significantly between subjects within strata of CYP2C19*2 genotypes (p-value = 0.11 for the Synergie Index for additive effectmodification, SI= 1.67 (95% CI 0.78 -3.57) (Hosmer and Lemeshow 1992)). In our population there were no users of esomeprazole or lansoprazole with the CYP2C19*2 AA genotype and therefore it was not possible to investigate additive effect modification by genetic variation in the CYP2C19 gene with a recessive model.

  However, the trials were only short term trials with small numbers of subjects and the case control study took dosage change as an outcome without regard to INR measurements. The highest number of case reports and .docC:\ScholarOne\conversions\9255027-5678283\$ASQ9255027_File000001_160903508.doc pharmacokinetic studies are available for omeprazole in combination with coumarins. Two case reports described an INR increase in users of phenprocoumon after initiation of omeprazole and return to normal after its cessation.[START_REF] Enderle | Drug interaction: omeprazole and phenprocoumon[END_REF] In one study for comedication with lansoprazole or rabeprazole in patients after open heart surgery,[START_REF] Hata | Proton pump inhibitors may increase the risk of delayed bleeding complications after open heart surgery if used concomitantly with warfarin[END_REF] lansoprazole enhanced the anticoagulation effects of warfarin whereas rabeprazole could be used concomitantly without increasing the risk of overanticoagulation. We are not aware of a clinical study on the effect of esomeprazole on coumarins to support our findings. Our results show a more pronounced association for esomeprazole than for omeprazole. Omeprazole is a racemic composition of its two optical isomers, (S-)omeprazole (esomeprazole) and (R-)omeprazole.(Andersson, et al 2001) The clearance of the two enantiomeric omeprazole-forms by CYP2C19 is stereo-selective. In vivo, CYP2C19 was responsible for 90% of the metabolism of (R-)omeprazole and 70% of (S-) omeprazole.(Aebeloe, et al 2000) In our study population, omeprazole and esomeprazole were used in similar defined daily doses (DDD), omeprazole at a mean DDD of 1.35 (range 0.21 -4.50) and esomeprazole at a mean DDD of 1.23 (range 0.33 -4.00). However, these DDDs do not represent equipotent dosages as 1 DDD omeprazole equals 20 mg and 1 DDD esomeprazole 30 mg. Consequently in our study population, esomeprazole was used in a much higher pharmacological dosage than omeprazole and this may explain a higher influence on acenocoumarol effectiveness.

  Page 9 of 19 C:\ScholarOne\conversions\9255027-5678283\$ASQ9255027_File000001_160903508.docC:\ScholarOne\conversions\9255027-5678283\$ASQ9255027_File000001_160903508.doc

  In brief, the Rotterdam Study is a prospective population-based cohort study, designed to study neurological, cardiovascular, locomotor and ophthalmologic diseases in a population of people of 45 years and older. During different periods, eligible inhabitants of Ommoord, a suburb of Rotterdam, were invited to participate. The RS-I cohort consisted of 7,983 subjects (response rate 78%), the RS-II cohort of 3,011 participants (response rate 67%) and the RS-III cohort of 3,932 subjects (response rate 65%). The RS-I cohort had baseline examinations during 1990 -1993 with completion of standardized questionnaires, blood sampling and DNA isolation. The RS-II was formed as an independent cohort in 1999 with baseline examinations between 2000 and 2001 and the RS-III cohort was examined between 2006 and 2008. A regional anticoagulation clinic, Star Medical Diagnostic Center, monitors all inhabitants of Ommoord with an indication for anticoagulant therapy. This clinic covers the patients from all three RS cohorts. From this clinic, since 1984, all data on dosing, laboratory and clinical data, including data on bleeding complications are fully computerized. The patients own treating physician decides about the type of anticoagulant. Prothrombin times are monitored every 1-6 weeks, depending on the target level and stability of the INR.

  Health Organization Collaborating Centre for Drug Statistics Methodology 2004), the filling date, the total amount of drug units per prescription, the prescribed daily number of units, and product name of the drugs. As PPIs we included omeprazole (selected from pharmacy dispensing data by the ATC code of the WHO (World Health Organization Collaborating Centre for Drug Statistics Methodology 2004), A02BC01), pantoprazole (A02BC02), lansoprazole (A02BC03), rabeprazole (A02BC04) and esomeprazole (A02BC05).

  be dependent on biotransformation via CYP2C19. (Koninklijke Nederlandse Maatschappij ter bevordering der Pharmacie; The Hague 2008) In further subanalyses, we studied effect modification of the CYP2C19*2 variant allele (rs4244285) on an INR ≥6 during use of acenocoumarol. To date, about 19 variant alleles of CYP2C19 have been identified.(Yang and Link 2010) The majority of individuals can be classified into three phenotype groups, homozygous and heterozygous extensive metabolizers and poor metabolizers. They are determined by the wild-type genotype and two mutated alleles, CYP2C19*2 and CYP2C19*3.The principal defect, CYP2C19*2 (rs4244285) with a G-allele for normal function and an Aallele as variant with decreased efficiency, leads to a truncated protein and was genotyped in our population. CYP2C19*3 is rare in Caucasians and was not genotyped in the Rotterdam Study.[START_REF] Cannegieter | Optimal oral anticoagulant therapy in patients with mechanical heart valves[END_REF] with available DNA were genotyped using Illumina Infinium II HumanHap BeadChips at the Department of Internal Medicine, Erasmus Medical Center following manufacturer's protocols. With this technology DNA is methylated and the assay can measure the methylation level at 14,495 genes.

  This model allows to study multiple events of an INR ≥6 within one subject and in association with PPI use as a time-varying covariable.[START_REF] Andersen | Cox's regression model for counting processes: a large sample study[END_REF] All subjects on acenocoumarol maintenance therapy were followed as of April 1 st , 1991, from their first INR assessment until the last INR assessment because of the end of their treatment, last INR ≥6 event or the end of the study period, whichever came first. The date on which an INR ≥6 occurred was taken as the index date. In order to exclude protopathic bias, patients were only considered as exposed to PPIs at the index date if they had started PPI at least 3 days before that date. Each case was compared for PPI exposure to all subjects who were on acenocoumarol maintenance treatment at the index date.[START_REF] Stricker | Analysis of individual drug use as a time-varying determinant of exposure in prospective population-based cohort studies[END_REF] Thus cases could serve as controls on other index dates when still being on acenocoumarol treatment. We computed hazard ratios (HRs) and their 95% confidence intervals (CIs) for all events of an INR≥6. Risk estimates were adjusted for age, sex, target INR, BMI and NSAID co-medication. To study effect modification by CYP2C19 genotype, patients were stratified according to their genotype as CYP2C19*2 homozygous G-alleles and variant type (CYP2C19*2 heterozygous G/A-alleles or CYP2C19*2 homozygous A-alleles). SPSS 15.0 was used for data management and SAS 9.20 for the Andersen-Gill analysis.

  management of GWAS genotype data for the Rotterdam Study is supported by the Netherlands Organisation of Scientific Research NWO Investments (nr. 175.010.2005.011, 911-03-012). This study is funded by the Research Institute for Diseases in the Elderly (014-93-015; RIDE2) and the Netherlands Genomics Initiative (NGI)/Netherlands Organisation for Scientific Research (NWO) project nr. 050-060-810. We thank Pascal Arp, Mila Jhamai, Dr Michael Moorhouse, Marijn Verkerk and Sander Bervoets for their help in creating the GWAS database. The authors are grateful to the study participants, the staff from the Rotterdam Study and the participating general practitioners and pharmacists.

Table 1 Characteristics of acenocoumarol users with INR≥6.0 and total cohort

 1

		Patients with INR>=6.0 (N=887)	Total cohort (N=2,755)	HR a	95% CI
	Gender
	Male (%)	405 (45.7)	1,192 (43.3)	1.00	Reference
	Female (%)	482 (54.3)	1,563 (56.7)	1.25	(1.15 -1.36)
	Start age (years)
	45 -54	10(1.1%)	114 (4.2%)	1.00	Reference
	55 -64	72 (8.1%)	383 (14.0%)	2.89	1.75 -4.78
	65 -74	214 (24.1%)	907 (33.2%)	2.74	1.69 -4.45
	75 -84	389 (43.9%)	998 (36.5%)	4.14	2.56 -6.68
	>85	202 (22.8%)	334 (12.2%)	6.10	2.63 -14.1
	BMI (kg/m 2 )
	15.0 -20.0	17 (1.9)	45 (1.6)	1.00	Reference
	20.1 -25.0	228 (25.7)	718 (26.1)	0.70	0.52 -0.93
	25.1 -30.0	492 (55.5)	1446 (52.6)	0.68	0.51 -0.91
	>30.0	150 (16.9)	539 (19.6)	0.56	0.42 -0.76
	Target level (INR)
	2.0 -2.5	1 (0.1)	118 (4.3)	1.00	Reference
	2.0 -3.5	389 (43.9%)	1728 (62.7%)	10.1	3.
	3.0 -4.0	482 (54.3%)	889 (32.3%)	27.9
	3.5 -4.5	15 (1.7%)	20 (0.7%)	65.4

23 -31.2 9.00 -86.7 17.3 -247

  Statistically significant values are printed in bold. a Univariate analysis of RR were performed with a n Andersen-Gill model. RRs cannot be calculated with the numbers in this table because controls may later become cases; significant values are printed in bold.

	Subjects with NSAID use	89 (10.0)	470 (17.0)	1.31	1.15 -1.50
	CYP2C19*2 genotypes b	N=648 (73.1% of all cases)	N=2059 (74.7% of the total group)
	GG	472 (72.8%)	1464 (71.1%)	1.00	Reference
	GA and AA	176 (27.2%)	595 (28.9%)	0.99	0.89 -1.11

Table 2 Association between overanticoagulation (INR≥6) under acenocoumarol maintenance treatment and PPIs

 2 Statistically significant values are printed in bold a In this time-dependent analysis, exposure in case patient and in the rest of the cohort is assessed at the time of the outcome in each case patient (index date). As control patients can be used multiple times , the number of assessments in the reference group is much larger than the number of individuals. Hence, crude RRs cannot be calculated from the data in this table. b adjusted for age, sex, BMI and target INR c adjusted for age, sex, BMI, target INR and use of non-steroidal anti-inflammatory drugs, ATC-group M01A.

	Andersen-Gill analysis

British Journal of Haematology

Acknowledgements

The Rotterdam Study is funded by Erasmus Medical Center and Erasmus University, Rotterdam, Netherlands Organization for the Health Research and Development (ZonMw), the Research Institute for Diseases in the Elderly (RIDE), the Ministry of Education, Culture and Science, the Ministry for Health, Welfare and Sports, the European Commission (DG XII), and the Municipality of Rotterdam.

Conflict of interest statement

All authors state that they do not have any conflicts of interest to declare.

F o r P e e r R e v i e w


The degrees of freedom of the Lasso in underdetermined linear regression models

Gabriel Peyré

Ceremade Université Paris-Dauphine

Abstract-In this paper, we investigate the degrees of freedom (df) of penalized ℓ 1 minimization (also known as the Lasso) for an underdetermined linear regression model. We show that under a suitable condition on the design matrix, the number of nonzero coefficients of the Lasso solution is an unbiased estimate for the degrees of freedom. An effective estimator of the number of degrees of freedom may have several applications including an objectively guided choice of the regularization parameter in the Lasso through the SURE or GCV frameworks.

Index Terms-Lasso, degrees of freedom, SURE.

I. INTRODUCTION

Consider the following linear regression model

y = Ax 0 + ε, (1)
where y ∈ R n is the response vector, A = (a1, • • • , ap) ∈ R n×p is a deterministic design matrix with n < p, x 0 ∈ R p is the unknown regression vector, and ε ∈ R n is the noise vector whose entries are i.i.d. N (0, σ 2 ). The goal is to solve (1) when the solution is assumed to be sparse. Towards this goal, a now popular estimator is the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF]. The Lasso estimate amounts to solving the following convex problem min

x∈R p 1 2 y -Ax 2 2 + λ x 1, (2)
where λ > 0 is the regularization or tuning parameter. In the last years, there has been a huge amount of work where efforts have focused on investigating the theoretical guarantees of the Lasso as a sparse recovery procedure from noisy measurements in the underdetermined case n < p. Degrees of freedom df is a familiar phrase in statistics. In overdetermined linear regression df is the number of estimated predictors. Degrees of freedom is often used to quantify the model complexity of a statistical modeling procedure (e.g. it corresponds to the penalty term of model selection criteria such as AIC and BIC). However, generally speaking, there is no exact correspondence between the degrees of freedom df and the number of parameters in the model. On the other hand, the degrees of freedom plays an important role for an objective selection of the tuning parameter.

Let us denote by x any estimator of x 0 which depends on y and let ŷ = Ax. Since y ∼ N (Ax 0 , σ 2 ), according to [START_REF] Efron | How biased is the apparent error rate of a prediction rule[END_REF], the degrees of freedom of ŷ is

df (ŷ) = n X i=1 cov(ŷi, yi) σ 2 . ( 3
)
If ŷ is almost differentiable, Stein's lemma [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] yields the following unbiased estimator of df

df (ŷ) = div ŷ = n X i=1 ∂ ŷi ∂yi . (4)
Contributions Let μλ = μλ (y) = Ax λ (y) be the Lasso response vector, where xλ (y) is a solution of the Lasso problem (2). In the overdetermined case, i.e. n > p, xλ (y) is unique, and the authors in [START_REF] Zou | On the "degrees of freedom" of the Lasso[END_REF] showed that for any given λ the number of non-zero coefficients of xλ is an unbiased estimator of the degrees of the freedom of the Lasso. Though their proof contains a gap. The contribution of this paper is to extend their result to the underdetermined case where the Lasso solution is not unique. To ensure the uniqueness of the solution, we introduce the condition (UC) on the design matrix.

II. MAIN RESULTS

Let z ∈ R p , S ⊆ {1, 2, • • • , p} and |S| its cardinality. We denote by AS the submatrix AS = [• • • , aj, • • • ]j∈S, where aj is the jth column of A and the pseudo-inverse

(A t S AS) -1 A t S of AS is denoted A + S .
Let zj be the jth component of z. Similarly, we define zS = (• • • , zj, • • • )j∈S for z. Let supp(z) = {j : zj = 0} be the support or the active set of z.

Definition 1 (Condition (UC) [START_REF] Dossal | A necessary and sufficient condition for exact recovery by l1 minimization[END_REF]): A matrix A satisfies condition (UC) if, for all subsets I ⊂ {1, • • • , p} with |I| ≤ n, such that (ai)i∈I are linearly independent, for all indices j ∈ I and all vectors

V ∈ {-1, 1} |I| , | aj, (A + I ) t V | = 1. (5)
Theorem 1: Suppose that A satisfies condition (UC). For any y ∈ R n , there exists a finite set of values λ, denoted by {λm}, for which we have max

j ∈I | aj, y -Ax λ (y) | = λ, (6)
where I = supp(x(λ)) and xλ (y) is the solution of the Lasso. Furthermore, if λ ∈ ˜0, A t y ∞ ˆ\ {λm}, then

max j ∈I | aj, y -Ax λ (y) | < λ. (7)
Theorem 2: Suppose that A satisfies condition (UC). For any y ∈ R n , and all values of λ for which (7) is satisfied, we have

That is, using Stein's lemma [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF] and the divergence formula (8), the number of non-zero coefficients of xλ is an unbiased estimator of the degrees of the freedom of the Lasso.

•

  The Lasso response μλ (y) = Ax λ (y) is a uniformly Lipschitz function of y; • The support and vector sign of the Lasso solution are locally constant with respect to y, and consequently div μλ (y) = | supp(x λ (y))|.


L A Inza

J I Mars

J P Métaxian

G S O'brien

O Macedo

Seismo-Volcano Source Localization with Triaxial-Broadband Seismic Array

Keywords: Spatial analysis, Volcan monitoring, Fourier analysis, Time series analysis, Seismic array analysis

Seismo-volcano source localization is essential to improve our understanding of eruptive dynamics and of magmatic systems. The lack of clear seismic wave phases prohibits the use of classical location methods. Seismic antennas composed of one-component (1C) seismometers provide a good estimate of the back-azimuth of the wavefield. The depth estimation, on the other hand, is difficult or impossible to determine. As in classical seismology, the use of three component (3C) seismometers is now common in volcano studies. In order to determine the source location parameters (back-azimuth and depth), we extend the 1C seismic antenna approach to 3C's. This article discusses a high-resolution location method using a 3C array survey (3C MUSIC algorithm) with data from two seismic antennas installed on an andesitic volcano in Peru (Ubinas volcano). One of the main scientific questions related to the eruptive process of Ubinas volcano is the relationship between the magmatic explosive and LP swarms. After introducing the 3C array theory, we evaluate the robustness of the location method on a full wavefield 3D synthetic dataset generated using a digital elevation model of Ubinas volcano and a heterogeneous velocity model obtained from a tomography study. Results show that the back-azimuth determined using the 3C array has a smaller error than a 1C array. Only the 3C method allows the recovery of the source depths. Finally, we applied the 3C approach to a seismic event recorded in 2009. Crossing the estimated back-azimuth and incidence angles, we find a source located 1000 ±250 m below the bottom the active crater. Therefore, extending 1C arrays to 3C arrays in volcano monitoring allows a more accurate determination of the source epicenter and now an estimate for the depth.

INTRODUCTION

Source location (back-azimuth and depth determination) is a fundamental goal in volcano monitoring. Long-period (LP) events and tremor, which are directly related to magma ascent, constitute the main classes of seismic events observed on andesitic volcanoes [START_REF] Chouet | Long-period volcano seismicity: Its source and use in eruption forecasting[END_REF]. Locating these events is therefore necessary to better understand the eruptive dynamics and to improve the knowledge of the magmatic system. The lack of clear body-wave phase arrivals and emergent onsets in LP events and tremor makes locating these events extremely difficult using classical hypocenter determination methods based on phase picking and calculation of travel times. Other location methods have been used in recent years including a method based on the spatial distribution of seismic amplitudes to locate eruptive tremor sources on the Piton de la Fournaise volcano [START_REF] Battaglia | Location of seismic events and eruptive fissures on the piton de la fournaise volcano using seismic amplitudes[END_REF] and to track tremors produced by lahars on Cotopaxi volcano [START_REF] Hiroyuki | Seismic tracking of lahars using tremor signals[END_REF]. Dense one-component (1C) array methods based on time delays between close sensors, hence an estimation of the slowness vector of the wavefronts propagating across the array, have been used by several authors and applied to a great variety of volcanoes. For example [START_REF] Saccorotti | A probabilistic approach to the inversion of data from a seismic array and its application to volcanic signals[END_REF] and [START_REF] Rocca | Probabilistic source location of explosion quakes at stromboli volcano estimated with double array data[END_REF] applied dense onecomponent array techniques to locate explosive activity at Stromboli. [START_REF] Almendros | Mapping the sources of the seismic wave field at kilauea volcano, hawaii, using data recorded on multiple seismic antennas[END_REF] characterized the spatial extent of a hydrothermal system at Kilauea volcano by using similar techniques. [START_REF] Metaxian | Locating sources of volcanic tremor and emergent events by seismic triangulation[END_REF] used several small dense arrays to locate LP events and tremor sources at Arenal volcano and (Di Lieto et al., 2007) used two dense one-component arrays of short period seismometers to track volcanic tremor at Etna. Dense onecomponent array methods usually consist of only the vertical seismic component and allow a good estimation of the back-azimuth of the wavefield. Unfortunately, the depth estimation is poorly resolved because the incidence angle is very difficult to determine. To overcome this problem, triaxial sensors (3C) can be used. In this work we focus on the gain in back-azimuth resolution and depth determination obtained by three-component rather than onecomponent seismometers. In our study two 3C arrays were installed on Ubinas volcano, Peru, in March 2009 in order to determine the back-azimuth and depth of the seismo-volcano sources. The two small-aperture cross-shaped seismic arrays consisted of 12 threecomponent broadband seismometers. One of the main scientific questions related to the eruptive process at Ubinas is the relationship between the magmatic explosions and the LP swarms preceding these events by several tens of minutes to a few hours [START_REF] Macedo | Seismicity associated with the 2006-2008 eruption ubinas volcano[END_REF]. Source location of the LP events with a higher resolution compared to the one component dense array methods, particularly in the depth determination, is the main objective of this work. To achieve this goal and before working with real wavefield data we performed several numerical simulations of seismic wave propagation using a 3D digital elevation map and heterogeneous velocity model determined from a tomography study, [START_REF] Monteiller | An efficient algorithm for double-difference tomography and location in heterogeneous media, with an application to the kilauea volcano[END_REF]. Sources were placed at different depths below the crater and receivers are situated at the same positions as the two experimental cross-shaped arrays. A high-resolution method based on the multiple signal classification (MUSIC technique) [START_REF] Bienvenu | Optimality of high resolution array processing using eigensystem approach[END_REF], [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], but adapted to the 3C case, is applied to the synthetic data to determine the back-azimuth, the apparent velocity and the incidence angle for both arrays and all sources. This procedure is then applied to the real data recorded in March 2009.

UBINAS VOLCANO

Ubinas volcano (16 22' S, 70 54' W; altitude 5672 m) began to erupt on March 25 th 2006 after nearly 40 years of quiescence. Situated in the Central Volcanic Zone (CVZ, southern Peru), Ubinas volcano is an active andesitic stratovolcano [START_REF] De Silva | Volcanoes of the central andes[END_REF] truncated in the upper part by a caldera 600 m in diameter. The caldera floor is a flat area lying approximately 5100 m above sea level. The active crater is situated in the southern section; the bottom is 300 m under the caldera floor (Figure 1a). Ubinas is considered the most active Peruvian volcano during the last 500 years, threatening 3,500 people living on the edge of the Ubinas rio [START_REF] Rivera | Ubinas, el volcan mas activo del sur del peru desde 1550: Geologia y evaluacion de las amenazas volcanicas[END_REF]. Arequipa airport, situated 60 km east of the volcano, has been closed several times since 2006 due to ash emissions. The Instituto Geofisico del Peru (IGP) with the cooperation of the Institut de Recherche pour le Developement (IRD) has started seismic monitoring of the volcano to understand the activity associated with this eruptive sequence under the 6 th framework EU project VOLUME (http://www.volume-project.net/). A network of four digital 1 Hz stations with a radio telemetry system has been operating since 2006. Data are transmitted to Arequipa observatory. At the time of the experiment, the eruption was characterized by an almost permanent ash emission. Two main types of degasification were observed: 1) exhalations rising a few hundred meters above the crater rim and 2) plumes produced by explosions that may reach 10 km, critical for aircraft duty. This activity is thought to be related to a magmatic plug positioned at the bottom of the south part of the caldera wall [START_REF] Macedo | Seismicity associated with the 2006-2008 eruption ubinas volcano[END_REF].

PROCESSING APPROACH

Three component broadband seismic array data is suitable for wideband multidimensional signal processing techniques. The antenna is designed (size, shape, aperture) to avoid spatial and temporal aliasing [START_REF] Mars | Advanced signal processing tools for dispersive waves[END_REF]. To illustrate the 3C MUSIC method, we assume an isotropic source radiation pattern and a homogeneous medium. Assuming K sources impinging at the antenna of N triaxial sensors, equation (1) represents the received signal by n th sensor in time-space (t -x) domain.

xn(t) = K X k=1 s k (t -τ n,k ) + bn(t); n = 0, 1, ..., N -1 (1)
where, s k (t) is the k th source signal received in the first sensor, the τ n,k is the relative propagation time delay of k th source for the n th sensor and the bn(t) represents the noise in n th sensor (we assume white noise and Gaussian with variance σ 2 and uncorrelated with the sources). The corresponding relative time delay is defined as

τ n,k = 1 v 0 [rn • u(θ k , φ k )]
where, v0 is the wave propagation apparent velocity beneath the array, the rn is the relative location for sensor n th respect to the first sensor located at (0,0,0) and the u(θ k , φ k ) is the slowness vector of the impinging wavefield from k th source. The dot product (rn • u) represents the propagation path between the n th sensor respect to the first sensor. The propagation direction u(θ k , φ k ) of the k th source arrived in the antenna can be described by back-azimuth angle θ k and vertical incidence angle φ k , these two parameters are combined in the slowness vector as is shown in equation (2).

u(θ k , φ k ) = 2 4 -sinθ k sinφ k -cosθ k sinφ k cosφ k 3 5 (2)
Taking into account the wavenumber k, that is related as k = 1 λw u = fo vo u, where, λw is the wavelength, fo the wave frequency and vo is the wave propagation apparent velocity. The Fourier transform in equation ( 1) can be represented in the frequencywavenumber (f -k) domain as:

Xn(fo) = K X k=1 S k (fo)exp(-2πjrn • k) + Bn(fo) (3)
where, j = √ -1 is the imaginary unit, and S k (f ) and Bn(f ) are the Fourier transform of s k (t) and bn(t) respectively. Then the antenna output X(fo) for narrowband signals can be represented as a Nx1 matrix form as:

X(fo) = A(fo) • S(fo) + B(fo) (4)
where, the NxK matrix A(fo) is called "array reponse" or "steering matrix", the Kx1 matrix S(fo) keeps the sources signals defined in (5), so the Nx1 matrix B(fo) represents the noise.

A(fo) = [a(θ1, φ1), a(θ2, φ2), ..., a(θK , φK)] (5) a (θ, φ) = 2 6 6 4 1 exp(-j2πr1 • k) ....... exp(-j2πrN-1 • k) 3 7 7 5 S(fo) = [s1(fo), ..., sK (fo)] T
Note that the array manifold a (θ, φ) depends on the frequency fo of the signal through the wavenumber k. The purpose of this paper is to develop and test a 3C MUSIC algorithm in order obtain a more robust estimation of the back-azimuth angle θ and a reliable incidence angle φ related to source depth determination.

3C MUSIC algorithm

MUltiple SIgnal Classification is an eigenstructure subspace analysis method, that is widely used in geophysics, particularly to enhance the signal-to-noise ratio. MUSIC algorithm tries to find the array manifold orthogonal to the noise subspace since the noise subspace is orthogonal to the signal subspace [START_REF] Mars | Advanced signal processing tools for dispersive waves[END_REF]. A time window array data is selected, correspondig to either explosion quakes or LP events. A way to find out the dominants frequencies in the wavefield is through average energy spectrum of each component, which the analysis frequency fo is determined. Observing M frequencies snapshots around the bin fo (energy peak) are used to estimate the cross-spectral matrix. The i th snapshot data is shown in equation ( 6) defined by a NxM matrix [START_REF] Paulus | New multicomponent filters for geophysical data processing[END_REF].

Xi = 2 6 6 4 W Ei,1 SNi,1 Zi,1 . . . . . . W Ei,N SNi,N Zi,N 3 7 7 5 ; i = 1, ..., M (6)
where, the triaxial sensor is represented by: the west-east component as W E, south-north as SN and vertical as Z. The crossspectral matrix b Γ can be estimated using the equation ( 7) [START_REF] Paulus | Wideband spectral matrix filtering for multicomponent sensors array[END_REF] and .

b Γx = 1 M M X i=1 Xi • X H i (7)
where H denotes conjugate transpose. The cross-spectral matrix in (7) can be expressed in terms of its N eigenvalues and eigenvectors (N=sensors number), as shown the equations ( 8).

b

Γx (ν) = N X n=1 λnvn • v H n (8)
where, λn and vn are the n th eigenvalue and eigenvector respectively. Anyway the sample cross-spectral matrix is the expected value of the outer product of the data vector. Using the orthogonality property between signal and noise subspaces, in the equation ( 4), the cross-spectral matrix Γx can be written as:

Γx = A • E{S • S H } • A H + E{B • B H } (9)
Γx (fo) = A (fo) Γs (fo) A H (fo) + σ 2 B (fo) I where, E{•} the expected operator, Γs (fo) is the signal crossspectral matrix, I is the identity matrix and σ 2 B is the variance of the Gaussian white noise. This decomposition can be split into two orthogonal subspaces such as signal and noise. In this case the cross-spectral matrix ( 9) is equal to the estimated cross-spectral matrix in ( 7) and ( 8). Due to errors in the estimated cross-spectral matrix, the eigenvalues of the noise subspace are no longer equal to σ 2 B . Assuming K sources impinging the antenna of N sensors (K < N ), the eigenvectors whose eigenvalues are the K largest eigenvalues define the signal subspace and the remaining N-K eigenvectors define the noise subspace. The largest eigenvalues are easily distinguishable from the rest in the eigenvalues λi, as λ1 ≥ ... ≥ λK ≥ λK+1 > ... > λN ≈ σ 2 shows the appropriate dimension of noise subspace on the N -K smallest. Then the noise projector is defined by [START_REF] Miron | Vector-sensor music for polarized seismic sources localization[END_REF]:

Π ⊥ = N X k=K+1 v k v H k
and the MUSIC estimator is [START_REF] Wong | Self-initiating music direction finding & polarization estimation in spatio-polarizational beamspace[END_REF]:

3CMUSIC (θ, φ) = 1 A H (θ, φ) Π ⊥ A (θ, φ) (10)
Then 3C MUSIC algorithm follows the procedure: 1) Selection by time windowing data for either a quake explosion or LP event. 2) Average power spectra density to determine the frequencies bins and entimate the cross-spectral matrix. 3) Eigen structure analysis to estimate the noise sub-space to build the noise projector. 4) Evaluate the 3C MUSIC estimator on slowness grid and velocity grid where the estimator is maximum.

DATA

A seismic experiment was carried out at Ubinas volcano between May and July 2009. We deployed two cross-shaped antennas with triaxial-broadband seismometers as depicted in Figure 1a. NUBI antenna was installed on a flat area 4632 m above sea level. It consists of twelve Guralp CMG-6TD seismometers (see reference for Guralp manual). WUBI antenna was installed on a sloping area from 4752 m to 4883 m above sea level. It consists of six Guralp CMG-6TD and six GEOMAX seismometers (see reference for Geomax manual) both with 24 bit seismic recorders (dynamic range 120 dB) and high accuracy time synchronization, GPS-based.

The acquisition performed continuous recording with a sampling frequency of a 100 Hz. The horizontal components of such seismometers were oriented to geographical coordinates West-East (E component) and South-North (N component) respectively. The distances between seismometers was set to approximately 50 meters. Each seismometer location was surveyed by Trimble GPS (GeoXH, accuracy 0.3 m). Twelve explosion quakes and hundreds of LP events were recorded during the experiment. A synthetic dataset was generated to test the accuracy of the location methods. The synthetic dataset was created using a digital elevation map of Ubinas topography and the 2009 experimental array locations. A 3D discrete numerical elastic lattice method (O'Brien & Bean, 2004) was used to propagate waves in the structure with six broadband isotropic sources (Figure 1b) located beneath the summit. The velocity structure was generated from a tomography study [START_REF] Monteiller | An efficient algorithm for double-difference tomography and location in heterogeneous media, with an application to the kilauea volcano[END_REF] assuming a density of 2300 kg/m 3 . The sampling frequency was fixed to 250 Hz.

Synthetic data analysis

To investigate the performance of the 3C MUSIC technique, we applied it to our synthetic dataset. For source number 3 (Figure 1b) and the synthetic array NUBI a one second time window including the first arrival was selected, Figure 2a. In order to find the dominant frequency, we averaged the power spectral density computed for the selected time window of each component and for each receiver. The dominant energy peak is located at 2.23 Hz (Figure 2b). The cross-spectral matrix was then calculated by using 32 windows around the bin. The eigendecomposition of the cross-spectral matrix gave us an overview of the signal-to-noise subspaces, by sorting the eigenvalues. In the first stage, we determined the eigenvectors corresponding to the noise subspace, which gives the projector.

In the second stage, we calculated the steering vector which consists of the slowness vector estimated between 0 and 2π with steps of 3.6 • . Given the projector and the steering vector, we obtained the 3C MUSIC spectrum represented in Figure 3a. The maximum amplitude of the spectrum gives the back-azimuth and the velocities. The largest estimator point is located at 181 ±3 • for the backazimuth and 2900±75 m/s for the velocity. The vertical and horizontal cross-section of the 3C MUSIC spectrum are shown in Figures 3b and 3c. The incidence angle is deduced from the estimated values, back-azimuth and velocity. We found 85.5±6 • as shown in Figure 3d. The errors are estimated by taking the peak width of 95% of the maximum amplitude. These values agree well with the theoretical values.

We applied the same analysis to the synthetic data generated for both arrays and for all six sources. In order to compare the 3C and 1C methods, we also analyzed the same synthetic dataset using only the vertical components. The results are shown in the Figure 4a and Figure 4b for NUBI and WUBI antennas respectively. Back-azimuth results obtained both with the 3C MUSIC and the 1C MUSIC fit well with the theoretical values for both antennas. Errors are ±3 • and ±6 • , with the 3C and the 1C MUSIC algorithm respectively. However, the incidence angles are significantly different when using the 3C MUSIC and the 1C MUSIC analysis. For NUBI, values vary from 75 • for the deepest source to 93 • for source 1. These incident angle variations are coherent with the variation of the sources depth. Errors are ±6 • for the 3C MUSIC analysis. These errors allow us to distinguish sources depths separated by several hundreds of meters but unable to separate close sources such as 1 and 2 sources 4 and 5. The 1C MUSIC analysis gives sources close to the theoretical values for the superficial sources (1, 2 and 3) while depths for sources 4, 5 and 6 are far away from the theoretical values. The incident angle remains approximately the same value around 90 • . The error obtained with the 1C MUSIC processing is ±12 • . None of the sources can be distinguished with the 1C MUSIC analysis. For WUBI antenna, the 3C MUSIC analysis gives incident angles between 63 • for source 6 and 100 • for source 2 (Figure 4b). As for NUBI, solutions follow the depth variation. Error are approximately 6 • . The 1C MUSIC analysis gives higher incident angles than the theoretical values (90 • for source 6 to 156 • for source 1). Errors are larger than the 3C algorithm at approximately ±15 • . The 1C MUSIC analysis does not give reliable solutions for any of the synthetic source depths. In this case, it seems that the solutions are affected by the topography below the antenna. Similar results are observed for the velocity. In summary, incident angles obtained by the 3C MUSIC algorithm are close to the theoretical values whereas those obtained with 1C MUSIC are not reliable.

Real Data Analysis

In this section, we discuss the performance of 3C MUSIC method by applying it to an explosion quake signal recorded at Ubinas volcano. The selected event, recorded at 01:27 PM, June 2 nd , 2009 is shown in Figure 2c (SHOW A LONGER TIME WINDOW OF THE SIGNAL). The data was corrected for the instrument response and bandpass filtered between 1 Hz -10 Hz. Figure 2d shows the antenna energy spectrum average obtained from two seconds of signal including the first arrival wave. The dominant peak is centred at 2.39 Hz. The results for NUBI antenna are presented in Figure 3e, 3f, 3g and 3h). (IF YOU WANT TO SEND TO DIFFERENT JOURNAL I'D SEPARATE THE SYNTHETIC AND DATA FIG-URES) The maximum of the 3C MUSIC spectrum gives a backazimuth of 184±5 • , a velocity of 2975±125 m/s and an incidence angle of 85.9±7 • . For WUBI, we obtained a back-azimuth of 119±6 • , a velocity of 3100±120 m/s and an incidence angle of 75±7 • . The back-azimuth and the incident angle can be repre-sented as gaussian variables with mean θ (0≤ θ ≤ 360 • ) and φ (0≤ φ ≤ 120 • ) and standard deviations σ θ and σ φ corresponding to the values found in the analysis. This allows the definition of a probability density function (PDF) ρ(θ k ) of the back-azimuth and a PDF ρ(φ k ) of the incident angle for each array k. ρ1(θ k ) and ρ1(φ k ) are shown in Figure 5 as rose diagrams in the horizontal plane for the back-azimuths and in the vertical planes oriented N-S and W-E for the incident angle. The last step is to locate the source by using the information obtained at each antenna. For each point with geographical coordinates (x,y,z) in the source region and each array k, the back-azimuths θ k (x, y, z) and the corresponding value of ρ1(θ k ), as well as the incident angle φ k (x, y, z) and the corresponding value of ρ1(φ k ) can be calculated. The PDF of the source position is derived from the different PDF's of the back-azimuth and the incident angle:

ρ2(x, y, z) = 2 Y k=1 ρ1(θ k (x, y, z)) • ρ1(φ k (x, y, z)) (11)
The maximum likelihood of the PDF ρ2 yields an estimate of the source location. We define the mean quadratic radius R = p (σ 2 1 + σ 2 2 + σ 2 3 )/3, where σ 2 1 , σ 2 2 and σ 2 3 are the principal variances of ρ2(x, y, z). Figure 5 shows the result obtained with the explosion quake analysis. The source area is situated 150 m West and 1000 m below the bottom of the crater at an altitude of 4200 m, the radius R is 660 m.

CONCLUSIONS

We have presented a source localization method, 3C MUSIC, based on the use of 3 component arrays and, for comparison, with the 1C MUSIC used in previous studies on volcanoes [START_REF] Saccoroti | Performances of the underground seismic array for analysis of seismicity in central italy[END_REF]). Synthetics have been generated for 6 sources with elevations from 3400 m to 4950 m a.s.l. The back-azimuths from the 3C MUSIC correspond to the theoretical values for both antennas with a resolution of ±3 • . 1C MUSIC gives equivalent results, but with higher errors (±6 • ). The incident angle varies with depth when it is determined with 3C MUSIC. The incident angle is determined with an error of ±6 • for NUBI and for WUBI. Knowing the distance from the center of the antennas and the hypocenter of the sources, depth resolution can be deduced for each antenna. It is 500 m for NUBI and 400 m and WUBI. On the other hand, the 1C MUSIC analysis does not allow the depth to be determined. It seems that the topography affects the results obtained with the 1C MUSIC. The velocity follows the same behaviour as the incident angle. The theoretical velocity has an accuracy of ±150 m/s using the 3C MUSIC algorithm. The 1C MUSIC measures higher velocities at NUBI and lower velocities at WUBI. Finally, we located an explosion quake recorded during the field experiment using the 3C MUSIC. This signal is characteristic of the explosive activity observed at Ubinas volcano [START_REF] Macedo | Seismicity associated with the 2006-2008 eruption ubinas volcano[END_REF]). We found a source located at 4200±500 m a.s.l. It is situated more than 1000 m below the summit of the intrusive conduit. We conclude that 3C MUSIC provides realistic values of the depth of volcanic sources, unlike the 1C MU-SIC or other antenna methods based on time delays measurements. Given the performance of the 3C MUSIC algorithm, we will apply it to other explosions and other types of volcanic signals (LP, tremors) recorded at Ubinas to better characterize the eruptive dynamics of this volcano. In addition, the 3C MUSIC will be tested with the IGP monitoring system to try to locate the seismic activity in real time. This algorithm is not restricted to volcanic sources but can be used to locate other types of non-volcanic signals.

Figure 1 .

 1 Figure 1. a) Map of Ubinas volcano showing the location of the north (NUBI) and the north-west (WUBI) antennas, b) Horizontal profile showing the position of the synthetic sources Figure 2. a) The vertical (Z), south-north (N) and west-east (E) synthetic seismograms calculated at the central receiver of NUBI antenna, b) Averaged energy spectrum calculated for all the receivers and all the components, c) The vertical (Z), south-north (N) and west-east (E) explosionquake event waveforms (real data) recorded by the central receiver of NUBI antenna, d) Averaged energy spectrum calculated for all the receivers and all the components. The vertical dash lines in a) and c) indicate the time window (P-wave) selected for the processing. The vertical dash lines in c) and d) represent the frequency windows used for the cross spectral matrix calculation.

Figure 3 .

 3 Figure 3. Results obtained for the synthetic data calculated for the six sources. Open triangles and open stars represent results obtained with the 3C MUSIC and 1C MUSIC, respectively. Results obtained for an explosion quake (real data) using the 3C MUSIC is represented by a closed circle. The abscissa represents the altitude. Sources are numbered as in figure 1b. a) Back-azimuth, incidence angle and velocity for NUBI antenna. b) Back-azimuth, incidence angle and velocity for WUBI antenna. Dash lines represent the theoretical values.

Figure 4 .

 4 Figure 4. a) Normalized 3C MUSIC sprectrum calculated with syntetic data generated at source 3 for NUBI antenna. b) Normalized back-azimuth profile (cross section at velocity 2900 m/s). c) Normalized velocity profile (cross section at back-azimuth at 181 degrees). d) Normalized Incidence angle. e) Normalized 3C MUSIC spectrum for explosion-quake event. f) Normalized back-azimuth profile (cross section at velocity 2925 m/s). g) Normalized velocity profile (cross section at back-azimuth at 184 degrees). h) Normalized Incidence angle. The vertical dot lines represent the error range.

Figure 5 .

 5 Figure 5. Probability density function of the source position ρ 2 . a) horizonthal view at 4200 m depth, b) vertical view oriented West-East crossing the maximum likelihood of ρ 2 , c) same as b) oriented North-South. The PDF ρ 1 (θ k ) and ρ 1 (φ k ) are represented as rose diagrams with an increment of 5 •