MSc Julie Samyde

MD Pierre Petit

MD Dominique Hillaire-Buys

MD Jean-Luc Faillie
email: jl-faillie@chu-montpellier.fr

Quinolone antibiotics and suicidal behavior: analysis of the World Health Organization's adverse drug reactions database and discussion of potential mechanisms

Keywords: psychopharmacology, adverse drug reactions, suicidal behaviors, suicide, quinolones

Rationale Several case-reports suggest that the use of quinolones may increase the risk of psychiatric adverse reactions such as suicidal behaviors.

Objectives The aim of this study is to investigate whether there is a safety signal for quinolone-related suicidal behaviors in a global adverse drug reactions database. Methods All antibiotic-related adverse reactions were extracted from VigiBase, the World Health Organization (WHO) global Individual Case Safety Report (ICSR) database. Disproportionality analyses were performed to investigate the association between reports of suicidal behavior and exposure to quinolones, in comparison with other antibiotics. Results From December 1970 through January 2015, we identified 992 097 antibiotic-related adverse reactions. Among them, 608 were quinolone-related suicidal behaviors including 97 cases of completed suicides. There was increased reporting of suicidal behavior (adjusted reporting odds ratios [ROR] 2.78, 95 % CI 2.51-3.08) with quinolones as compared to other antibiotics. Candidate mechanisms for quinolone-induced suicidal behaviors include GABAA antagonism, activation of NMDA receptors, decreased serotonin levels, oxidative stress and altered microRNAs expression. Conclusions We found a strong safety signal suggesting an increased risk of suicidal behaviors associated with quinolones use. Plausible psychopharmacological mechanisms could underlie this association. Further investigations are urgent to confirm and better understand these findings.

Introduction

Quinolones antibiotics are frequently used for the treatment of a variety of infectious diseases. This drug class is popular because of its relatively broad spectrum of activity, multiple approved indications for use, and high bioavailability, resulting in comparable blood concentrations when given intravenously or orally (Melhorn and Brown, 2007). The mechanism of action of quinolones involves the inhibition of the enzymes topoisomerase type II (also known as DNA gyrase) and topoisomerase type IV. These enzymes are involved in bacterial DNA synthesis, and are essential for bacterial DNA replication. Thus, quinolones are both specific and bactericidal antibiotics. More precisely, DNA gyrase tends to be the primary target of quinolones in Gram-negative organisms, whereas topoisomerase IV is typically the primary target in Gram-positive bacteria but these two mechanisms can be combined (Zhanel and Noreddin, 200;[START_REF] Blondeau | Fluoroquinolones: mechanism of action., classification., and development of resistance[END_REF]. Well established adverse reactions related to quinolones include allergies, phototoxicity, tendinopathies, cardiac arrhythmias or peripheral neuropathy. In the recent past, several quinolones were withdrawn from the market because of the occurrence of serious adverse drug reactions: temafloxacin for haemolytic anemia, trovafloxacin for hepatotoxicity, grepafloxacin for QTc interval prolongation, clinafloxacin for phototoxicity or gatifloxacin for dysglycemia [START_REF] Tomé | Quinolones: review of psychiatric and neurological adverse reactions[END_REF][START_REF] Lewis | Dysglycaemias and fluoroquinolones[END_REF][START_REF] Dukes | Fluoroquinolones. In Meyler's side effects of drugs: the international encyclopedia of adverse drug reactions and interactions[END_REF][START_REF] Sprandel | Safety and tolerability of fluoroquinolones[END_REF][START_REF] Smith | Drugs withdrawals and the lessons within[END_REF][START_REF] Rubinstein | History of quinolones and their side effects[END_REF]. Other quinolones such as moxifloxacin or norfloxacin are only available in certain countries or are used with restrictions (European Medicines Agency, 2008).

Several case reports have raised concerns that exposure to quinolones antibiotics may increase the risk of suicidal behavior i.e. suicidal ideations, suicide attempt or completed suicide. In a review of the literature performed in 2011, Tomé and Filipe identified 111 case reports describing psychiatric adverse drug reactions (ADRs) associated with the use of quinolones, including two depressions, two manic depressions and one report of suicidal ideations. Other published case reports described violent suicidal behaviors like stabbing and cutting with a kitchen knife [START_REF] Behera | Antitubercular drug-induced violent suicide of a hospitalised patient[END_REF], defenestration [START_REF] Tomé | Quinolones: review of psychiatric and neurological adverse reactions[END_REF], or cutting the both wrists and the neck before electrocution (Barthélemi and Lebrun, 2015) in patients exposed to quinolones. One resulted in the death of the patient [START_REF] Behera | Antitubercular drug-induced violent suicide of a hospitalised patient[END_REF]. A summary of published case reports is provided in Table S1. Recently, a review of the adverse drug reactions spontaneously reported to the French pharmacovigilance system also suggested a possible association between completed suicide and the use of quinolones. Among the reported cases, three completed suicide by defenestration, two to five days after initiation of quinolones. Patients were men aged from 62 to 82 years and two of them had a history of depression (Barthélemi and Lebrun, 2015). In summary, the case reports suggest violent suicidal behaviors, a predominance of male gender, a relatively short time to onset after exposure to quinolones and suspected drugs which include ofloxacin, levofloxacin and ciprofloxacin [START_REF] Behera | Antitubercular drug-induced violent suicide of a hospitalised patient[END_REF][START_REF] Labay-Kamara | Fluoroquinolone-induced suicidal ideation and suicidality[END_REF][START_REF] Barthélémi | Completed suicide and fluoroquinolones: a possible association[END_REF][START_REF] Ahmed | A man who wanted to commit suicide by hanging himself: an adverse effect of ciprofloxacin[END_REF][START_REF] Lasalvia | Fluoroquinolone-induced suicidal ideation[END_REF]. According to the European Medicine Agency, a safety signal is generated from several sources such as spontaneous reports, clinical studies and the scientific literature and informs on an adverse event that is potentially caused by a medicine and that warrants further investigation.

In order to investigate a safety signal for quinolone-related suicidal behaviors in a large adverse event database, we compared the reporting of adverse events of suicidal behaviors among quinolone users with users of other antibiotics by performing a disproportionality analysis within VigiBase, the World Health Organization (WHO) global individual case safety report (ICSR) database [START_REF] Lindquist | VigiBase, the WHO global ICSR database system: basic facts[END_REF]. We also discussed the potential mechanisms involved in quinolone-related suicidal behaviors.

Methods

Data source

The study was based on information extracted from VigiBase, the largest and most comprehensive pharmacovigilance database which collects ICSR from 110 countries over the five continents. The National Pharmacovigilance Centres (NPCs), who are members of the WHO Programme for International Drug Monitoring [START_REF] Lindquist | The WHO programme for international drug monitoring, its database., and the technical support of the Uppsala monitoring center[END_REF], submit the data on their spontaneous reports of ADRs to VigiBase, which is maintained by the Uppsala Monitoring Centre (UMC) in Sweden. The data, which origins from a variety of sources (health care professionals, pharmaceutical companies and patients) is evaluated by the relevant NPCs prior to submission and registration in VigiBase. Over 10 million ICSRs are currently recorded in VigiBase, and about 1 million ICSRs are added each year. The reports usually contain reporter information, patient characteristics (age, gender, and medical history), clinical description of the ADR, its seriousness and its evolution, drug exposures with dates, dosages and indications. In VigiBase, the WHO Adverse Reaction Terminology (WHO-ART) or the Medical Dictionary for Regulatory Activities (MedDRA) [START_REF] Brown | The medical dictionary for regulatory activities (MedDRA)[END_REF] are used to code ADRs (WHO-ART coded reports are assigned a corresponding MedDRA code on the Preferred Term level). Drugs are coded according to the WHO Drug Dictionary Enhanced and classified using the Anatomical Therapeutic Chemical (ATC) classification (WHO, 2013). Details about each suspected adverse reaction received by the UMC may vary from one country to another. Thus, the information in these reports is not homogenous, at least with regard to origin, completeness of documentation, or likelihood that the suspected drug caused the adverse events.

In this study, we extracted from VigiBase all the reports of adverse drug reactions (identified with MedDRA version 18.1) associated with an antibiotic exposure (ATC code J01) recorded from December 9 th , 1970 to January 19 th , 2015. This selection includes psychiatric and non-psychiatric events as well as serious and nonserious events. In our study, several ICSRs may refer to the same patient.

Study design

We studied the disproportionality of reports for suicidal behaviors associated with quinolones by performing a case-non case study. This design allows the calculation of Reporting Odds Ratios (ROR), which compares the odds of exposure to quinolones between cases and non-cases [START_REF] Moore | The history of disproportionality measures, reporting odds ratio, proportional reporting rates in spontaneous reporting of adverse drug reactions[END_REF][START_REF] Van Puijenbroek | A comparison of measures of disproportionality for signal detection in spontaneous reporting systems for adverse drug reactions[END_REF]. Cases of suicidal behaviors were reports of an ADR defined by the MedDRA High Level Term (HLT) "Suicidal and self-injurious behaviors". All other reports associated with an antibiotic exposure were considered as noncases.

Exposure definition Exposure to quinolones drugs was identified in the ICSR by the use of at least one quinolone (ATC code J01M) preceding the onset of the adverse reaction. The most frequently reported quinolones (ciprofloxacin, levofloxacin, moxifloxacin and ofloxacin) were also identified and analyzed separately. Comedications with antidepressant or antipsychotic drugs were retrieved by the ATC codes "N06A" and "N05A", respectively.

Statistical analysis

Descriptive statistics were used to summarize the baseline characteristics retrieved from the case reports of suicidal behaviors and as well as completed suicides. A temporal year-by-year presentation of the reported cases was also provided.

Among all adverse drug reactions reports associated with antibiotic use, our primary analysis consisted in calculating the ROR of suicidal behavior (and corresponding 95% confidence interval) for quinolones compared to other antibiotic drugs using logistic regression models adjusted for age and gender. The ROR for the frequently reported quinolones (ciprofloxacin, levofloxacin, moxifloxacin, ofloxacin) were also calculated.

Our first secondary analyses studied the variation in reporting rate among different regions and periods: We performed a geographical analysis among the various regions and countries, and the presence of a notoriety effect [START_REF] Pariente | Impact of safety alerts on measures of disproportionality in spontaneous reporting databases: the notoriety bias[END_REF] was studied by a year-by-year analysis of the cumulative crude RORs during the more recent years (from 2000 to 2015). Second, in order to study the influence on ROR estimations of a potential event-competition bias by well-established ADRs associated with quinolones [START_REF] Salvo | A potential event-competition bias in safety signal detection: results from a spontaneous reporting research database in France[END_REF], we repeated the disproportionality analyses when excluding well-known quinolone-related adverse reactions such as musculoskeletal and dermatological events. Third, we used the same design than the primary analysis to separately study the disproportionality of the reporting of depressive disorders on one hand and completed suicides on the other hand, in association with quinolone use. Cases of depressive disorders were defined as reports of ADRs corresponding to the MedDRA High Level Group Term (HLGT) "Depressed mood disorders and disturbances" and cases of completed suicide were retrieved by selecting the reports of ADRs including the MedDRA Preferred Term (PT) "Completed suicide". Finally, to explore a potential effect modification by age or gender, we performed stratified analyses according to age categories and gender. Analyses were conducted using SAS software version 9.2 (SAS Institute, Cary, NC, USA).

Identification of potential mechanisms

A separate literature survey used the Medline database to search for data on possible underlying mechanisms of psychiatric disorders induced by the use of quinolones.

Results

Descriptive statistics

We identified 992 097 ICSRs associated with antibiotic use reported to VigiBase from 117 countries between December 1970 and January 2015. Reported cases of suicidal behavior were 1627 (0.2%). Among the ADRs with exposure to quinolones, there were 608 cases of suicidal behavior including 97 cases of completed suicide and 171 927 other adverse reactions (Figure 1). The characteristics of reported cases of suicidal behavior and completed suicide that involved quinolones are presented in Table S2. The majority of cases of suicidal behavior exposed to quinolones were reported after 2000 (93.4%) and originated from America and Europe (61.3% and 36.5% respectively). The two countries which reported the most were the United States (57.1%) and Germany (23.5%). Mean age of exposed cases of suicidal behavior was 40.4 and 53% were women. The majority of the cases of completed suicide exposed to quinolones occurred in subjects aged from 45 to 64 (50.5%) and mean age was 44.1. Exposed cases of suicidal behavior were most frequently reported by health care professionals (46.4%). Ciprofloxacin was the most frequently used quinolone in the reports of suicidal behavior (43.2% of the cases), followed by levofloxacin (26.6%), moxifloxacin (14.9%) and ofloxacin (11.4%). In total, 13.7% of the quinolone-related cases of suicidal behavior were also exposed to antidepressants and 5.8% to antipsychotics. Comedications with antidepressants and antipsychotics were more frequent for the quinolonerelated cases of completed suicide (37.2% and 26.7% respectively).

Annual counts of reported suicidal behavior and completed suicide exposed to quinolones are presented in Figure 2. The gap in reports for the year 2007 and the following peak in 2008 is an artifact due to the implementation by the UMC in 2007 of a new method to enter ICSRs into VigiBase resulting in delayed import of the ICSRs of 2007 (Uppsala Monitoring Centre, 2008). Nonetheless, the annual number of quinolone-related cases recorded in VigiBase is clearly higher after 2007, for both suicidal behaviors and completed suicides.

Primary analysis

Table 1 presents the disproportionality analysis of reported suicidal behaviors associated with quinolones. Overall, the use of quinolones was associated with a significantly increased reporting of suicidal behavior compared with other antibiotics (adjusted reporting odds ratio [ROR] 2.78, 95% confidence interval [CI] 2.51-3.08). Among the most frequently reported drugs, ciprofloxacin was the most strongly associated with increased reporting of suicidal behaviors (adjusted ROR 4.01, 95% CI 3.50-4.59). The ROR for the other quinolones levofloxacin, moxifloxacin and ofloxacin also showed significant increase in reporting, with greater values than the one observed for the whole class (Table 1).

Secondary analyses

The use of quinolones was associated with a statistically significant increase in reporting of suicidal behavior compared with other antibiotics in all the regions and countries where cases were notified, except in Germany where the safety signal was not statistically significant (ROR 0.84, 95% CI 0.69-1.03). Otherwise, among the countries which reported more than 10 cases, ROR ranged from 1.91 (in United Kingdom) to 5.32 (in Canada). The detailed geographic analysis of reported suicidal behaviors associated with quinolones in the main regions and countries is provided in Table S3.

The year-by-year analysis of cumulative ROR of suicidal behaviors associated with quinolones is presented in Figure 3. The cumulative crude ROR in the period from 1970 to 2000 was relatively high but calculated in a small proportion of the total cases exposed to quinolones (40 cases, 6.6%). As of 2000, ROR progressively decreased and, since 2010, stabilized at a value around 3.

When excluding well-known quinolone-related adverse reactions such as dermatological events, musculoskeletal events and both, ROR for quinolone related suicidal behaviors were slightly increased (values ranging between 3.01 and 3.83) but showed a similar significant disproportionality signal than the primary analysis (Table S4).

Disproportionality analyses of reported depressive disorders and completed suicide associated with quinolones are shown in Tables S5 andS6. Compared to other antibiotics, exposure to quinolones was associated with an increased reporting of depressive disorders (adjusted ROR 4.15, 95% CI 3.87-4.45) and completed suicide (adjusted ROR 1.56, 95% CI 1.24-1.97). Among the most frequently reported drugs, ofloxacin was the most strongly associated with increased reporting of depressive disorders (adjusted ROR 5.99, 95% CI 5.20-6.89) (Table S5) and ciprofloxacin was the most strongly associated with increased reporting of completed suicide (adjusted ROR 3.49, 95% CI 2.67-4.55) (Table S6).

Stratified analyses according to age categories and gender are presented in Tables S7 andS8. Slightly increased ROR of suicidal behaviors were observed for the categories <17 years and 65 years -74 years (adjusted ROR 4.10, 95% CI 2.42-6.91 and 5.85, 95%CI 3.64-9.41, respectively). ROR for completed suicide were clearly increased for the younger category <17 years (adjusted ROR 8.96,. Stratification according to gender did not provide relevant differences compared to the primary analysis.

Discussion

Our study used a case non-case design to evaluate the association between the use of quinolones and the reporting of suicidal behaviors as well as depressive disorders and completed suicides, recorded in the WHO drug adverse event database from 1970 to 2015. Strong significant disproportionality signals were found, showing an increased reporting of suicidal behaviors, completed suicides and depressive disorders with quinolone drugs and particularly with ciprofloxacin, levofloxacin, moxifloxacin and ofloxacin, the most frequently reported drugs. Although stratified analyses were based on fewer exposed cases, the increased ROR observed for completed suicide in the category <17 years suggests that young patients may be more at risk when using quinolones. Although quinolones are usually not recommended in children, this might be related to a greater vulnerability of younger patients to psychiatric harms of drugs, as it was suggested for suicidality among young patients exposed to antidepressants [START_REF] Stone | Risk of suicidality in clinical trials of antidepressants in adults: analysis of proprietary data submitted to US Food and Drug Administration[END_REF].

The potential risk of suicidal behaviors and more generally for neuropsychiatric adverse effects associated with quinolones has been suggested by several case reports [START_REF] Tomé | Quinolones: review of psychiatric and neurological adverse reactions[END_REF]Behera et al, 2013;[START_REF] Labay-Kamara | Fluoroquinolone-induced suicidal ideation and suicidality[END_REF][START_REF] Barthélémi | Completed suicide and fluoroquinolones: a possible association[END_REF][START_REF] Ahmed | A man who wanted to commit suicide by hanging himself: an adverse effect of ciprofloxacin[END_REF][START_REF] Lasalvia | Fluoroquinolone-induced suicidal ideation[END_REF][START_REF] Montané | Multiple drug interactions -induced serotonin syndrome: a case report[END_REF][START_REF] Mulhall | Ciprofloxacin induced acute psychosis[END_REF][START_REF] Feinberg | Fluoroquonolones-induced depression[END_REF]. Quinolones were associated with various neuropsychiatric adverse drug reactions, including delirium, agitation, confusion, psychosis, seizures/convulsions, anxiety and major depression (Melhorn and Brown, 2007;[START_REF] Dukes | Fluoroquinolones. In Meyler's side effects of drugs: the international encyclopedia of adverse drug reactions and interactions[END_REF][START_REF] Feinberg | Fluoroquonolones-induced depression[END_REF][START_REF] Fish | Fluoroquinolone adverse effects and drug interactions[END_REF]. Suicidal behavior is a heterogeneous and multifactorial entity involving complex interactions between biological, psychological and socio-environmental factors. Even though suicide is now considered a condition on its own and not just a consequence of other psychiatric disorders, the presence of a psychiatric disorder, particularly depression and psychosis, represents a main risk factor for suicide.

Potential mechanisms

A number of neurobiological actions that could explain the psychiatric adverse effects of quinolones have been reported. Since the neurobiology of suicide may be, at least in part, independent of other psychiatric disorders, the question whether quinolone-induced suicidal behavior is independent or secondary to the psychopathological context remains to be clarified.

Impairment in the GABAergic/glutamatergic balance in the brain has been suggested as a factor contributing to the neurotoxic effects of quinolones. It has been suggested that the Central Nervous System stimulating effects of quinolones could be related to their inhibitory effect on γ-amino-butyric acid (GABA) binding [START_REF] Tsuji | Inhibitory effects of quinolone antibacterial agents on gamma-aminobutyric acid binding to receptor sites in rat brain membranes[END_REF], with subsequent decrease of GABAA-mediated inhibitory neurotransmission. The GABA antagonistic effects of quinolones are possibly related to the GABA-like structures of their substituent at position 7 of the heterocycle, and may explain, at least in part, their epileptogenic activity [START_REF] Akahane | Structure-epileptogenicity relationship of quinolones with special reference to their interaction with gamma-aminobutyric acid receptor sites[END_REF]. Quinolones were also shown to interact with the benzodiazepine-GABAA-receptor complex with increased Central Nervous System activity [START_REF] Unseld | Possible interaction of fluoroquinolones with the benzodiazepine-GABAA-receptor complex[END_REF]. Even though the affinity of quinolones for the GABAergic system is weak, a recent study reported that repeated administration of ciprofloxacin in rats could induce anxiety-like behavioral effects in the elevated plus-maze model [START_REF] Ilgin | Ciprofloxacin-induced neurotoxicity: evaluation of possible underlying mechanisms[END_REF]. In that study, the brain GABA content in ciprofloxacintreated animals was decreased when compared with the control group. These effects are in agreement with the potential induction of anxiety by quinolones, which might represent an additional facilitating factor for suicidal behavior in some individuals at risk.

Quinolones do not bind to glutamate receptors [START_REF] Dodd | Neurochemical studies on quinolone antibiotics: effects on glutamate, GABA and adenosine systems in mammalian CNS[END_REF]. However, considering their Mg 2+ chelating properties (the ketone and carboxylic acid functional groups tend to chelate divalent ions), it has been speculated that the excitatory potency of quinolones might be based on activation of the N-methyl-d-aspartate (NMDA) receptor by abolishing the Mg 2+ block in the ion channel [START_REF] Schmuck | Determination of the excitatory potencies of fluoroquinolones in the central nervous system by an in vitro model[END_REF]. The direct role of a glutamatergic alteration in suicide is still unclear but its potential implication cannot be excluded, considering the putative anti-suicidal (suppression of suicidal ideations) properties of the NMDA receptor antagonist ketamine [START_REF] Reinstatler | Ketamine as a potential treatment for suicidal ideation: a systematic review of the literature[END_REF].

On the other hand, the quinolone ciprofloxacin has been shown to significantly decrease serotonin levels in the brain with respect to control values and to induce a depression-like behavioral response in the forced swimming test in rats [START_REF] Ilgin | Ciprofloxacin-induced neurotoxicity: evaluation of possible underlying mechanisms[END_REF], which is in agreement with many studies linking altered serotonin transmission not only to depression but also to impulsive aggressive behavior [START_REF] Oquendo | The biology of impulsivity and suicidality[END_REF] and suicide risk [START_REF] Oquendo | Toward a biosignature for suicide[END_REF].

In addition, oxidative stress inducing capacity of quinolones has been discussed as a possible mechanism leading to neuropsychiatric adverse effects. It has been reported that glutathione levels and catalase activities were decreased and malondialdehyde levels increased, after repeated administration of ciprofloxacin in rats, suggesting an enhanced oxidative stress with weakened antioxidant defence system [START_REF] Ilgin | Ciprofloxacin-induced neurotoxicity: evaluation of possible underlying mechanisms[END_REF]. Major depression is accompanied by decreased antioxidant status and induction of oxidative pathways [START_REF] Maes | A review on the oxidative and nitrosative stress., O&NS pathways in major depression and their possible contribution to the., neurodegenerative processes in that illness[END_REF]. It was also reported that oxidative stress and lowered antioxidant levels may play a role in the pathophysiology of suicidal behaviour [START_REF] Vargas | Oxidative stress and lowered total antioxidant status are associated with a history of suicide attempts[END_REF].

Finally, microRNAs (miRNAs) may have a potential role in some of the central effects of quinolones. MicroRNAs are small non-coding RNAs that regulate gene function, particularly by repressing the translation of target mRNA; these molecules are also involved in cell survival, synaptogenesis and neural plasticity. Alterations of various miRNAs have been implicated in depressive disorders. A global down-regulation of miRNA expression in prefrontal cortex of depressed suicide subjects has been shown [START_REF] Smalheiser | Micro RNA expression is down-regulated and reorganized in prefrontal cortex of depressed suicide subjects[END_REF]. Interestingly, miR-132 and miR-182 negatively regulate Brain-Derived Neurotrophic Factor (BDNF) expression, and their serum levels were found to be increased in patients suffering from depression [START_REF] Li | Alterations of serum levels of BDNF-related miRNAs in patients with depression[END_REF]. Given at doses within the usual anti-bacterial therapeutic range, enoxacin was shown to be quite effective in raising miRNA levels in the cortical area of the brain in a dose-dependent manner [START_REF] Smalheiser | Enoxacin elevates microRNA levels in rat frontal cortex and prevents learned helplessness[END_REF]; however, in these animals, the drug decreased the proportion of rats exhibiting learned helplessness behavior following inescapable shock, which is a common rat model of depressive behavior. These results are not unequivocal and their relevance to the mechanism of suicidal behavior linked to quinolones still remains to be elucidated.

Overall, the many and diverse neurobiological actions of quinolones discussed above may strengthen the biological plausibility of the link between their use and the onset of suicidal behavior.

Strengths and limits of the study

Our study has several strengths. First, we used the global adverse drug reaction database VigiBase which collects information on adverse drug reactions from nearly all national pharmacovigilance centers in the world, and since nearly the origin of their pharmacovigilance systems. This large pharmacovigilance database allowed us to estimate the RORs with precision and good applicability to the global population. Adjusting for age and gender, ROR for suicidal behaviors were consistently statically significant for the main used quinolones. Second, ROR estimations were getting more precise over the years and showed stability around the value of 3 in the last 4 to 5 years, without evidence for a notoriety effect. Third, the results remained consistent after analyzing events without dermatological and/or musculoskeletal effects suggesting the absence of modifying effect by eventcompetition bias. Finally, the significant disproportionality signal for suicidal behavior was consistent across countries, except in only one, Germany. Complementary queries to VigiBase showed that, whereas Germany represented 16% of all adverse reactions reported in Europe, it represented 46% of European reports of suicidal behavior (and 63% for suicidal behaviors associated with antibiotics). Moreover, Germany was also the European country with the most numerous reported adverse effects for moxifloxacin and ofloxacin. We could hypothesize that the absence of signal observed in Germany may be either due to a specific effect of competition between quinolones and other drugs associated with suicide [START_REF] Pariente | A potential competition bias in the detection of safety signals from spontaneous reporting databases[END_REF], or due to different practices for reporting of voluntary drug poisoning, or may represent a chance finding.

The limits of this study include the issue of under-reporting of drug adverse reactions, which is a major drawback inherent to studies using reported adverse events databases [START_REF] Martin | Underreporting of suspected adverse drug reactions to newly marketed, "black triangle" drugs in general practice: observational study[END_REF][START_REF] Bégaud | Rates of spontaneous reporting of adverse drug reactions in France[END_REF][START_REF] Van Der Heijden | On the assessment of adverse drug reactions from spontaneous reporting systems: the influence of underreporting on odds ratios[END_REF]. Under-reporting maybe all the more important than this study concerned non-recent and commonly used drugs. However, suicidal behaviors or depressive disorders are generally considered serious events (which are usually less under-reported) and there is no apparent reason to suggest that these adverse reactions were differentially reported with quinolones in comparison to other antibiotics drugs. Moreover, the absence of safety information by regulatory agencies, the absence of published evidence for this specific risk and our temporal analysis suggest that no notoriety effect could have affected our results. Second, the large amount of ICSRs selected from VigiBase (nearly one million) limited the extraction of all detailed information. For instance, we could not perform our analyses on serious adverse reactions only, it was not possible to analyze the time to onset of the effect after exposure to quinolone, and we could not adjust our analyses for potential confounders of suicidal behaviors or depressive disorders such as personal and family medical history and other drug exposures. We only adjusted our analyses for age and gender, nonetheless, a previous disproportionality analysis of suicidal behaviors performed in VigiBase only in psychiatric events showed that additional adjustment for reporting year, reporting region, reporter type, use of antidepressant and antipsychotic drugs did not substantially change the ROR estimations (data not shown). Third, the high proportion of missing data in spontaneous reports and the substantial risk of data misclassification could have led to information bias. In particular, the analysis of completed suicides may suffer from misclassification as the selection using the MedDRA terms corresponding to completed suicide did not include the cases of suicidal behaviors not coded as completed suicides but resulting in the death of the patient.

Figure 2 .Figure 3 .

 23 Figure 2. Evolution of reported suicidal behaviors (including completed suicide) and completed suicides associated with quinolones from 2000 to 2015

Table 1 . Disproportionality analysis of suicidal behaviors associated with quinolones compared to other antibiotic drugs. Drug Cases of suicidal behaviors Non-cases Crude ROR (95% CI) Adjusted ROR* (95% CI)

 1 Adjusted for age and gender ** This group includes the following drugs: alatrofloxacin/trovafloxacin, enoxacin, fleroxacin, gatifloxacin, lomefloxacin, norfloxacin, pefloxacin, sparfloxacin, temafloxacin, nalidixic acid, and the generic term « fluoroquinolones »

	All quinolones	608	171 927	2.84 (2.57-3.14)	2.78 (2.51-3.08)
	Ciprofloxacin	269	53 792	4.02 (3.51-4.60)	4.01 (3.50-4.59)
	Levofloxacin	166	38 678	3.47 (2.95-4.09)	3.44 (2.91-4.07)
	Moxifloxacin	93	22 498	3.32 (2.68-4.11)	2.84 (2.28-3.53)
	Ofloxacin	71	16 831	3.39 (2.66-4.31)	3.26 (2.56-4.16)
	Other quinolones**	24	41 915	0.46 (0.31-0.69)	0.45 (0.30-0.68)
	Other antibiotics	1019	818 543	1 [Reference]	1 [Reference]

*

Figures Figure 1. Numbers of reported cases of suicidal behaviors and corresponding non-cases according to exposure to quinolone

		Reports involving at least one
		antibiotic recorded in VigiBase
		between 1970 and 2015
		(n = 992,097)
	Cases of suicidal behavior	Other adverse drug reactions
		(n = 1627)	(non cases)
			(n = 990,470)
			)
	Exposed to	Exposed to other	Exposed to	Exposed to other
	quinolones	antibiotics	quinolones	antibiotics
	(n = 608;	(n = 1019;	(n = 171,927;	(n = 818,543;
	37.4%)	62.6%)	17.4%)	82.6%)

Acknowledgments

The Uppsala Monitoring Centre has provided the data but the study results and conclusions are those of the authors and not necessarily those of the Uppsala Monitoring Centre, National Centres, or WHO.

Research ethics

Since this study did not include individual patients and only studied the data of an existing and anonymously recorded database, ethical approval was not needed.

Conclusions

Although uncommon, suicidal behaviors may be observed during treatment with quinolones. Even though disproportionality signals could not be automatically associated with increased incidence of the ADR of interest in the population, the results of this study and the potential mechanisms which are discussed clearly support the hypothesis of an increased risk of suicidal behavior as well as depressive disorders with the use of quinolones drugs when compared to other antibiotics. They highlight the need for psychiatrists and physicians to thoroughly consider recent treatment history in their assessment of patients with suicidal behaviors. As quinolones are commonly used and sometimes overused drugs [START_REF] Werner | Unnecessary use of fluoroquinolone antibiotics in hospitalized patients[END_REF], further investigations, such as large pharmacoepidemiologic studies, are urgent to confirm and better understand this safety signal.

Contributors

All authors participated in the study design. JLF acquired the data. JS did the statistical analyses and wrote the initial draft. All authors critically revised the manuscript.


Jie Yang

Dr Zhen-Yu Yin
email: zhenyu.yin@polyu.edu.hk

Farid Laouafa

Pierre-Yves Hicher

Three-dimensional hydro-mechanical modelling of internal erosion in dikeon-foundation

Keywords: granular media, dike, internal erosion, critical state, fines content, stabilized finite element method

Currently, numerical studies at the real scale of an entire engineering structure considering internal erosion are still rare. This paper presents a 3D numerical simulation of the effects of internal erosion within a linear dike located on a foundation. A 2D finite element code has been extended to 3D in order to analyze the impact of internal erosion under more realistic hydro-mechanical conditions. The saturated soil has been considered as a mixture of four interacting constituents: soil skeleton, erodible fines, fluidized fine particles, and fluid. The detachment and transport of the fine particles have been modeled with a mass exchange model between the solid and the fluid phases. An elasto-plastic constitutive model for sand-silt mixtures has been developed to monitor the effect of the evolution of both the porosity and the fines content induced by internal erosion upon the behavior of the soil skeleton. An unsaturated flow condition has been implemented into this coupled hydromechanical model to describe more accurately the seepage within the dike and the foundation. A stabilized finite element method was used to eliminate spurious numerical oscillations in solving the convection-dominated transport of fluidized particles. This numerical tool was then applied to a specific dike-on-foundation case subjected to internal erosion induced by a leakage located at the bottom of the foundation. Different failure modes were observed and analyzed for different boundary conditions, including the significant influence of the leakage cavity size and the elevation of the water level at the upstream and downstream sides of the dike.

Introduction

Internal erosion occurs when fine particles are plucked off by hydraulic forces and transported through a coarse porous matrix. The known causes are either a concentration of leak erosion, backward erosion, soil contact erosion, or suffusion [START_REF] Bonelli | On the modelling of internal soil erosion[END_REF][START_REF] Fell | Time for development of internal erosion and piping in embankment dams[END_REF][START_REF] Wan | Investigation of rate of erosion of soils in embankment dams[END_REF] . This paper focuses on suffusion, which corresponds to the detachment and migration of fine particles within the voids (pores) of coarse particles by means of seepage flow. Due to the loss of a part of the solid phase, the mechanical properties of the soil are progressively modified and degraded, and with time, earthen structures may be subjected to risk of failure [START_REF] Fell | Time for development of internal erosion and piping in embankment dams[END_REF][START_REF] Xu | Breaching parameters for earth and rockfill dams[END_REF][START_REF] Zhang | Seepage failure mechanism of the Gouhou rockfill dam during reservoir water infiltration[END_REF][START_REF] Zhang | Analysis of earth dam failures: A database approach[END_REF][START_REF] Rönnqvist | On the use of empirical methods for assessment of filters in embankment dams[END_REF] . Foster et al. [START_REF] Foster | The statistics of embankment dam failures and accidents[END_REF] showed that nearly 46% of the damage of 128 embankment dams could be attributed to internal erosion. Moreover, sinkholes and cavities triggered by internal erosion are frequently observed within dams and dikes [START_REF] Wood | The magic of sands-the 20th Bjerrum Lecture presented in Oslo, 25 November 2005[END_REF][START_REF] Sterpi | Effects of the erosion and transport of fine particles due to seepage flow[END_REF] . Internal erosion can also trigger slope failures and generate significant landslides [START_REF] Crosta | On slope instability induced by seepage erosion[END_REF][START_REF] Hu | Seismic precursor to instability induced by internal erosion in loose granular slopes[END_REF] .

Laboratory experimental tests [START_REF] Sibille | A description of internal erosion by suffusion and induced settlements on cohesionless granular matter[END_REF][START_REF] Rochim | Effects of Hydraulic Loading History on Suffusion Susceptibility of Cohesionless Soils[END_REF][START_REF] Moffat | Spatial and temporal progression of internal erosion in cohesionless soil[END_REF][START_REF] Chang | A stress-controlled erosion apparatus for studying internal erosion in soils[END_REF][START_REF] Marot | Assessing the susceptibility of gap-graded soils to internal erosion: proposition of a new experimental methodology[END_REF] and numerical analyses [START_REF] Vardoulakis | Hydro-mechanical aspects of the sand production problem[END_REF][START_REF] Stavropoulou | Coupled wellbore erosion and stability analysis[END_REF][START_REF] Papamichos | Volumetric sand production model and experiment[END_REF][START_REF] Cividini | Finite-element approach to the erosion and transport of fine particles in granular soils[END_REF][START_REF] Lominé | Modeling of fluid-solid interaction in granular media with coupled lattice Boltzmann/discrete element methods: application to piping erosion[END_REF][START_REF] Mansouri | Numerical simulation of the quicksand phenomenon by a 3D coupled Discrete Element-Lattice Boltzmann hydromechanical model[END_REF][START_REF] Zhang | Three-dimensional DEM modeling of the stress-strain behavior for the gap-graded soils subjected to internal erosion[END_REF] at the scale of the representative elementary volume have undoubtedly played an important role in understanding the physical mechanisms behind internal erosion phenomena. However, these studies have not provided a clear understanding of where and when seepage flow could trigger internal instability in a granular soil mixture at larger space scales for real boundary value problems. Several physical tests have been carried out for this purpose 12, 25, 26 . However, few numerical studies at the scale of an engineering structure are available in the literature [START_REF] Vardoulakis | Hydro-mechanical aspects of the sand production problem[END_REF][START_REF] Stavropoulou | Coupled wellbore erosion and stability analysis[END_REF][START_REF] Cividini | Finite-element approach to the erosion and transport of fine particles in granular soils[END_REF][START_REF] Zhang | A thermodynamics-based model on the internal erosion of earth structures[END_REF][START_REF] Uzuoka | Hydro-mechanical analysis of internal erosion with mass exchange between solid and water[END_REF] . In the available studies, the mechanical responses of an eroded soil were often bypassed or simplified by using an elastic model with induced damage parameters. Yang et al. [START_REF] Yang | Internal erosion in dike-on-foundation modeled by a coupled hydromechanical approach[END_REF] proposed a coupled hydro-mechanical approach based on porous continuous medium theory to assess how internal erosion impacts the safety of earthen structures. They developed an elasto-plastic constitutive model for sand-silt mixtures to monitor the effect of the evolution of both the porosity and the fines content induced by internal erosion upon the behavior of the soil skeleton. This model was then applied to analyse internal erosion within a fully saturated dike under plane strain condition. However, a three-dimensional (3D) condition has long been thought to be a significant factor in the internal erosion process. Under 3D conditions, the constraint of the surrounding soil may limit or delay the development of the soil deformation, which may in turn introduce a more 'local' failure-collapse mode (i.e., a sinkhole) instead of a landslide-type instability.

The aim of this study is to simulate the 3D mechanical consequences of internal erosion within a linear dike on a foundation under more general and realistic conditions. First, the unsaturated flow condition was implemented into the coupled hydro-mechanical model of Yang et al. [START_REF] Yang | Internal erosion in dike-on-foundation modeled by a coupled hydromechanical approach[END_REF] in order to describe more accurately the seepage within the dike and foundation due to the difference in water pressure at the upstream and downstream sides of the dike-onfoundation. Afterwards, the enhanced model was implemented into a 3D finite element code using a stabilized finite element method. Then it was applied to simulate the effect of internal erosion induced by a seepage flow within a dike, and a leakage at the bottom of the foundation due to the presence of a cavity. The influence of the leakage cavity size and the elevation of the water level at the upstream and downstream sides of the dike-on-foundation were thus analysed.

Coupled hydro-mechanical model 2.1. Mass balances and transported particles

According to Vardoulakis et al. [START_REF] Vardoulakis | Hydro-mechanical aspects of the sand production problem[END_REF] and Schaufler et al. [START_REF] Schaufler | Infiltration processes in cohesionless soils[END_REF] , it is possible to consider a saturated porous medium as a material system composed of 4 constituents in 2 phases: the stable fabric of the solid skeleton, erodible fines, fluidized particles and pure fluid. The fines can behave either as a fluid-like (described as fluidized particles) or as a solid-like (described as erodible fines) material. The hydraulic process of suffusion can be described by a system of non-linear partial differential equations deduced from the mass balance of the mixture system (More details can be found in [START_REF] Yang | Internal erosion in dike-on-foundation modeled by a coupled hydromechanical approach[END_REF] ):

    div div 0 w s   q v ...
        div div c c c s c s f f f f n t t           v v ..
      div s c c c n t t          w v q ..

Erosion law

The variable n in Eqs. ( 2)-( 4) is the volumetric mass exchange at any point in time, which determines the rate of erosion. It should be defined by an appropriate erosion law. Since the proposed numerical approach is modular and independent of the erosion law, it is possible to choose the most appropriate expression of the erosion law for different cases. Due to the lack of experimental data, the present paper is mainly a theoretical study. Its focus is not upon a comparison with a real case and it does not discuss the choice of the erosion law. A simple erosion law [START_REF] Uzuoka | Hydro-mechanical analysis of internal erosion with mass exchange between solid and water[END_REF] has therefore been adopted as follow:  is a material parameter, experimentally determined. Furthermore, in the following analyses, a residual fines content has been artificially introduced for the sake of simplicity.

ˆe c n f    w q ....

Unsaturated flow in dike-on-foundation

For isotropic partially saturated porous continua, the classical expression for the Cauchy effective stress tensor  σ consistent with continuum principles of thermodynamics can be expressed as follows [START_REF] Borja | On the mechanical energy and effective stress in saturated and unsaturated porous continua[END_REF] : where σ is the total Cauchy stress tensor, B is the Biot coefficient, I is the second-order identity tensor and p is the mean pore pressure; w p and a p are the pore water and pore air pressures, and w s and a s are the degrees of water saturation and air saturation, respectively.

Bp   
The tensile stresses are considered positive, which explains the sign in Eq. ( 9). We consider in the following a Biot coefficient B equal to one [START_REF] Uzuoka | Dynamics of unsaturated poroelastic solids at finite strain[END_REF] . In a transient and quasi-static analysis, the model can be simplified by assuming that the air pressure is constant throughout the domain and is small enough so that its value can be disregarded. , the saturation condition lies within certain limits [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF][START_REF] Mckee | Flow-testing coalbed methane production wells in the presence of water and gas[END_REF] . In this study, the empirical soil-water characteristic curve (SWCC) proposed by van Genuchten [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] , which has been proved to be suitable for silty sand [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] , was adopted to explicitly determine the degree of water saturation for a given capillary pressure.

In the van Genuchten model, the SWCC is expressed as follows [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] : s was set to be 1.0.

    , , , , , 1 vg
The water permeability coefficients within the van Genuchten model were assumed to be dependent on effective water saturation as: Therefore, the unsaturated flow could be introduced into the hydro-mechanical model for internal erosion by combining Eqs. ( 11)- (13) with Eqs. ( 1)-(4). A critical state based constitutive model was adopted to determine the Cauchy effective stress. This model has the capacity to consider the impact of the change of porosity and fines content upon the mechanical behavior of the soil during the erosion process, as described below. Note that in this study, the foundation soil was considered as a sandy soil with coarse grains, and thus, the above modification to the effective stresses was suitable when applied directly to the saturated sand model.

      2 1 , , 1

Critical state based soil model considering fines effect

To assess the mechanical response to the loss of fines content induced by internal erosion in an earthen structure, it was decided to adopt a non-associated elastoplastic constitutive model for sand-silt mixture for the solid skeleton [START_REF] Yin | Elastoplastic modeling of sand-silt mixtures[END_REF] . This model was extended from the basic SIMSAND model [START_REF] Jin | Investigation into MOGA for identifying parameters of a critical-state-based sand model and parameters correlation by factor analysis[END_REF][START_REF] Wu | A straightforward procedure of parameters determination for sand: a bridge from critical state based constitutive modelling to finite element analysis[END_REF] under the framework of the critical state concept and elastoplasticity theory (summarized in Appendix) by defining the position of the critical state line as a function of the fines content in order to unify the mechanical behavior of a sand-silt mixture from silt to sand or sand to silt (see ( 14) and ( 15)). Therefore, this model takes into account the influence of the porosity and the fines content evolution on the mechanical behavior of the soil during the erosion process. It should be noted that for sandy soils, only the shear yield surface is needed.

A non-linear expression of the critical state line is adopted here, and has the following expression: Given the evolution of the fines content during internal erosion, the constitutive model needs to take into account such internal soil mass evolution. Yin et al. [START_REF] Yin | Elastoplastic modeling of sand-silt mixtures[END_REF][START_REF] Yin | A micromechanics-based model for sand-silt mixtures[END_REF] considered 0 cr e as a function of the fines content in order to unify the mechanical behavior of a sand-silt mixture for different fines contents:

            0 , 0 , 0 1 tanh 1 tanh 1 1 2 2 c t h c t h c cr hc cr c c hf cr c m d f f f f f e e f af e f R                         ....... ( 15
)
The function was mathematically structured based on the hyperbolic tangent function

    2 2 tanh 1 1 x x
x e e    to ensure a smooth transition between a silty sand and a sandy silt.

, 0 hc cr e and , 0

hf cr e are the reference critical void ratios for pure sand and pure silt, respectively.

a ,  , d R , m and th f are material constants. More details can be found in Yin et al. [START_REF] Yin | Elastoplastic modeling of sand-silt mixtures[END_REF][START_REF] Yin | A micromechanics-based model for sand-silt mixtures[END_REF] .

A finite difference code was developed to investigate the validity of the proposed model [START_REF] Yang | Modeling coupled erosion and filtration of fine particles in granular media[END_REF][START_REF] Yang | Hydro-mechanical modeling of granular soils considering internal erosion[END_REF][START_REF] Yang | Analysis of suffusion in cohesionless soils with randomly distributed porosity and fines content[END_REF] .

The comparison between experimental results and simulations demonstrated that, with an appropriate erosion law and well calibrated model parameters, the numerical model was able to reproduce with success the initiation and development of internal erosion.

Finite element approach

Standard Galerkin-weighted residual method

The following governing equations were obtained by combining the balance and constitutive x . For the sake of legibility, the time t and the space x variables have been omitted in the equations.

      div div + 0 s s ec w v v f q t          ....
This coupled mechanical-erosion process is a non-linear transient problem. Weak forms of the governing partial differential equations were implemented into the finite element code ABAQUS [START_REF] Hibbitt | ABAQUS/standard User's Manual[END_REF] . ABAQUS [START_REF] Hibbitt | ABAQUS/standard User's Manual[END_REF] where AMATRX and RHS are the Jacobian matrix and the residual nodal fluxes or forces needed to be defined corresponding to Eqs. ( 16)- (20); N d is the nodal vector of the DOFs   u , u , u , , , , q , q , q

Stabilized finite element method

Several problems were noted in the treatment of the convection (advective) dominated flows (Eq. ( 19)). Numerical oscillation occurred when the standard Galerkin method was applied to convection dominated flows, affecting the computation accuracy and efficiency [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF][START_REF] Codina | Comparison of some finite element methods for solving the diffusionconvection-reaction equation[END_REF] . The non-linear convection terms render the governing differential equations non-self-adjoint [START_REF] Sampaio | A Petrov-Galerkin/modified operator formulation for convectiondiffusion problems[END_REF] .

The matrix associated with the convective term is non-symmetric. As a result, the "best approximation" [START_REF] Surana | The finite element method for boundary value problems: mathematics and computations[END_REF] property is lost, resulting in numerical oscillations.

In practice, the oscillations become more significant as the Péclet number increases, which is a function of the properties of the porous medium governing diffusion and velocity of the fluid, as well as the finite element size [START_REF] Donea | Finite element methods for flow problems[END_REF][START_REF] Al-Khoury | Computational modeling of shallow geothermal systems[END_REF][START_REF] Zienkiewicz | The finite element method for solid and structural mechanics[END_REF] . Solutions are often corrupted by spurious oscillations when a downstream boundary condition forces a rapid change in the solution.

One way to eliminate such numerical oscillations is to refine the mesh, such that convection no longer dominates at an element level [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] . However, the prohibitive computational cost in 3D condition becomes an obstacle whenever real soil structures such as embankment dams are studied. This has provided motivation for an alternative to the Galerkin formulation that precludes spurious oscillations regardless of mesh refinement.

In the last decades, considerable effort has been made towards developing a stabilized finite element method to improve the numerical solutions of the heat or solute transport problems

Converged ? No

Yes [START_REF] Yu | Petrov-Galerkin methods for the time-dependent convective transport equation[END_REF][START_REF] Yu | Petrov-Galerkin method for multidimensional, time-dependent, convectivediffusion equations[END_REF][START_REF] Dick | Accurate Petrov-Galerkin methods for transient convective diffusion problems[END_REF][START_REF] Hughes | A new finite element formulation for computational fluid dynamics: II. Beyond SUPG[END_REF][START_REF] Tezduyar | Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations[END_REF][START_REF] Hughes | A new finite element formulation for computational fluid dynamics: VIII. The galerkin/least-squares method for advectivediffusive equations[END_REF][START_REF] Codina | The intrinsic time for the streamline upwind/Petrov-Galerkin formulation using quadratic elements[END_REF][START_REF] Al-Khoury | Efficient finite element formulation for geothermal heating systems. Part II: transient[END_REF][START_REF] Cui | Numerical analysis of coupled thermo-hydraulic problems in geotechnical engineering[END_REF] . The present paper does not focus upon a comparison of different algorithms. The Petrov-Galerkin weighted residual method proposed by Yu and Heinrich [START_REF] Yu | Petrov-Galerkin methods for the time-dependent convective transport equation[END_REF][START_REF] Yu | Petrov-Galerkin method for multidimensional, time-dependent, convectivediffusion equations[END_REF] was employed for solving Eq. ( 19). The method involves perturbed weighting functions in the weighted residual formulations that depend on two parameters ( ,  ) calculated according to a local analysis. The perturbations to the weighting functions can be interpreted as an added local anisotropic balancing diffusion and an added dispersion in the direction of the convective motion. The integration of the weak form of Eq. ( 19) from time t to t t   provides an average equilibrium statement for the increment, which gives: where N N are shape functions, v is the average fluid velocity over the element; v is its Euclidean norm; h is a characteristic element length measure, defined as: where  h is the  isoparametric line across the element, passing through its centroid. h  is the projection of  h in the direction of the fluid velocity vector at the element centroid.

      , , , , , , 1 2 1 2 2 4 ˆ0 2
N N M M t t M t N M N M t tM t N M M t t M t N N h dV c c t h dV c c h dV c c h n c dV t                                                                      v N N N v x v N N N v v x x v N N v v x x v N N v
In Eq.( 22),  and  are control parameters. The  term in the weighting function was introduced to eliminate artificial diffusion of the solution, while the  term was introduced to avoid numerical dispersion, given as [START_REF] Yu | Petrov-Galerkin methods for the time-dependent convective transport equation[END_REF][START_REF] Yu | Petrov-Galerkin method for multidimensional, time-dependent, convectivediffusion equations[END_REF] : It should be noted that the deformation terms have been omitted from the equilibrium equation ( 22) for the sake of simplicity, which assumes that the deformation of the specific domain will not change the concentration of the fluidized fine particles. The boundary conditions were taken as symmetry plane (rollers) on the left and right boundary surfaces of the foundation and as fully fixed at the base of the foundation. The gravity load was initially applied to generate the initial stresses and water pressure fields.

2 coth 2      and 2

Three

Afterwards, the water head at the cavity was set to 11m in order to initiate water flow in and out the cavity. The alluvium and the dike were considered to be made of the same soil: a mixture of sand and silt. The physical properties and the material parameters, summarized in . The residual fines content was artificially fixed equal to 0.05 for the sake of simplicity in the 3D analysis.

Due to the lack of experimental data, the set of the van Genuchten model parameters used by Uzuoka et al. [START_REF] Uzuoka | Hydro-mechanical analysis of internal erosion with mass exchange between solid and water[END_REF] was adopted, as shown in Table 3. Figure 6 shows the SWCC and the relationship between the relative permeability of water and the degree of water saturation given by the van Genuchten model with the parameters in Table 3.

Plane of symmetry

D=0.2m

Plane of symmetry [START_REF] Jin | Multi-objective optimization-based updating of predictions during excavation[END_REF][START_REF] Yin | Optimization techniques for identifying soil parameters in geotechnical engineering: comparative study and enhancement[END_REF][START_REF] Yin | An efficient optimization method for identifying parameters of soft structured clay by an enhanced genetic algorithm and elastic-viscoplastic model[END_REF][START_REF] Jin | Selection of sand models and identification of parameters using an enhanced genetic algorithm[END_REF] or Bayesian based methods [START_REF] Jin | Bayesian model selection for sand with generalization ability evaluation[END_REF][START_REF] Jin | Identifying parameters of advanced soil models using an enhanced transitional Markov chain Monte Carlo method[END_REF] .

Hydro-mechanical responses in 3D condition

Due to the symmetry of the domain and the mechanical and hydraulic boundary conditions, only half the domain was modelled. Figure 7 shows the initial space distribution of the pore pressure in the dike and foundation under the gravity load. The phreatic surface corresponds to the locus of the points at which the pore water pressure is zero. Figure 8 shows the space distribution of the fines content after erosion for Case 1. As internal erosion took place, the fines content which was initially spatially uniform throughout the whole soil volume began to decrease with time. From a numerical point of view, the erosion process is more active where the hydraulic gradient changes more drastically. These analyses in the present paper were theoretical, they were carried out at constant water head and the materials of the dike and the foundation were assumed to be initially homogeneous. The water heads at the upstream and downstream were pre-fixed. The flux was therefore larger along the phreatic surface, especially at the two ends of the phreatic surface. Thus, for this theoretic case, it was reasonable to find a greater loss of fines at the top of the dike and at the bottom of the system along the phreatic surface. The water head may vary slowly or rapidly as a result of rainfall or human intervention. The seepage path and the location of the phreatic surface may become more complex and could evolve with time within a heterogeneous medium. The eroded zone may change accordingly. It can be concluded that the shear sliding and the volumetric settlements are two major consequences of the internal erosion. These phenomena depend, to a certain extent, on the location of the eroded zones. This aspect, under 2D condition, was discussed in a previous work [START_REF] Yang | Internal erosion in dike-on-foundation modeled by a coupled hydromechanical approach[END_REF] . In plane strain condition, a significant settlement of the dike was observed when the cavity was located below the top of the dike. In contrast, the results of the 3D simulations showed that the constraint of the surrounding soil limited the development of the deformation, which in turn introduced a local failure (i.e. the sinkhole shown in Figure 11-Figure 13) before the global sliding failure occurred.

It should be noted that, in this study, a large value of the erosion parameter e  was artificially chosen to reduce the simulation time. However, the time of failure is closely related to the erosion parameters and the permeability of the soil, both of which need to be determined experimentally [START_REF] Marot | Assessing the susceptibility of gap-graded soils to internal erosion: proposition of a new experimental methodology[END_REF][START_REF] Vardoulakis | Hydro-mechanical aspects of the sand production problem[END_REF] with care.

Influences of cavity size and water head

Effect of the size of the leakage cavity

In Case 3, a larger leakage cavity of 1.0 m in diameter was established at the bottom of the symmetry plane of the model, as shown in Figure 14, in order to investigate how the size of the leakage cavity influences the development of internal erosion within the dike and foundation. The observations of the Val d'Orléans dike showed that the sinkholes developed at different locations (Figure 22). The formation of these sinkholes was assumed to be due to the presence of karstic cavities below the soil forming the foundation of the dike. It can be seen that the numerical simulations were able to reproduce qualitatively the initiation and the development of these sinkholes whose locations within the dike are strongly related to the seepage conditions. will surely improve the understanding of the origin of damages or failures induced by internal erosion and serve as a reference to update the specifications for design, construction, and adaptation in order to mitigate and prevent the highest risks.

Conclusion

This study is an attempt to provide a novel method to identify the potential consequences of internal erosion developed within a dike under different boundary conditions. Even if the analyzed cases were mainly theoretical, they bear a close resemblance to real case problems.

The newly developed hydro-mechanical model for internal erosion was enhanced by introducing unsaturated flow in order to describe the seepage within the dike and foundation due to the difference in water pressure at the upstream and downstream sides. This enhanced model was implemented into a 3D finite element code. A stabilized finite element method was used to eliminate spurious numerical oscillations to solve the convection dominated flow of fluidized particles.

Then, a series of 3D simulations of the development of internal erosion within the dike and the foundation were carried out. It was observed that the decrease of the fines content close to the toe of the dike could stimulate a global sliding failure of the slope. The elevation of the water surface within the foundation at the downstream side increased the loss of the fines close to the toe. On the other hand, the presence of a leakage cavity could lead to the initiation of a sinkhole formed as a consequence of the rigid movement of the soil above the eroded area in the vicinity of the cavity. The location of the sinkhole was found to be related to the location of the phreatic surface. Moreover, the enlargement of the cavity accelerated the formation of the sinkhole. For a larger leakage cavity, the internal erosion within the dike and foundation developed faster, leading to the earlier appearance of a sinkhole at the top of the dike. In terms of failure mode, the presence of a leakage cavity introduced a more 'local' failure-collapse mode (i.e. a sinkhole) instead of a global sliding failure which was obtained by the numerical simulations performed under plain strain condition. Furthermore, the elevation of the water surface within the foundation at the downstream side increased the risk of a global sliding of the slope.

It is worth mentioning that the erosion model was developed for a saturated medium. In the unsaturated zone, internal erosion is assumed to be negligible due to the absence of seepage forces. However, the development of capillary forces in the unsaturated zone increases the mechanical resistance of the soil. This phenomenon will affect the global mechanical response of the dike. Besides, these analyses were carried out at constant water head and the materials of the dike and the foundation were assumed to be initially homogeneous. In nature, the water head may vary, slowly or rapidly, as a result of rainfall or human activities. The seepage path and the location of the phreatic surface may become more complex and could evolve with time within a heterogeneous medium. The eroded zone may change accordingly.

These aspects should be analyzed in future studies.

  ....................................................................................... are the porosity, the amount of erodible fines and the concentration of the fluidized particles, respectively; n is the source term describing the exchange between the erodible fines and the fluidized particles; s v is the velocity of the soil skeleton and   , w x t q is the total discharge of the pore fluid assumed to be governed by Darcy's law: ......................................................................................................... (5) with the intrinsic permeability of the medium k and the density of the mixture   ...................................................................................................................... (6) ............................................................................................................. (7) where k  is the kinematic viscosity of the fluid, w p is the pore fluid pressure, s  and f  are the density of the solid and the fluid, respectively, and 0 k is the permeability parameter of Kozeny-Carman relationship 31 .

  ............................................................................................................. (14) where 0 cr e is the reference critical void ratio corresponding to the void ratio e at the mean effective stress ' 0 p  . 0 cr e determines the position of the CSL in the ' pressure.  and  are material constants controlling the non-linearity of the critical state line.

  equations in order to formulate the Initial Boundary Value Problem (IBVP) of internal erosion: ........................................................................................................................ (16) ......................................................................................................... (17)

  provides a user subroutine option (UEL) that allows users to define a new type of element with user's defined governing equations and degrees of freedom (DOFs). The subroutine makes it possible to use advantageously the powerful pre-and postprocessing and highly efficient matrix solver of ABAQUS. The global calculate flow of ABAQUS with the call of UEL subroutine is shown in Figure1where the part of "from Call UEL to update RHS" was developed in this study. In general, Abaqus/UEL solves the overall system of non-linear equations by Newton's method: N   AMATRX d RHS ........................................................................................................ (21)

  isoparametric brick element with 9 DOFs has been developed and implemented into ABAQUS via the UEL subroutine. The detailed vector of residual and the Jacobian matrix can be found in Yang et al.[START_REF] Yang | Internal erosion in dike-on-foundation modeled by a coupled hydromechanical approach[END_REF]

Figure 1 .

 1 Figure 1. Global calculate flow of ABAQUS with the call of UEL subroutine

  ................................................................................................ (23)

  ..................................................................................... (24) ................................................................................................... (25) where  is the local Péclet number in an element and C is the local Courant number.

-dimensional modeling of internal erosion within a dike 4 . 1 . 1 S in Figure 2 )

 4112 Description of the problemThe characteristics of the studied configuration are presented in Figure2according to the case study of the dike of Val d'Orléans[START_REF] Alboresha | Evaluation of the impact of a cavity upon an earth dike (analytical and numerical approaches): Application to the Val d'Orléans area[END_REF] . The dike was built on sandy-silt alluvium formation over limestone. The upstream face of the dike-on-foundation (Face is exposed to the water within the river behind the dike. Likewise, the downstream face of the dike-on-foundation (Face 4 S in Figure2) is exposed to the atmosphere. If the soil is partially saturated, the phreatic surface in the dike and foundation (Face 2 S in Figure2) corresponds to the locus of the points at which the pore fluid pressure, w p , is equal to zero. Above this surface, the pore fluid pressure is negative, which represents the capillary tension that causes the fluid to rise against gravitational forces, creating a capillary zone. The saturation associated to specific values of the capillary pressure for the absorption of the fluid by the porous medium is a physical property of the material and has been defined by the van Genuchten model. The bottom of the alluvium layer (Face 3 S in Figure2) was assumed to rest on an impermeable limestone layer. The seepage due to the water pressure gradient at the upstream and downstream faces of the dike-on-foundation could lead to the detachment of fine particles that will move from the upstream to the downstream side, eventually leading to the degradation of the dike. Furthermore, in the case of the Val d'Orléans dike, it has been reported that karst collapses underneath the foundation might lead to a leakage cavity through which underground water and fine particles could flow out. Internal erosion could occur due to local hydraulic gradients in the vicinity of the cavity. The computations aimed to show how internal erosion may weaken the foundation and the dike under different boundary conditions.

Figure 2

 2 Figure 2 Scheme of the dike-on-foundation near Orléans (France) with a cavity located at the bottom.

Figure 3 Figure 4 Figure 5 .(a) Case 1 (b) Case 2 Figure 5

 345125 Figure 3 Schematic configuration of the dike and mesh of the numerical model

Figure 6

 6 Figure 6 (a) SWCC and (b) relative water permeability using van Genuchten modelHowever, it should be underestimated that the present work is a theoretical study. The aim is to provide a novel method to identify theoretically the potential consequences of internal erosion developed within a dike under different conditions. The numerical results were analyzed in a qualitive way. With an appropriate erosion law and well calibrated parameters, the numerical results could be compared with the in-situ data in the future. At the present

Figure 9 and

 9 Figure 9 and Figure 10 show the displacement and the deviatoric plastic strain fields within the dike at the end of the erosion ( time 30  days) for Case 1. The increase of the void ratio accompanied by a downward shift of the CSL in the e-p' plane induced by the decrease of the fines content made the soil increasingly looser. Consequently, the soil strength decreased.This led to the onset of a shear sliding surface in the most eroded area near the phreatic surface at the downstream side. The plastic deformation continued to develop with time until the dike ruptured by a sliding of the slope.

Figure 7 Figure 8 Figure 9 Figure 10

 78910 Figure 7 Spatial distribution of initial pore pressure in the dike and foundation under gravity load

Figure 11 Figure 12 Figure 13

 111213 Figure 11 Displacement magnitude field after erosion for Case 2 (Unit: m)

Figure 14 hydraulic 3 .

 143 Figure 14 hydraulic boundary conditions for Case 3 with the cavity of leakage (Diameter D=1.0 m) Figure 15 shows the displacement of the dike-on-foundation at the end of erosion (

Figure 15 Figure 16 (Figure 17 5 . 2 . 5 Figure 18 Figure 19 -Figure 19 Figure 20 Figure 21

 1516175251819192021 Figure 15 Displacement magnitude field after erosion for Case 3 (Unit: m)

Figure 22

 22 Figure 22 Photos of sinkholes within the Val d'Orléans dike (France) according to Gombert at al. 67 : (a) at the toe of the dike, (b) at the head of the dike, (c) on the slope face of the dike, (d) behind the dike It can be concluded that the development of internal erosion and the induced deformation of the dike-on-foundation are closely related to its hydraulic boundary conditions. Systematic studies about the influence of the hydraulic boundary conditions (i.e. non-constant water head)

Figure 1 .Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11 Figure 12 Figure 13 Figure 15 3 Figure 16 3 Figure 17 Figure 18 Figure 19 5 Figure 20 Figure 21 Figure 22

 123456789101112131531631718195202122 Figure 1. Global calculate flow of ABAQUS with the call of UEL subroutine Figure 2 Scheme of the dike-on-foundation near Orléans (France) Figure 3 Schematic configuration of the dike and mesh of the numerical model Figure 4 Mesh of numerical model (1/2 model due to Figure 5 Boundary hydraulic conditions for: (a) Case1-without the cavity; and (b) Case 2with the cavity (diameter D=0.2 m) Figure 6 (a) SWCC and (b) relative water permeability using van Genuchten model Figure 7 Spatial distribution of initial pore pressure in the dike and foundation under gravity load Figure 8 Spatial distribution of fines content after erosion for Case 1 (1/2 model) Figure 9 Displacement magnitude field after erosion for Case 1 Figure 10 Deviatoric plastic strain field after erosion for Case 1 (1/2 model) Figure 11 Displacement magnitude field after erosion for Case 2 (Unit: m) Figure 12 Deviatoric plastic strain field after erosion for Case 2 (1/2 model) Figure 13 Spatial distribution of fines content after erosion for Case 2 (1/2 model) Figure 14 hydraulic boundary conditions for Case 3 with the cavity of leakage (Diameter D=1.0 m) Figure 15 Displacement magnitude field after erosion for Case 3 Figure 16 Time variations of the maximum settlement (Point O in Figure 3) during erosion for Case 2 and Case 3 Figure 17 Spatial distribution of fines content after erosion for Case 2 and Case 3 (1/2 model) Figure 18 Hydraulic boundary conditions for Case 4 and Case 5 with different water levels downstream and upstream Figure 19 Displacement magnitude field after erosion for Case 2, Case 4 and Case 5 Figure 20 Deviatoric plastic strain field after erosion for Case 2, Case 4 and Case 5 (1/2 model) Figure 21 Spatial distribution of fines content after erosion for Case 2, Case 4 and Case 5 (1/2 model) Figure 22 Photos of sinkholes within the Val d'Orléans dike (France) according to Gombert at al. 67 : (a) at the toe of the dike, (b) at the head of the dike, (c) on the slope face of the dike, (d) behind the dike

  σ σ I ............................................................................................................................ (9) .....................................................................................................................(10)

	w w p s p s p a a  

  ....................................................................................................................(11)

						This simplification reduces
	the effective stress expression to:
	' with   c p σ σ I Therefore, the medium is fully saturated for , w c w s p p  	0 p  ; negative values of w w p represent the
	capillary effects. For	0 p  , it is known that at a given value of the capillary pressure, w
	c p		 	w p	

  1 ................................................................................ (13)

	vg m ...where s vg vg m s we we k k s s           k is the saturated (maximum) coefficient of water permeability and vg  is a material
	parameter.

   ..........................................................................................................(20)

	  c t   w  k  k q   where i div grad   w cq w p w denotes the body force per unit volume. Together with Eqs. (6)-(15), a system of 5   div 0 s ec w c v f q      0 .  equations with 5 unknowns   , , , ,

........................................................................ (18) ....................................................................... (19)



  44,

	Beginning of Analysis
	Define the mesh and
	initial/boundary conditions
		Start of iteration
	Calculate	Call UEL and update ,    	   
	Update Jacobian matrix AMATRX
	Define loads and update RHS
	Solve	N   AMATRX d	RHS
		End of iteration

Table 1 and

 1 Table 2, corresponded to the Ottawa 50/200 sand-silt characteristics 36 with

	0 e 	0.6	, 0 0.01 c 	and 0 0.4

c f 

Table 1

 1 Physical properties of soils

	Density of fluid Density of solids	 f  s	1.0 g/cm 3 2.65	Kinematic viscosity of fluid Initial intrinsic permeability	 k 5.0E-6 m 2 s -1 k 2.5E-9 m 2

Table 2

 2 Material constants of soils

	Elastic parameters G 0 /MP  n	Plastic parameters k p A d  u	CSL-related parameters e hc,cr0 e fc,cr0 		Fines parameters a m f th  	Erosion parameters  e
	a
	150	0.2	0.5	0.0017	1 32	0.805 1.03 0.196 0.081	0 0.7 0.3	0.02

Table 3

 3 Material parameters in van Genuchten model

, w s

Acknowledgement

The financial supports provided by a GRF project (Grant No. 15209119) from Research Grants Council (RGC) of Hong Kong and the National Institute for Industrial Environment and Risks of France (INERIS) are gratefully acknowledged.

Appendix: Brief introduction to the SIMSAND model