Brice Reoyo-Prats

Dominique Aubert

Amélie Sellier

Benoit Roig

Carmen Palacios
email: carmen.palacios@univ-perp.fr.

Dynamics and sources of pharmaceutically active compounds in a coastal Mediterranean river during heavy rains

Keywords: pharmaceutically active compounds, water quality, multicontamination phenomena, coastal rivers, Mediterranean climate, storm events, multiple stressors

Concentrations of Pharmaceutically Active Compounds (PACs) in freshwater systems

depend on numerous factors such as land use and hydrometeorological conditions. In the Mediterranean, heavy rain events are of particular importance as they highly influence the concentration of micropollutants found in freshwater and are a source of recurrent first foul flushes due to Combined Sewer Overflows (CSOs). In this study, we seek to assess the dynamics of pharmaceuticals during storm events in coastal Mediterranean rivers at a fine scale and to determine their contribution to multicontamination phenomena owing to CSOs. Our results showed that, while dissolved PACs followed the same trend as other contaminants i.e. they increased significantly during CSOs, PACs in the total fraction did not peak yet maintained their already high concentrations for slightly longer due to their release via CSOs. Pharmaceutical concentrations for both the dissolved and the total fraction were dramatically diluted during the peak river flow. A fine-scale follow-up of PACs dynamics in the total fraction, including the differentiation of sewer overflows from both the right and left river banks, as well as the analyses of a large amount of PACs molecules, allowed us to clearly identify their major sources. While domestic inputs were dominated by nicotine and caffeine, the use of gadolinium (an MRI contrast agent) as a marker, attributed the main source of medical drugs such as tramadol, ibuprofen and diclofenac to the major public hospital of the region. Thus, identifying major sources of PACs and implementing adapted water treatments directly at those sources would be the most cost-efficient alternative to cope with pharmaceutical drugs in coastal Mediterranean aquatic environments. Moreover, PACs behavior differed depending on the molecules considered and the source of these molecules, but we could not establish a direct link between their behavior and their chemical or physical properties. Our study highlights the importance of monitoring at strategic locations and with a high frequency sampling in order to better understand fate, sources and behavior of pharmaceuticals in aquatic environments.

Introduction

Two decades ago the European Environment Agency, the scientific community and the public started to be concerned about the environmental and health effects of chemically bioactive substances used in human and veterinary medicine found at trace levels (micropollutants) in aquatic environments (EEA 2010). Still, these Pharmaceutically Active Compounds (PACs) were not monitored in the European Union until 2015 and today, only three substances are regularly monitored including diclofenac, 17α-ethinylestradiol, 17β-estradiol and estrone [START_REF] Ec | Establishing a watch list of subtances for Union-wide monitoring in the field of water policy pursuant to Directive[END_REF]. PACs, like other micropollutants such as pesticides or metal elements, are released in the aquatic environment via point sources of pollution such as WasteWater Treatment Plants (WWTP) and hospitals [START_REF] Pal | Impacts of emerging organic contaminants on freshwater resources: Review of recent occurrences, sources, fate and effects[END_REF]. Furthermore, human and veterinary-use drugs accumulate in superficial soils through the use of sewage sludge or manure as fertilizers in farm fields. The remobilization of PACs from soils through runoff and leaching is therefore a potential diffuse source for these compounds [START_REF] Halling-Sorensen | Occurence, fate and effects of pharmaceuticals substance in the environment -A review[END_REF][START_REF] Gielen | Factors impacting on pharmaceutical leaching following sewage application to land[END_REF]). Thus, the concentration of PACs in surface waters will depend on both watershed land use and rainfall regimes in the study area.

In the Mediterranean, long dry periods alternate with short but intense seasonal rainfalls [START_REF] Miller | Canopy structure of mediterranean-type shrubs in relation to heat and moisture[END_REF][START_REF] Cowling | Rainfall reliability, a neglected factor in explaining convergence and divergence of plant traits in fire-prone mediterranean-climate ecosystems[END_REF]. During the summer, contaminants from human activities tend to concentrate in soils and in river sediments as water resources decrease through evaporation [START_REF] Gasith | Streams in Mediterranean climate regions: Abiotic influences and biotic responses to predictable seasonal events[END_REF][START_REF] Olías | Seasonal water quality variations in a river affected by acid mine drainage: The Odiel River (South West Spain)[END_REF]). During heavy rain events in the Mediterranean, flow augmentation, surface soil runoff, river sediments and the resuspension of sewer solids from Combined Sewer Overflows (CSOs) lead to a modification in the level of contaminants in surface waters [START_REF] Coynel | Scenario of particulate trace metal and metalloid transport during a major flood event inferred from transient geochemical signals[END_REF][START_REF] Sakrabani | Biodegradability of organic matter associated with sewer sediments during first flush[END_REF][START_REF] Taghavi | Flood event impact on pesticide transfer in a small agricultural catchment (Montoussé at Auradé, south west France)[END_REF][START_REF] Dumas | Riverine transfer of anthropogenic and natural trace metals to the Gulf of Lions (NW Mediterranean Sea)[END_REF] and in particular, an increase in the level of pharmaceutical products [START_REF] Escher | Environmental toxicology and risk assessment of pharmaceuticals from hospital wastewater[END_REF]Osorio et al. 2012a;[START_REF] Backhaus | Medicines, shaken and stirred: a critical review on the ecotoxicology of pharmaceutical mixtures[END_REF][START_REF] Lalone | Leveraging existing data for prioritization of the ecological risks of human and veterinary pharmaceuticals to aquatic organisms[END_REF][START_REF] Osorio | Hydrological variation modulates pharmaceutical levels and biofilm responses in a Mediterranean river[END_REF]. Drought/flood hydrodynamics of coastal Mediterranean rivers have been identified as the major factor affecting the presence of micropollutants in these environments [START_REF] Pailler | Behaviour and Fluxes of Dissolved Antibiotics, Analgesics and Hormones During Flood Events in a Small Heterogeneous Catchment in the Grand Duchy of Luxembourg[END_REF]Osorio et al. 2012a;[START_REF] Fairbairn | Sources and transport of contaminants of emerging concern: A two-year study of occurrence and spatiotemporal variation in a mixed land use watershed[END_REF].

While [START_REF] Ashley | The erosion and movement of sediments and associated pollutants in combined sewers[END_REF] predicted the release of high concentrations of contaminants during first foul flushes from CSOs at the start of a wet weather flow, studies of pollutant dynamics during heavy rains and high flows are scarce in the literature [START_REF] Pailler | Behaviour and Fluxes of Dissolved Antibiotics, Analgesics and Hormones During Flood Events in a Small Heterogeneous Catchment in the Grand Duchy of Luxembourg[END_REF]. [START_REF] Reoyo-Prats | Multicontamination phenomena occur more often than expected in Mediterranean coastal watercourses: Study case of the Têt River (France)[END_REF] performed a high frequency follow-up of six families of contaminants, including pesticides, nutrients, heavy metals and fecal indicators, during a flood event. They demonstrated that floods and perhaps even more importantly, recurrent heavy Mediterranean rain events, are responsible for chronic multicontamination phenomena under this climate. Such micropollutant mixtures are potentially highly toxic for living organisms of the land-sea continuum, even at small concentrations [START_REF] Schwarzenbach | The Challenge of Micropollutants in Aquatic Systems[END_REF]. Determining the importance of PACs compared to other pollutants is therefore fundamental to assess their environmental risk [START_REF] Boxall | Review Pharmaceuticals and Personal Care Products in the Environment : What Are the Big Questions ?[END_REF], particularly during first flush multicontamination phenomena. Moreover, few environmental studies detail the fate and the way pharmaceuticals are transported along well-characterized river watersheds during different seasons [START_REF] Pal | Impacts of emerging organic contaminants on freshwater resources: Review of recent occurrences, sources, fate and effects[END_REF], particularly in conjunction with other micropollutants [START_REF] Biales | Complex watersheds, collaborative teams: Assessing pollutant presence and effects in the San Francisco Delta[END_REF][START_REF] Fairbairn | Sources and transport of contaminants of emerging concern: A two-year study of occurrence and spatiotemporal variation in a mixed land use watershed[END_REF][START_REF] Garrido | Monitoring of emerging pollutants in Guadiamar River basin (South of Spain): analytical method , spatial distribution and environmental risk assessment[END_REF]) and during storm events and floods [START_REF] Pailler | Behaviour and Fluxes of Dissolved Antibiotics, Analgesics and Hormones During Flood Events in a Small Heterogeneous Catchment in the Grand Duchy of Luxembourg[END_REF]). The lack of information regarding their sources, seasonal occurrence and potential environmental risk might be the reason why medical and veterinary drugs are not yet properly regulated and thoroughly included in monitoring lists of priority pollutants.

In this study, our first goal was to finely assess the dynamics of PACs concentrations during first foul flushes from CSOs from heavy rainfalls on the Têt River, a coastal Mediterranean watercourse situated in the Southeast of France. This river was previously used as a model to study the transfer of suspended matter, associated nutrients and pollutants to the Mediterranean Sea [START_REF] Roussiez | Background levels of heavy metals in surficial sediments of the Gulf of Lions (NW Mediterranean): An approach based on 133Cs normalization and lead isotope measurements[END_REF][START_REF] Garcia-Esteves | Predicting the impact of land use on the major element and nutrient fluxes in coastal Mediterranean rivers: The case of the Têt River (Southern France)[END_REF][START_REF] Dumas | Riverine transfer of anthropogenic and natural trace metals to the Gulf of Lions (NW Mediterranean Sea)[END_REF][START_REF] Reoyo-Prats | Multicontamination phenomena occur more often than expected in Mediterranean coastal watercourses: Study case of the Têt River (France)[END_REF]. Our second goal was to identify the different sources of PACs and their contribution to river water contamination. Because no animal farming activities exist along the Têt River, all drugs are derived from human consumption and their sources are thus easier to identify. Our final goal, though less significant that the others, was to determine how concentrations of PACs changed relative to a variety of other contaminants, particularly during typical Mediterranean heavy rainfall events.

Material and Methods

Sampling

Stations

To follow the dynamics of medical residues, we focused our efforts on the most urbanized section of the Têt River, from Perpignan, the biggest city in the Pyrénées-Orientales department, to the river mouth (Fig. 1). Stations were chosen according to two criteria: for safety reasons, sampling had to be done from bridges during flood risk events and sampling locations had to be close enough to the laboratory to be sampled at high frequency. The first station, R2, was located in the city of Perpignan, downstream the Basse River (Conseil Général des Pyrénées Orientales (CG66) 2009; Conseil Général des Pyrénées Orientales (CG66) 2012), a tributary highly polluted by agriculture and urban activities and downstream Perpignan city's right bank sewage-overflow system, but upstream main WWTP inputs (Fig. 1). The second station, R0, was located 9 km downstream from R2 at the first bridge downstream from all Perpignan's sewage-overflow system and WWTP (350,000 Population Equivalent, PE) and the WWTPs for the Sainte-Marie-la-Mer and Canet-en-Roussillon vacation resorts (24,000 PE and 66,000 PE, respectively). Both R0 and R2 received inputs from extensive gardening and vineyard culturing activities from the last portion of the catchment area (Fig. 1, main figure and inlet).

Campaigns

A campaign comprised three sampling periods including one autumnal storm event with high frequency sampling, and two punctual samplings on either side of this event, during summer and winter droughts.

The first campaign was a pilot study that aimed to describe pharmaceutical trends for a limited number of molecules to determine if a more thorough study would be relevant. Sampling was conducted at the R0 station of the Têt River during periods of drought in summer 2013 (17 th of September) and winter 2014 (13 th of February). Eleven additional samples were collected from the 16 th to the 21 st of November 2013 during an autumnal flood event, with a total follow-up time of 109 h. This campaign also compared the behavior of PACs with that of other contaminants' reported in [START_REF] Reoyo-Prats | Multicontamination phenomena occur more often than expected in Mediterranean coastal watercourses: Study case of the Têt River (France)[END_REF] during the same flood event.

Based on the results from the first campaign, a second campaign was carried out in 2015-16 to follow the temporal dynamics of a larger number of drugs, upstream and downstream WWTPs, at R2 and R0 respectively (Fig. 1). Punctual samplings during periods of drought were conducted in summer 2015 (18 th of August) and winter 2016 (10 th of February). In addition, a total of 10 samplings were done from the 2 nd to the 4 th of November 2015 at R2 and R0 stations, with a total of 5 samplings per station and a total follow-up of 53 h.

Heavy rains and flood sampling strategy

Perpignan city's sewer system is designed to overflow according to the 2-month return period as defined by the maximum precipitation distribution curves for the Têt River (Online resource 1). Therefore, in order to determine the level of contaminants before and after the Perpignan sewer system overflowed, sampling began when precipitation approached the 2-month return period distribution curve and continued at a high frequency thereafter. Sampling time frequency during first heavy rains was controlled by the "live" monitoring of precipitation on the Perpignan hydrological station's webpage, which is located within the Ministère de la Transition Écologique et Solidaire website (http://www.rdbrmc.com/hydroreel2/station.php?codestation=459).

Once we estimated that the sewer system had overflowed, sampling frequency was determined by either changes in precipitation or river water flow levels, which are available at this same webpage.

Reported final river flow was retrieved at the Banque Hydro (http://www.hydro.eaufrance.fr).

Medical, veterinary and emerging drugs analyses

River water was sampled with 10 L capacity polyethylene terephthalate (PET) tanks previously cleaned with 1.5 L of 1 M HCl and rinsed with 2 L of distilled water. Plastic tanks were required for field sampling during the heavy rainfalls and for the flood, as crude weather conditions did not allow for sampling with independent bottles for each contaminant. In the field, tanks were rinsed three times with river water before they were filled. Once in the laboratory, 0.8 L of sampled water was transferred to an uncontaminated 1 L bottle and put at -20°C for future pharmaceutical analyses. The use of plastic bottles instead of glass bottles was necessary for safety reasons as glass containers explode when stored at freezing temperatures. This constraint requires that results are interpreted with caution, particularly with respect to comparative studies of absolute concentrations of molecules.

PACs analyses from the 2013-14 pilot campaign were conducted in the dissolved fraction using internally validated methods. The pharmaceuticals chosen for analyses were four medical residues commonly found in waters: sulfamethoxazole, carbamazepine, ibuprofen and diclofenac.

The Solid Phase Extraction (SPE) was performed using PolyCleanTM 2H (60 µm, 70 Å), 200 mg / 6 mL (Interchim, Montluçon, France) and a vacuum manifold from phenomenex (Torrance, CA, USA). After cartridge conditioning with 20 mL of methanol followed by 6 mL of reagent water, samples were loaded onto the cartridge at a flow rate of 5-10 mL/min. The cartridge was then dried under vacuum for approximately 5 min. The cartridge was eluted with 5 mL of methanol and the extract was filtered through a 0.45 µm PolyTetraFluoroEthylene (PTFE) filter (VWR, Radnor, PA, USA), concentrated until dry under a gentle stream of nitrogen in a dry bath at 50 ± 5 °C and reconstituted in a volume of 200 µL of methanol. HPLC analyses were performed by Liquid Chromatography-Mass Spectrometry (LC-MS). The chromatography was done in a Waters Alliance 2695 separation module (Waters, Milford, MA, USA) equipped with a reversed phase C18 analytical column of 100 mm × 2.2 mm and 3.5 µm particle size (XTerra MS C18). The mobile phases used were ultrapure water and acetonitrile at different percentages (Table 1). The injected sample volume and flow-rate was 5 µL and 0.25 mL/min respectively. The mass spectral analysis was performed on a Waters micromass ZQ (Waters, Milford, MA, USA) equipped with a Z-Spray™ ElectroSpray Ionization (ESI) source (Table 2). Detection and quantification limits are also reported in Table 2.

During the 2015-16 campaign, a more thorough monitoring of pharmaceutical products was conducted on a total of 89 molecules. For this purpose, the total fraction was analyzed by COFRAC accredited laboratory Groupe CARSO, via Liquid Chromatography-tandem Mass Spectrometry (LC-MS/MS) using a Thermo Scientific TSQ Quantum (Waltham, MA, USA) connected to an API 5500 mass spectrometer (Applied Biosystems-Sciex, Foster City, CA, USA) or an Agilent 6460 (Santa Clara, CA, USA). Ionization was performed in the switch positive/negative mode for all runs. Chromatographic analyses were achieved using a C18 analytical column of 100 mm × 2.2 mm and 3.5 µm particle size (XTerra MS C18) with a non-polar elution gradient at a flow rate of 0.4 mL/min while the injected sample volume varied from 30 to 5,000 µL depending on the method. Direct injection was performed for diclofenac. Direct injection and EDTA treatment were used for tetracycline, chlortetracycline, doxycycline, oxytetracycline, ciprofloxacin, danofloxacin, norfloxacin, ofloxacin, enrofloxacin, marbofloxacin and sulfamethoxazole. The same instrument was used with on-line SPE injection with methanol on a C18 cartridge for trimethoprim, lincomycin, clindamycin, chloramphenicol, clenbuterol, erythromycin, clarithromycin, spiramycin, monensin, virginiamycin M1 (pristinamycin IIA), roxithromycin, irbesartan, indomethacin, fenoprofen, antipyrine (phenazone), propyphenazone, dexamethasone, naproxen, ibuprofen, ibuprofen 1-hydroxy, ibuprofen 2-hydroxy, piroxicam, ketorolac (Macril), hydrocortisone, fenofibric acid, clofibric acid, simvastatine, aminopyrine, tramadol, desmethyltramadol, dihydrocodeine, oxycodone, morphine, cyclophosphamide, methotrexate, ifosfamide, prilocaine, mepivacaine, bupivacaine, atenolol, propranolol, acebutolol, metoprolol, bisoprolol, sotalol, timolol, diltiazem, ticlopidine, pentoxifylline, fluoxetine, doxepine, imipramine, amitriptyline, zolpidem, norfluoxetine, lorazepam, oxazepam, bromazepam, diazepam, alprazolam, carbamazepine, phenytoin (Dilantin), caffeine, salbutamol, progesterone, testosterone, 4-androstenedione, norethindrone, ranitidine, cotinine and nicotine. Again, the same instrument was used but with off-line SPE injection with methanol on an ion exchange cartridge according to ISO 25101 (International Organization for Standardization 2009) for estrone, ethinylestradiol, estradiol 17 alpha, estradiol 17 beta, estriol, dienestrol and diethylstilbestrol. Detection limit and the method used for different analytes analysis are summarized in Table 3.

Mass fluxes (g/h) were calculated for the 2015-16 campaign as mass concentrations (g/L) multiplied by river flow (L/h). Flow data were only available at R2 station but as R0 station was only 9 km downstream from R2 and is at the very end of the catchment area, we assumed the same flow for both.

In order to better discern the relationship of PACs with positive concentrations at R0 station during the rain event of the 2015-16 campaign, a multivariate analyses using Principal Component Analysis (PCA) was conducted using rda function from vegan package [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2[END_REF]) in R freeware (R Core Team 2017).

Rare Earth Element (REE) analyses

REE concentrations were measured by Inductively Coupled Plasma Mass Spectrometry (ICP-MS, Agilent 7700X) at CEFREM laboratory. REE concentrations followed an intercalibration exercise conducted on river water certified reference material SLRS-6 developed for the analysis of trace metals and other constituents [START_REF] Yeghicheyan | A New Compilation of Element Concentrations in the Natural River Water Standard SLRS-6 (NRC-CNRC)[END_REF]. This exercise demonstrated that REE concentrations provided by CEFREM were consistent with other participants. REE patterns were determined by using the Upper Continental Crust (UCC) elemental composition of rock standards in order to normalize sample concentrations of these elements [START_REF] Piper | Normalized Rare Earth Elements in Water, Sediments, and Wine: Identifying Sources and Environmental Redox Conditions[END_REF].

Results

Pilot study on pharmaceutical trends in the Têt River during the 2013-14 sampling campaign

Our qualitative study during the 2013-14 campaign sought to estimate major trends of dissolved pharmaceuticals for the R0 station (Fig. 1) located in the Têt River downstream of Perpignan, Canet and Sainte-Marie-la-Mer WWTPs, during a typical autumnal flood event as compared to summer and winter droughts. Four dissolved pharmaceuticals (carbamazepine, sulfamethoxazole, ibuprofen and diclofenac) were studied as representatives of three families of pharmaceuticals (anticonvulsants, antibiotics and anti-inflammatories, respectively) that are frequently found in river waters. Pharmaceutical levels for all seasons, in parallel with changes in the river flow during the storm event, are summarized in figure 2. Maximum levels of dissolved PACs, with the exception of ibuprofen which was not detected, were attained at 20 h sampling time.

Exhaustive monitoring of medical, veterinary and emerging drugs on the total fraction during the 2015-16 sampling campaign over several seasons

Given the encouraging results of our pilot study, a more extensive study was conducted during the 2015-16 campaign. A total of 89 molecules were analyzed in the total fraction at both R0 and R2 sampling stations located upstream and downstream of the primary WWTPs in urbanized areas of the Têt River (Fig. 1). Twenty-one molecules were identified and their concentrations are plotted together with river flow for the 2015 autumnal storm event, and for the summer and winter droughts at both stations, R2 (Fig. 3a) and R0 (Fig. 3b). These results are also summarized in the Online Resource 2. In parallel, figures 3c and d show discharges of the main Sewer-Overflows (SO) for the city of Perpignan, river flow and precipitations along sampling time during the same storm event.

All molecules detected in 2013-14 (Fig. 2) were also found in 2015-16, except for sulfamethoxazole, detected only in 2013-14 with a peak concentration of 4.4 ng/L, well beneath the 2015-16 quantification limit (20 ng/L, Table 3). For both stations, the highest drug concentrations were found for the winter drought (right axes in Fig. 3a&b) and at the start of sampling during the autumnal storm event (0 h sampling time in Fig. 3a&b). At R2, seven molecules were detected but nicotine and caffeine were dominant at the beginning of the autumnal storm event (Fig. 3a), with concentrations of 217 ng/L and 175 ng/L, respectively (Online Resource 2). These values were at least 2.5 times higher than those of other molecules for all seasons for this station. At R0, however, a total of 20 molecules were detected, corresponding to 11 PACs families (Fig. 3b and Online Resource 2). Top concentrations for nicotine and caffeine were less important with 149 ng/L and 59 ng/L, respectively.

But analgesics, represented by tramadol and desmethyltramadol, had the highest cumulated concentrations with 382 ng/L at 0 h and 289 ng/L at 18 h sampling time for the storm event and with 297 ng/L and 224 ng/L for the winter and summer droughts respectively, with tramadol as the most abundant molecule in this family (Online Resource 2). Anti-inflammatories (ibuprofen and diclofenac)

were also dominant in our study with 93 ng/L and 293 ng/L cumulated concentrations at 0 h and 18 h respectively and 171 ng/L in winter while only 30 ng/L in summer. Anti-hypertensives, only represented by irbesartan, also had important concentrations, 195 ng/L at 0 h and 85 ng/L at 18h during the autumnal storm event, and 182 ng/L and 132 ng/L in winter and summer droughts respectively. With a total of five (see Fig. 3), beta blockers were the most significant family in terms of the number of molecules found in the river, with cumulated concentrations of 146 ng/L at 0 h, 130 ng/L at 18 h, 112 ng/L in the summer drought and 144 ng/L in the winter drought.

Pharmaceuticals were also monitored at a high frequency during the autumnal heavy rain event. At R2, all six molecules detected at 0 h sampling time decreased remarkably after this first sampling and were at their lowest concentration at the flow peak (144 m 3 /s) that occurred at around 23 h sampling, except for irbesartan, which increased slightly at this point (Fig. 3a andOnline Resource 2). In contrast to R2, at R0, most PACs remained at high concentration during the first rains until 18 h (Fig. 3b and Online Resource 2), when the major overflows of the sewer systems occurred (Fig. 3c&d). They then decreased sharply and found their lowest levels around the flow peak.

Exceptions to this rule existed for irbesartan which showed a linear dramatic decrease similar to emerging drugs in R2, and ibuprofen and fenofibric acid whose concentrations peaked sharply at 18 h and then decreased rapidly at the flow peak. A principal component analysis to explore the relationships among PACs concentrations throughout the rain event in the R0 station was also conducted. The ordination plot showed that all PACs were negatively correlated with respect to the first axis and therefore to initial moments of the rain event, at 0 h and 18 h sampling times (PCA, Online Resource 3). The second axis discriminated amongst PACs correlated to 0 h and those correlated to 18 h. While fluoxetine and irbersartan were very positively correlated to axis 2 and had indeed the highest concentrations at 0 h sampling time (see also Online Resource 2), fenofibric acid and ibuprofen were very negatively correlated to axis 2 and had their highest concentrations at 18 h sampling time.

In order to consider the influence of the river flow, we calculated instantaneous mass fluxes at both stations, R2 (Fig. 4a) and R0 (Fig. 4b). Considering seasonal variations, most PACs fluxes were higher in winter than in summer at both R2 and R0 stations, although some were 0 g/L for both seasons. Exceptions to this rule were nicotine at R2, and acebutolol at R0 with the opposite behavior.

During the autumnal heavy rain event at R2 station, cotinine and bisoprolol fluxes were always around 0 g/h and so was the case for caffeine except from a flux of 2.1 g/h at the very beginning of the event (Fig. 4a). Tramadol and irbesartan fluxes were also near 0 g/h throughout the rain event with the exception of peaks of 1.8 and 4.1 g/h, respectively, but at different times. While tramadol peaked at 18 h, irbesatan did so at 23 h. Nicotine had two high peak fluxes, one of 14.8 g/h at 18 h, followed by a dramatic decrease to 0 g/h at 23 h and a second peak at 11.2 g/h at 29 h before a soft decrease to 5.8 g/h at 51 h. At R0 station, all top flux peaks were attained at 22 h with the highest fluxes of 30.6 g/h for caffeine, 21.3 g/h for diclofenac and nicotine, 19.2 g/h for ibersartan and 18.7 g/h for tramadol. All these four molecules decreased sharply after 22 h but never reached 0 g/h with the exception of diclofenac. All other molecules that peaked at 22 h showed a rapid decrease to 0 g/h after 22 h. Finally, ibuprofen, sotalol, clarithromycin, atenolol, acebutolol, diltiazem, clarithromycin, erythromycin and trimethoprim had their peak fluxes at 18 h followed by a sharp decrease immediately after, and remained at near 0 g/h during the rest of the rain event.

Rare earth element patterns during a heavy rain in autumn 2015

Rare Earth Element (REE) patterns were defined at R2 and R0 stations to determine the contribution of hospitals as sources of PACs in the Têt River. Upstream station R2 showed a generally flat pattern of REEs along different sampling times, except for europium (Eu) with small peaks at 0 and 18 h sampling times (Fig. 5a), probably due to the mineral origin of this element in the Têt River catchment area. However, the downstream station R0 showed mostly large peaks of gadolinium (Gd) during the same sampling times (Fig. 5b). Indeed, at 0 and 18 h, Gd values were respectively 6 and 18 times higher at R0 than at R2.

Discussion

Dissolved PACs contribute to chronic risk multicontamination phenomena for coastal

Mediterranean rivers

In 2013-14, we conducted an exploratory qualitative study whose main purpose was to compare pharmaceutical dynamics with those of other pollutants during heavy rainfalls and floods in the Têt River, a typical coastal Mediterranean watercourse. We studied the dissolved fraction because it allows for comparisons to be made with the only other known published study on PACs during storm events and sewer overflows [START_REF] Pailler | Behaviour and Fluxes of Dissolved Antibiotics, Analgesics and Hormones During Flood Events in a Small Heterogeneous Catchment in the Grand Duchy of Luxembourg[END_REF]). Although our results were only internally validated and should be interpreted with caution, they showed a simultaneous increase of dissolved pharmaceuticals just before the increase in river flow during the first heavy rains at 20 h sampling time (Fig. 2). [START_REF] Pailler | Behaviour and Fluxes of Dissolved Antibiotics, Analgesics and Hormones During Flood Events in a Small Heterogeneous Catchment in the Grand Duchy of Luxembourg[END_REF] [START_REF] Gromaire | Contribution of different sources to the pollution of wet weather flow in combined sewers[END_REF][START_REF] Mclellan | Distribution and Fate of Escherichia coli in Lake Michigan Following Contamination with Urban Stormwater and Combined Sewer Overflows[END_REF][START_REF] Phillips | Wastewater effluent, combined sewer overflows, and other sources of organic compounds to Lake Champlain[END_REF][START_REF] Weyrauch | Contribution of combined sewer overflows to trace contaminant loads in urban streams[END_REF][START_REF] Passerat | Impact of an intense combined sewer overflow event on the microbiological water quality of the Seine River[END_REF]. Because heavy rainfalls are recurrent under a Mediterranean climate, Reoyo-Prats et al. 2017 argued that first flush multicontamination phenomena constitute a chronic risk in Mediterranean coastal environments and, according to our results, dissolved PACs are no exception.

A fine-scale study of PACs dynamics during the 2015 rain event allowed for identification of sources, fate and transport of PACs in surface waters

A higher number and concentration of pharmaceutical compounds in the total fraction was found downstream (R0 station, Fig. 3b) than upstream WWTPs (R2, Fig. 3a) at all seasons, indicating their constant release by WWTPs, which are recognized as major sources of pharmaceuticals in the environment [START_REF] Santos | Occurrence and risk assessment of pharmaceutically active compounds in wastewater treatment plants. A case study: Seville city (Spain)[END_REF][START_REF] Fernández | Occurrence of pharmaceutically active compounds in surface waters of the henares-jarama-tajo river system (madrid, spain) and a potential risk characterization[END_REF][START_REF] López-Roldán | Advanced monitoring of pharmaceuticals and estrogens in the Llobregat River basin (Spain) by liquid chromatography-triple quadrupoletandem mass spectrometry in combination with ultra performance liquid chromatography-time of flight-mass spectrometry[END_REF][START_REF] Vazquez-Roig | Assessment of the occurrence and distribution of pharmaceuticals in a Mediterranean wetland (L'Albufera, Valencia, Spain) by LC-MS/MS[END_REF][START_REF] Barber | Persistence and Potential Effects of Complex Organic Contaminant Mixtures in Wastewater-Impacted Streams[END_REF]). The highest concentrations of PACs in the Têt River were found at the beginning of the autumnal heavy rain event (0 h sampling time on Fig. 3) and during the winter drought. This result could be explained by the higher river flow during these seasons leading to a higher resuspension of suspended matter, in contrast to the summer drought. Other studies of Mediterranean rivers also reported higher PACs levels during this season (Osorio et al. 2012a;Osorio et al. 2012b). Another possible explanation is the decrease in the attenuation of PACs during winter and autumn due to the decline of phototransformation processes at these seasons [START_REF] Daneshvar | Seasonal variations in the occurrence and fate of basic and neutral pharmaceuticals in a Swedish river-lake system[END_REF], which largely contribute to PACs elimination in WWTPs [START_REF] Vieno | Seasonal variation in the occurrence of pharmaceuticals in effluents from a sewage treatment plant and in the recipient water[END_REF]. However, in order to compare seasonal variations, it is best to use mass fluxes of organic pollutants such that river flow can be accounted for [START_REF] Bernot | Human and veterinary pharmaceutical abundance and transport in a rural central Indiana stream influenced by confined animal feeding operations (CAFOs)[END_REF][START_REF] Moreno-González | Seasonal input of regulated and emerging organic pollutants through surface watercourses to a Mediterranean coastal lagoon[END_REF]. Our flux results indicated that for most PACs, levels were slightly higher in winter than in summer. This result differs from the study of [START_REF] Bernot | Human and veterinary pharmaceutical abundance and transport in a rural central Indiana stream influenced by confined animal feeding operations (CAFOs)[END_REF], which only found caffeine and paraxanthine more aboundant in winter whereas the abundance of some personal care products actually increased in summer because they were consumed at this time (i.e. DEET from insect repellents). This result is not captured in Fig. 4, because summer and winter drought fluxes were negligible in comparison to the high fluxes of PACs observed for the rain event. To the best of our knowledge no study on total PACs mass flux loads during a rain event has ever been conducted. Due to the high recurrence of heavy rains under a Mediterranean climate, these events largely contribute to PACs release into aquatic environments.

The most abundant PACs families encountered at R0 were analgesics, anti-inflammatories and beta blockers (Fig. 2), which are known to be the most resistant to wastewater treatments with a removal rate of only 30-40% (reviewed in [START_REF] Deblonde | Emerging pollutants in wastewater: A review of the literature[END_REF]. Analgesics had the highest cumulated concentrations. In France, analgesics are the most popular prescription drugs with tramadol, a drug found in high concentrations in the Têt River, as the fourth most sold analgesic (Agence Nationnale de la Sécurité du Médicament et des produits de santé (ANSM) 2014). Analgesics are also found in high abundance in other freshwater systems in Europe [START_REF] Kasprzyk-Hordern | The occurrence of pharmaceuticals , personal care products , endocrine disruptors and illicit drugs in surface water in South Wales[END_REF][START_REF] Rúa-Gómez | Upgrading of wastewater treatment plants through the use of unconventional treatment technologies: Removal of lidocaine, tramadol, venlafaxine and their metabolites[END_REF][START_REF] Loos | EU-wide monitoring survey on emerging polar organic contaminants in wastewater treatment plant effluents[END_REF][START_REF] Osorio | Concentration and risk of pharmaceuticals in freshwater systems are related to the population density and the livestock units in Iberian Rivers[END_REF]. Anti-inflammatories were also predominant in our study, a result similar to other Mediterranean European countries such as Italy [START_REF] Patrolecco | Simultaneous determination of human pharmaceuticals in water samples by solid phase extraction and HPLC with UV-fluorescence detection[END_REF] and Spain [START_REF] Moreno-González | Input of pharmaceuticals through coastal surface watercourses into a Mediterranean lagoon (Mar Menor, SE Spain): Sources and seasonal variations[END_REF][START_REF] Osorio | Concentration and risk of pharmaceuticals in freshwater systems are related to the population density and the livestock units in Iberian Rivers[END_REF], as well as in other developed countries around the world [START_REF] Stumpf | Polar drug residues in sewage and natural waters in the state of Rio de Janeiro, Brazil[END_REF][START_REF] Tixier | Occurrence and Fate of Carbamazepine, Clofibric Acid, Diclofenac, Ibuprofen, Ketoprofen and Naproxen in Surface Waters[END_REF][START_REF] Bendz | Occurrence and fate of pharmaceutically active compounds in the environment, a case study: Höje River in Sweden[END_REF][START_REF] Roberts | The occurrence of selected pharmaceuticals in wastewater effluent and surface waters of the lower Tyne catchment[END_REF][START_REF] Palmer | Sources of pharmaceutical pollution in the New York City Watershed[END_REF]. Beta blocker levels were also significant in the Têt River. Among them, sotalol was predominant at our study site as well as at another nearby coastal river (Osorio et al. 2012a), although atenolol, metoprolol and propranolol are typically the most common beta blockers found in aquatic environments (revised in [START_REF] Godoy | Occurrence, ecotoxicological effects and risk assessment of antihypertensive pharmaceutical residues in the aquatic environment -A review[END_REF].

At R2 station, however, only nicotine and caffeine were found in abundance at the beginning of the storm event. As they are common urban contaminants coming from human consumption of cigarettes and coffee, their presence was either due to discharges from WWTPs from small towns situated upstream R2, or to illicit sanitary inputs, i.e. houses not connected to the wastewater collection system that directly and illegally release their untreated wastewaters in rivers.

Illegal sanitary inputs were reported to us by locals at the Basse River tributary, which is directly upstream R2. The hypothesis of a nearby discharge of these drugs by illicit sanitary inputs is supported by the presence of nicotine and caffeine at higher concentrations at R2 compared to R0, where these molecules must be buffered by the WWTPs. If this hypothesis is confirmed, the rapid release of these emerging drugs during the first hours of a heavy rain could therefore serve as markers of illicit sanitary inputs.

Regarding PACs dynamics during the storm event at 2015-16 campaign, only two samplings were done at the beginning of the rain event, at 0 and 18 h respectively, but we are confident that our second sampling point is fully representative of the CSO phenomena and no major changes in pollutants occurred in between for two major reasons. Firstly, although a CSO occurred at 10 h, it was minor, with ~2,000 m 3 of untreated water released into the river, while the major CSO at 16 h represented a volume of approximately 7,000 m 3 . Secondly, a higher frequency follow-up during the flood event of the 2013-14 campaign included a sampling point during first CSO of approximately 18,000 m 3 , which did not lead to other contaminants increase. Only at the second CSO (~50,000 m 3 ) pollutants peaked [START_REF] Reoyo-Prats | Multicontamination phenomena occur more often than expected in Mediterranean coastal watercourses: Study case of the Têt River (France)[END_REF] including dissolved PACs (see above). Coming back to the 2015 storm event, at R0 station the great majority of PACs remained at a constant concentration during first flush events (18 h sampling time in Fig. 3b), resembling their concentration in winter (vertical right axis in Fig. 3b). This result can be explained by CSO during first flushes (Fig. 1 and Fig. 3d). However, right bank sewer overflows SO11, SO64 and SO50 situated upstream R2 (Fig. 1) also overflowed during this flood (Fig. 3c), but concentrations of medical drugs decreased rapidly at this station (Fig. 3a) as is usually expected because of the rapid dilution of molecules with increased water flow during first rains. First, the absence of pharmaceuticals other than caffeine and nicotine at 0 h sampling time at R2 indicate that, in spite of numerous agricultural areas around our sampling stations (Fig. 1), diffuse sources of PACs via runoff from potential sewage sludge used as fertilizers did not occur during intense rainfalls. Our study therefore supports previous findings that diffuse sources are less important than point sources for pharmaceutical products in surface waters [START_REF] López-Roldán | Advanced monitoring of pharmaceuticals and estrogens in the Llobregat River basin (Spain) by liquid chromatography-triple quadrupoletandem mass spectrometry in combination with ultra performance liquid chromatography-time of flight-mass spectrometry[END_REF][START_REF] Vazquez-Roig | Assessment of the occurrence and distribution of pharmaceuticals in a Mediterranean wetland (L'Albufera, Valencia, Spain) by LC-MS/MS[END_REF][START_REF] Barber | Persistence and Potential Effects of Complex Organic Contaminant Mixtures in Wastewater-Impacted Streams[END_REF]. Second, the prevalence of prescribed drugs at first flushes at R0 but not at R2 demonstrate that major sources of prescribed medical drugs along the Têt River are situated on the left bank (Fig. 1). The public hospital of Perpignan city lies on this bank and is, with 1,183 bed-places, the largest hospital in the Pyrénées-Orientales department and the third of its category in the Languedoc-Roussillon region.

To determine if hospitals are the source of pharmaceuticals, we conducted an analysis of the trace metal gadolinium (Gd). This element is used as a contrast agent in Magnetic Resonance Imaging (MRI) and is not removed by sewage water treatments, and thus produces positive anthropogenic anomalies in Rare Earth Element (REE) patterns in rivers [START_REF] Rabiet | Positive gadolinium anomalies in wastewater treatment plant effluents and aquatic environment in the Hérault watershed (South France)[END_REF][START_REF] Kulaksiz | Anthropogenic gadolinium as a microcontaminant in tap water used as drinking water in urban areas and megacities[END_REF][START_REF] Piper | Normalized Rare Earth Elements in Water, Sediments, and Wine: Identifying Sources and Environmental Redox Conditions[END_REF]. A general flat pattern of REEs was observed at R2 station (Fig. 5a), which received discharges of only right bank SOs, i.e. SO11, SO64 and SO50 (Fig. 1). In contrast, R0 station, which also received left bank discharges from SO09 and SO10 (Fig. 1), showed sharper peaks of Gd REE element in autumn (Fig. 5b). Both right and left bank SOs collect urban water, but left bank SOs discharged vast quantities of pharmaceutical compounds compared to the right bank.

Worth noting is that the left bank collects public hospital wastewaters. These results indicate that the public hospital is responsible for the origin of the high concentration of medical drugs detected in the Têt River at R0 station.

When looking at fluxes at R0 station (Fig. 4.b), three different behaviors were observed. Firstly, ibuprofen, sotalol, clarithromycin, atenolol, acebutolol, diltiazem, clarithromycin, erythromycin and trimethoprim had their peak fluxes at 18 h during the CSOs while they remained near 0 g/h during the rest of the rain event. Therefore, these compounds seemed to come mostly from CSOs while under normal circumstances they are eliminated by WWTPs. Secondly, diclofenac, desmethyltramadol, oxazepam, fenofibric acid, carbamazepine propranolol and bisoprolol reached their highest fluxes during the flow peak at 22 h but their fluxes were nearly as high during the CSOs at 18 h while they decreased sharply to 0 g/h just after the flow peak at 27 h. These second group of compounds were therefore released through direct discharges from CSOs as well as by the remobilization of the riverbed sediments. Thirdly, caffeine, nicotine, ibersartan and tramadol increased during the CSOs at 18 h before they attained their highest values (top levels among all PACs) during the flow peak at 22 h, then they sharply decreased until 27 h before they decreased softly until the end of the rain event, never attaining a null flux value. This dynamic seemed to characterize compounds predominantly remobilized from river sediments at the flow peak, but also those punctually discharged by CSOs and constantly released into the river by WWTPs. Interestingly, these compounds are those that are predominantly found at R2 station (Fig. 4.b) and are not released by the hospital but are entirely derived from urban discharges. At the flow peak in R2, all compounds had a null flux except for irbesartan. This proves that river sediments in R2 are not as heavily contaminated as those in R0, which is the station downstream main urban WWTPs.

Irbesartan's unique behavior can be explained by its physicochemical characteristics that are discussed hereafter.

As the 2015-16 pharmaceutical campaign was conducted on the total fraction, corresponding to both the dissolved and the suspended water fractions, hydrophilic and hydrophobic properties of PACs had an important impact on their fate and transport. While at R0 most molecules remained at a constant concentration during CSOs first flush event at 18 h (Fig. 3b), several of them behaved differently, especially fluoxetine, irbesartan, fenofibric acid and ibuprofen, which were also more highly correlated to axis 2 in the PCA ordination plot (Online resource 3). Fluoxetine had a null concentration throughout the rain event, and before the event it had a concentration of 1 ng/L (Online resource 2), which explains its particular position in the PCA (Online resource 3). Irbesartan was also highly positively correlated with axes 2 and its highest concentration occurred at the beginning of the rain event, and then rapidly decreased (Fig. 3.b). This rapid dilution behavior is typical of hydrophobic molecules that are primarily attached to suspended matter. Irbesartan has a logarithm of octanol-water partition coefficient (log K ow ) of 5.31 (Table 4) and a very high organic carbon-based sorption coefficient (K oc = 1,300,000), and together these properties confirmed its hydrophobic nature. A log K ow higher than 5 and a high molecular weight (irbesartan weights 428.53 g/mol) are properties typical of molecules easily sorbed to sediments according to [START_REF] Pal | Impacts of emerging organic contaminants on freshwater resources: Review of recent occurrences, sources, fate and effects[END_REF] and sorption has been identified as one of the predominant attenuation mechanisms for PACs [START_REF] Riml | Evaluating the fate of six common pharmaceuticals using a reactive transport model: Insights from a stream tracer test[END_REF][START_REF] Acuña | Occurrence and in-stream attenuation of wastewater-derived pharmaceuticals in Iberian rivers[END_REF]. These characteristics also explain why irbesartan peaked slightly at R2 at the flow peak (23 h, Fig. 3a), when suspended solids reached their maximum load (data not shown).

On the contrary, ibuprofen and fenofibric acid dramatically increased at first flushes, peaking at 18 h at R0 (Fig. 3b), just as PACs of the dissolve fraction behaved in our 2013-14 study. They were also highly negatively correlated with axis 2 in the PCA plot. These pharmaceuticals therefore behaved as highly hydrophilic molecules, even though their molecular weights are 206.28 g/mol and 360.83 g/mol, respectively, their log K ow are 5.19 and 3.97, and their K oc are higher (3,800 and 3,400) than those of the molecules that maintained their concentration levels during first flush (Table 4). Thus, it does not seem PACs behavior can be fully explained by the parameters described by [START_REF] Pal | Impacts of emerging organic contaminants on freshwater resources: Review of recent occurrences, sources, fate and effects[END_REF] particularly for hydrophilic compounds which are known to be affected principally by the variation of environmental conditions, such as temperature or dissolved oxygen [START_REF] Acuña | Occurrence and in-stream attenuation of wastewater-derived pharmaceuticals in Iberian rivers[END_REF]. In this case, other attenuation processes besides sorption should be taken into consideration, such as biotransformation or photolysis [START_REF] Radke | Dynamics and attenuation of acidic pharmaceuticals along a river stretch[END_REF][START_REF] Riml | Evaluating the fate of six common pharmaceuticals using a reactive transport model: Insights from a stream tracer test[END_REF][START_REF] Baena-Nogueras | Degradation kinetics of pharmaceuticals and personal care products in surface waters: photolysis vs biodegradation[END_REF].

Conclusions and Perspectives

First foul flush multicontamination events are known to occur during heavy rainfalls due to Combined Sewer Overflows (CSO), which are recurrent phenomena largely impacting the water quality of coastal Mediterranean rivers. As expected, dissolved Pharmaceutical Active Compounds (PACs) follow the same trend as other contaminants increasing rapidly during first flushes. However, we have demonstrated for the first time that CSOs will actually have little impact on PACs in the total fraction. First flush events during heavy Mediterranean rainfalls will maintain high drug concentrations for slightly longer, counteracting their dilution due to increased river flow. Our high frequency sampling at two strategically located river stations during a heavy rainfall event allowed for the precise identification of pharmaceutical sources, fate and transport. Our results support previous studies on the major contribution of point sources (WWTPs, illicit discharges and hospitals) compared to diffuse sources (runoffs) of PACs. Moreover, we could attribute the origin of medical residues downstream urban areas to CSOs from the left river bank and, with the use of gadolinium (a MRI contrast agent) as a specific marker, and pinpoint the heavy discharge of medical drugs on this side of the river to a nearby public hospital. PACs behavior could not be fully explained by their chemical properties or their molecular weight, highlighting the need to conduct studies at a fine scale on well-characterized river models to better track the fate of PACs in natural environments. A comparison of fluxes at different seasons highlighted the comparatively high contribution of recurrent heavy rain events in the Mediterranean to pharmaceuticals in aquatic environments.

Finally, if our hypothesis of illicit sanitary inputs upstream from Perpignan rejects is confirmed, release of emerging drugs during heavy rainfalls could be a good marker for these inputs.

River water quality could clearly be improved by implementing quaternary water treatment processes in WWTPs; such processes include ozonation [START_REF] Ternes | Ozonation: A tool for removal of pharmaceuticals, contrast media and musk fragrances from wastewater?[END_REF][START_REF] Rúa-Gómez | Upgrading of wastewater treatment plants through the use of unconventional treatment technologies: Removal of lidocaine, tramadol, venlafaxine and their metabolites[END_REF][START_REF] Antoniou | Required ozone doses for removing pharmaceuticals from wastewater effluents[END_REF], photocatalysis [START_REF] Dalrymple | Review -Removing pharmaceuticals and endocrinedisrupting compounds from wastewater by photocatalysis[END_REF] or activated carbon [START_REF] Rúa-Gómez | Upgrading of wastewater treatment plants through the use of unconventional treatment technologies: Removal of lidocaine, tramadol, venlafaxine and their metabolites[END_REF]. Another solution would be to use membrane bioreactors as the primary treatment [START_REF] Verlicchi | Occurrence of pharmaceutical compounds in urban wastewater: Removal, mass load and environmental risk after a secondary treatment-A review[END_REF][START_REF] Fernandez-Fontaina | Biodegradation kinetic constants and sorption coefficients of micropollutants in membrane bioreactors[END_REF], but this would need to be implemented from the time of WWTP construction. According to our results, identifying major PACs sources and enforcing adapted water treatments directly at those sources would be a more cost-efficient alternative in coastal Mediterranean aquatic environments. Although further monitoring at strategic points with high frequency sampling is necessary to better understand the fate, sources and behavior of PACs in aquatic ecosystems, governments must begin to accept responsibility for the presence of at least prescription drugs in the environment, as they can no longer be considered "emerging contaminants".

  obtained the same results for dissolved PACs dynamics during storm events in Luxembourg. Reoyo-Prats et al. 2017 demonstrated that multiple pollutants (nutrients, metals, pesticides and fecal indicators) reached concentrations never before recorded for the Têt River at exactly the same time during this rain event in 2013. They linked this pollutant cocktail to first foul flush of Combined Sewer Overflows (CSO). First flushes of wastewater and insewer cumulated sediments are a source of major pollutants

Acknowledgments

Experiments have been funded by BQR 2015 UPVD and DEBi2Micro (EC2CO 2016 and 2017 CNRS INSU) projects to CP and supported by a doctoral grant to BRP from Ecole Doctorale Energie et Environnement (E2 -ED 305 UPVD). We are thankful to O. Verneau, C. Menniti, J. Sola, B. Charriere and N. Delsaut (CEFREM) for helping us during sampling and to O.

Verneau, W. Ludwig and S. Heussner (CEFREM) for their aid and encouragement on this project. We also thank N. Faure and A. Telouk (Groupe Carso, France) for their assistance in writing protocols analyses and J-F. Lluch and J. Truffery (Veolia-Eau Perpignan) for data on sewer system volumes during the 2015 rain event. We are in debt to four anonymous reviewers for their critical review of an early version of this manuscript that highly improved its final version.


N Levernier

H Rouault

Analytical properties of UMAP dimensionality reductions

Keywords:

Dimensionality reduction techniques are essential tools for simplifying and interpreting highdimensional datasets by mapping them to a lower-dimensional space. Among existing methods, Uniform Manifold Approximation and Projection (UMAP) has recently gained a huge popularity, being applied in contexts as diverse as transcriptomics, machine learning, image processing or thermodynamics. However, understanding of this method is still sparse. Here we gives analytical prediction for the UMAP projection of well defined datasets including pure noise, binary data and one-dimensional signal. For the latter, we uncover a phase transition in the lateral spreading of the projected points. We hope that our results could help improving data analysis using non-linear dimensionality reduction in various fields.

I. INTRODUCTION

Dimensionality reduction techniques are essential tools for simplifying and interpreting high-dimensional datasets by mapping them to a lower-dimensional space. Several popular algorithms exist, including principal component analysis (PCA), t-distributed stochastic neighbor embedding (t-SNE) [1], isometric mapping (ISOMAP) [2], and Uniform Manifold Approximation and Projection (UMAP) [3]. Among them, UMAP (Uniform Manifold Approximation and Projection) has gained popularity due to its ability to process highdimensional data faster and more accurately than other techniques.

UMAP is a newer algorithm that has gained popularity due to its ability to process highdimensional data faster and more accurately than other non-linear dimensionality reduction techniques. It uses a combination of graph theory and manifold learning to identify the underlying structure in the data and project it into a lower-dimensional space. This makes it an ideal tool for visualizing complex data and identifying patterns or clusters that may be difficult to detect otherwise.

The power and limitations of UMAP have been studied in several research papers, and it has been shown to be highly effective in a variety of applications, us diverse as transcriptomic profiles [4,5], evolution of microglial morphology [6] or pan-viral interactome of Sars-COV2 [7] in biology, and free-energy landscapes [8] and phase transitions in physics [9]. However, there is still much to learn about the properties of UMAP, especially in well-defined conditions. Therefore, in this article, we propose to bridge this gap by studying UMAP's properties analytically under three well-defined conditions: pure Gaussian noise, discrete clusters, and one-dimensional datasets.

Through this study, we aim to provide a deeper understanding of the properties of UMAP's dimensionality reduction, and how it can be used to improve data analysis in various fields. By investigating the systematic properties of the dimensionality reduction obtained through UMAP, we can better understand its strengths and limitations and how to optimize its performance in different applications. Among our analytical results, we uncover phase transitions in the UMAP projection of one-dimensional signal, of different nature from transitions observed when using PCA [START_REF] Lesieur | IEEE International Symposium on Information Theory (ISIT)[END_REF].

A. Optimization function

As stated in the article presenting UMAP [3], the algorithm can be seen as an optimization problem with the following function:

L UMAP = - i =j v ij log(w ij ) + (1 -v ij ) log(1 -w ij ) (1)
In the latter, the weights v ij measures the proximity of the points in the initial space, whereas w ij measures the proximity in the space of embedding. In [3], the authors have chosen

v ij = 2 exp - ||x i -x j || -ρ σ -exp -2 ||x i -x j || -ρ σ (2)
w ij = 1 1 + a||y i -y j || 2b (3)
where x and y denote respectively the position of point i in the initial and embedding space.

The optimisation is then performed on the embedded positions y.

Later, it has been seen that the actual optimization problem UMAP solves is slightly different [START_REF] Damrich | Advances in Neural Information Processing Systems 34: Annual Conference on Neural Information Processing Systems 2021[END_REF]. Here however, we will rely on the optimization function 1.

II. DIMENSIONALITY REDUCTION OF PURE NOISE DATA AND CONNEC-TION TO RANDOM MATRIX THEORY

A. No disorder in the coupling constants

First, we consider a dataset in which v ij = v is the same for all pairs of points (all pairs are at the same distance relative to each other, which is asymptotically true for a gaussian cloud in high dimension). In such case, the function to minimize is:

L UMAP = - i =j v log(w ij ) + (1 -v) log(1 -w ij ) (4)
= i =j v log 1 -w ij w ij -log(1 -w ij ) (5)
We focus on the large N limit, so that we can replace the sum by an integral with the introduction of the density ρ of points in the low dimensional space. With such a manipulation, one should add a self-interaction energy which turns out to be negligible in the thermodynamic limit in the presence of long-range interactions [START_REF] Dyson | [END_REF]13].

Let assume a projection in dimension d = 2. Using complex notation z = y x + iy y , one has to minimize over ρ:

L UMAP = dzdz ρ(z)ρ(z ) (v -1) log(a|z -z | 2b ) + log(1 + a|z -z | 2b ) + λ dzρ(z) -1 (6)
where the last term is a Lagrange mupltiplier enforcing the normalization of the density ρ.

We assume that the integral runs on a compact support, that should be a disk by symmetry.

Cancelling the functional derivative leads to:

0 = dz ρ(z ) (v -1) log(a|z -z | 2b ) + log(1 + a|z -z | 2b ) + λ (7)
Let now apply the Laplacien ∂ z ∂ z :

0 = (v -1)bπ dz ρ(z )δ(z -z ) + dz ρ(z ) ab|z -z | 2b 1 + a|z -z | 2b-2 (8) = (v -1)bπρ(z) + dz ρ(z ) ab|z -z | 2b-2 1 + a|z -z | 2b (9)
We now assume that we can neglect the a|z -z | 2b expression in the denominator in the second term (we explain the validity afterwards) and we choose b = 1:

0 = (v -1)πρ(z) + a (10)
So, on the compact support, ρ is uniform with density ρ(z) = a π(1-v) . We thus get the radius of the disk R = (1 -v)/a (Fig. 1A-B). Interestingly, this problem turns to be the same than the distribution of eigenvalues of complex matrices with independent and identically distributed entries (the Girko's circular law) [START_REF] Mehta | Random matrices[END_REF]. Surprisingly enough, the way we have recovered the law above is, to the best of our knowledge, new. Neglecting the term in the denominator is thus valid for v close to 1. If not, the attracting force is not strong enough and the cut-off in the repulsion plays a role. For b = 1, we expect the density to be not perfectly flat, but the size of the disk should scale the same way.

B. Disordered coupling constants

We now turn to the case of disordered coupling v ij . This is the usual regime of use of UMAP, which substract the typical large distance between pairs of points -coming from dimensionality curse-to each distance. More precisely, we assume that points of the initial space are i.i.d, with any component being drawn from a normal distribution (up to a rescaling, we can choose the variance to be 1). In the large D limit, the distribution of radial squared distances is then normal with mean D and standard variance 3D. We note r the squared radial distance. We also take ρ uniform. To compute v(x, x ), we first compute:

||x -x || 2 = D i=1 x 2 i + x 2 i + 2x i x i r + r (11)
where the crossed term is negligible in the high D limit. Then, we get

v(r, r ) = 2 exp - √ r + r -ρ σ -exp -2 √ r + r -ρ σ (12) = 1 - ( √ r + r -ρ) 2 σ 2 (13)
after expanding for large σ sot that weights are close to 1. Noting u = r -m, u = r -m and ρ 2 = 2m -δ we get (see Appendix):

1 -v(u, u ) = 4m + δ 2 -2δ(u + u ) + (u + u ) 2 8mσ 2 = A + B(u + u ) + C(u + u ) 2 (14)
Indexing each point by its initial radial distance r, one has to minimize (after expansion of the repulsive term and a = b = 1):

L UMAP = dzdz drdr ρ(z, r)ρ(z , r ) (v(r, r ) -1) log(|z -z | 2 ) + |z -z | 2 (15)
where ρ stands for the joint distribution in the initial and final space (we forget the Lagrange multiplier because we will directly enforce the normalisation). Note that the polar integration on radial distances r has been replaced by a cartesian integral because the distribution of r is peaked very far from r = 0 due to the large-D hypothesis.

In this disordered case, points are not equivalent anymore: points with large initial r will tend to be further of all the other points, and thus less attracted than average, whereas small r will tend to be attracted by more points (Fig 1C). Hence, it is quite understandable that after optimization in the low dimensional space, points with small initial r will be put in the center of the cloud and points with large initial r on the periphery. We thus expect a strong coupling between initial and final radial distance. To make the calculation feasible, we make the following ansatz:

ρ(z, r) = ρ 0 (r)δ(|z| -R(r))(2πR(r)) -1 (16)
meaning that any point with initial radial distance r will be placed at a radial distance R(r) in the space of embedding. ρ 0 (r) denotes the distribution of the initial distances r. The integration on z is now straightforward using 2π 0 log(x 2 + y 2 -2xy cos θ)dθ = 2π log(max(x, y) 2 ):

L UMAP = drdr 2ρ 0 (r)ρ 0 (r ) (v(r, r ) -1) log(max((R(r)), R(r ))) + R(r) 2 + R(r ) 2 = 4 ∞ -∞ r -∞ drdr ρ 0 (r)ρ 0 (r )(v(r, r ) -1) log(R(r)) + 2 ∞ -∞ drρ 0 (r)R(r) 2 (17)
We can now differentiate with respect to the unknown function R, so that the optimal function satisfies:

r -∞ dr ρ 0 (r)ρ 0 (r )(v(r, r ) -1) 1 R(r) + ρ 0 (r)R(r) = 0 ( 18
)
The solution is thus

R(r) = r -∞ (1 -v(r, r ))ρ 0 (r )dr (19)
It is straightforward to check that this expression is perfectly consistent with the previous case of uniform v(r, r ) = v, that is a uniform distribution of points stacked on a disk of radius √ 1 -v. Let now compute this for the case of v(r, r ) given by ( 14). We get:

R(u) 2 = u -∞ du e -u 2 2s 2 √ 2πs (A + B(u + u ) + C(u + u )) = (A + Bu + C(u 2 + s 2 ))Φ(u/s) - s(B + 3Cu) √ 2π e -u 2 2s 2 ( 20
)
where Φ is the cumulative distribution function of standard gaussian variable. This result is in striking contrast with the previous one: the fact that v's are disordered lead to a cloud which is not anymore of compact support (even though most points are densely packed inside a circle of radius B + Cs).

We can finally compute the local density of the embedded cloud as a function of the radial distance d:

f (d) = dr |z|=d dzρ(r, z)/(2πd) = ρ 0 (u) 2πR(u)R (u) (21)
with u = R -1 (d) (see Fig. 1E). On the other hand, the cumulative function (the fraction of points at radial distance smaller than d) is simply:

C(d) = R -1 (d) 0 ρ 0 (r) dr (22)
Let now turn to the case where ρ is not uniform, but depends on the closest neighbour of each point. ρ 2 i is correlated to r 2 i and we write ρ

2 i = r 2 i + m -α + δ i = 2m -α + u i + δ i
, where δ i is approximately normally distributed of mean zero and with zero correlation with r 2

i . The calculation is mostly similar to the previous one and one obtain:

1 -v(u i , u j ) = 4m + α 2 -α(u i + u j -δ i -δ j ) + 3/2(u i + u j -δ i -δ j ) 2 8mσ 2 (23)
Hence, the problem is now turned into the mapping of the variable v i = 2u i -δ i from high to low dimension (See Fig. 1D andF).

III. DIMENSION REDUCTION OF BINARY NOISY DATA

We now turn to the use of UMAP as a clustering tool. We focus on the case of data without any noise. Hence, let consider the case where there are two populations of points, with attractive weights v 11 = v 22 = v and v 12 = v -δv. For δv = 0, all points are equivalent and we are back to the situation of the previous part without noise: the embedding lead to a disk of radius R = √ 1 -v. For non-zero δv, this disk will break in two hemispheres H 1 and H 2 , one for each community, separated by a distance ∆ (see Fig. 2A-C). We expect that at first order, the diameter of both hemispheres is the same than for δv = 0. We want to estimate this distance as a function of δv. One has to minimize the following quantity over ∆:

L = H 1 H 2 dzdz (v 12 -1) log |z -z | 2 + |z -z | 2 (24)
Let note H * 2 the hemisphere shifted by -∆ so that the reunion of H 1 and H * 2 is a disk. We can now write:

L = H 1 H * 2 dzdz (v 12 -1) log(|z -z | 2 -2∆(x -x ) + ∆ 2 ) + |z -z | 2 -2∆(x -x ) + ∆ 2 = (v -δv -1) H 1 H * 2 dzdz log |z -z | 2 - 2(x -x ) |z -z | 2 ∆ + 1 |z -z | 2 + 2(x -x ) 2 |z -z | 4 ∆ 2 + H 1 H * 2 dzdz |z -z | 2 -2∆(x -x ) + ∆ 2 (25)
with z = (x, y), z = (x , y ) and x chosen along the axis of shift. We know that for δv = 0, ∆ = 0 is a minimum. Thus:

H 1 H * 2 dzdz (1 -v) (x -x ) |z -z | 2 -(x -x ) = 0 (26)
Finally, at second order in ∆:

L = C∆ 2 + 2δv H 1 H * 2 dzdz (x -x ) |z -z | 2 ∆ (27) with C = H 1 H * 2 dzdz v-1 |z-z | 2 1 + 2 (x-x ) 2 |z-z | 2 + 1 > 0.
This leads to a trancritical transition:

∆ = 4πR 3 3C δv = 4πR 3C δv (28)
where we have used (26) to simplify the numerator and C = C/R 2 is a pure number (see

Fig 2D)

. Thus, as soon as δv is of order of unity, both clouds are very well separated, with distance between them comparable to their size. For the sake of completeness, we also give the scaling behaviour for very dissimilar communities (δv 1): ∆ ∼ 1/ √ δv.

IV. DIMENSION REDUCTION OF 1D CONTINOUS SIGNAL AND PHASE TRANSITION A. No noise

Let first investigate the case of perfect signal, with no noise. Let {x N } be the set of points in the ambient space of dimension D, and let note {y N } their UMAP-projection in the space of dimension d < D. Here, we focus on the case where the points {x N } lie perfectly on a segment of length 1, and are equally spaced. We can also assume d = 1, and we note

y i = g(x i ) = g(i/N ) for 0 ≤ i < N .
To compute the length of the obtained projection, we look for a solution g(x) = αx that minimizes L UMAP , in the limit of small σ (i.e. for large α). Performing a change of variable and using the symmetry of the integrand, we have to compute:

L UMAP = -2 1 0 dv v 0 duK(u) log(aα 2b u 2b ) + 2 1 0 dv v 0 du log 1 + 1 aα 2b u 2b (29)
The first integral is easy to evaluate: where C does not depend on α.

I 1 = -2 1 0 dv v 0 duK(u) log(aα 2b u 2b ) = -4b log α 1 0 dv v 0 duK(u) + C = (6bσ + O(σ 2 )) log α + C (30)
The second part gives, for large α:

I 2 = 2 αa 1 2b +∞ 0 du log(1 + 1/u 2b ) = 2K αa 1 2b (31)
For UMAP default parameters, this integral is numerically evaluated to K = 3.196. Putting all together, deriving with respect to α to find the extremum, we obtain the length of the UMAP projected line:

α = K 3ba 1/(2b) σ (32)

B. With noise

We now assume that uncorrelated gaussian noise of amplitude / √ 2 is added on each point (the √ 2 is only here to make the variance of the difference between two realisations of the noise to be ). We are in interested in the regime D 1 and D 2 1 (so that the 1-d signal can be extracted from the noise). The Umap projection typically creates an elongated cloud of length α and thickness δ. Let first assume general σ and ρ. We focus in the limit of large α (which corresponds to small sigma, similarly to the no noise case).

To get an analytical estimate of α and δ, we will make some assumptions and approximation. The first one is that we can minimize the annealed functional, which is expected to be exact in the thermodynamic limit. The second one is that the density of the cloud along the transverse axis is gaussian. This is not exact but should give the same scaling with the

problem parameters. Let note m = D 2 , s = √ 3D 2 .
First, after introducing η, the radial noise in the initial space η = i=1...D-1 η 2 i , which is a gaussian v.a. of mean m and variance s in the large D limit, we get:

L = - 1 0 dx x 0 dy ∞ -∞ e -(η-m) 2 2s 2 N dη ∞ -∞ e -ξ 2 4δ 2 N dξ [(-2e -y 2 +η-ρ σ + e -2 y 2 +η-ρ σ ) log(a(α 2 y 2 + ξ 2 ) b ) -log a(α 2 y 2 + ξ 2 ) b 1 + a(α 2 y 2 + ξ 2 ) b (33)
For large α, we can neglect side effects and obtain a translation-invariant problem:

L = - ∞ 0 dy ∞ -∞ e -(η-m) 2 2s 2 N dη ∞ -∞ e -ξ 2 4δ 2 N dξ (-2e -y 2 +η-ρ σ + e -2 y 2 +η-ρ σ ) log(a(α 2 y 2 + ξ 2 ) b ) -log( a(α 2 y 2 + ξ 2 ) b 1 + a(α 2 y 2 + ξ 2 ) b ) (34)
Now, we develop the expression in the limit of large noise. We have to compute I = ), but the following calculations is more general. In the limit we are interested in, σ is always much smaller than m, so that: We see that after averaging over the initial noise, we get an expression close to the one without noise : simply a change from σ to σ = √ mσ and a gaussian kernel in place of an exponential one with prefactors c in front of exponential terms. Finally, one has to minimize:

I = ∞ -∞ dη e -η 2
L = ∞ 0 dy ∞ -∞ dξ e -ξ 2 4δ 2 N b(2ce -y 2 2σ 2 -c 2 e -y 2 σ2 ) log(α 2 y 2 + ξ 2 ) -log a(α 2 y 2 + ξ 2 ) b 1 + a(α 2 y 2 + ξ 2 ) b (37)
Unfortunately, we cannot get explicit expressions with a gaussian kernel. We simplify a bit by replacing by a square of size σ and amplitude v = 2c -c 2 . We also restrict to b = 1 and a = 1 (which amounts to rescale length in low dimenion space). Thus:

L = 1 α ∞ 0 dy δ 0 dξ δ vH(y/ασ) log(y 2 + ξ 2 ) -log y 2 + ξ 2 1 + y 2 + ξ 2 (38)
Let pose S = ασ and divide the former by σ. We now optimize on S and δ the function:

L = 1 S ∞ 0 dy δ 0 dξ δ vH(y/S) log(y 2 + ξ 2 ) -log y 2 + ξ 2 1 + y 2 + ξ 2 (39)
Now we can compute everything. The first term with the Heaviside gives:

L att = v log(S 2 + δ 2 ) + δ S arctan S δ + S δ arctan δ S (40)
The second one gives

L rep = - π 2 δ S 1 -1 + 1/(2δ 2 ) -argsh( √ 2δ)/(2δ 2 ) (41)
Now, it is not hard to see that if v > 1, then δ opt = 0. This is due to the fact that at short distance points always attract more then they repulse so that points with same signal collapse. Then, we can easily compute S = π 2 √

2v Then, we can get the critical behaviour close to v = 1. Let introduce e = 1 -v. At first order we get δ = 3e 2 √ 2-2/(πS) . At this order, the value of S is unchanged. We have thus obtained a phase transition in the UMAP projection, with a transcritical bifurcation (due to the fact that there is no δ ↔ -δ symmetry): for strong attraction v, one gets a line of strictly zero thickness, whereas for large attraction one obtains a thick elongated cloud.

We can estimate the aspect ratio of the cloud. Our results above are for general ρ and σ, but we will focus only to their values with UMAP default settings. Then, c is only controlled by ∆ρ, so that:

a = δ/α = cst.σe = cst. 3/2 D 1/4 σ 1/2 (42)
As for σ defined as the distance to the k-th neighbour minus ρ, we have two regimes:

σ = d k -ρ = (k/N ) 2 + D 2 - √ D + ∆ρ = ∆ρ + (k/N ) 2 √ D for √ D k/N (43) = k/N - √ D for √ D k/N (44)
The large noise regime can be split in two: either ∆ρ

(k/N ) 2 √ D and σ = (k/N ) 2 √ D or ∆ρ (k/N ) 2 √ D
and σ = ∆ρ k . Finally, we get:

a = cst. D 1/4 σ 1/2 for √ D k/N D 1/4 (45) = cst. 2 D 1/2 k/N for k/N √ D k/N D 1/4 (46) = cst. 3/2 D 1/4 k/N for √ D k/N (47)
V. DISCUSSION

Our analytical study of UMAP dimensionality reductions brings insight into the structure of the obtained patterns. By considering the limit of a large number of points, we were able to use mean-field approaches and explore the continuous density of points. For instance, a pure Gaussian noisy input (without any structure) produces a cloud of size of order one with nearly uniform density. This should prompt users of this technique to pay attention to the size of the obtained reductions.

At the same time, the discrete cluster analysis shows that the technique is highly sensitive to minor differences between clusters. This effect results from the long-range interactions between dimension reduction points, which are specific to UMAP. While the entropic effect has a negligible effect on the average density, it remains present and has a strong impact on cluster segregation.

Finally, our mean-field approximation can predict the shape (length and thickness) of clouds resulting from the dimensionality reduction of high-dimensional signals containing only a one-dimensional manifold structure. We were able to characterize a transcritical phase transition as a function of the interaction strength, v, that explains the effectiveness of UMAP in revealing such structures.

We believe that our analysis provides a strong foundation for assessing the significance of UMAP dimensionality reduction in the presence of high levels of noise, which is often the case in realistic data analysis scenarios. Furthermore, this mean-field approach can extend beyond the simple scenarios that we have considered here.

FIG. 1 :

 1 FIG. 1: UMAP dimensionality reduction for pure noise datasets. A. Example of a two dimensional UMAP dimensionality reduction for constant v. B. Radius of the dimensionality reduction as a function of √ 1 -v, with simulation (dots) and our prediction (line). C. Same as A but for a gaussian input in D=512 dimensions with non-uniform ρ. D. Prediction of the mapping function R in the case of non uniform ρ. E. Comparison of our prediction of the radial density of the cloud (line) with simulations (histogram). F. Same as E but for uniform ρ.

FIG. 2 :

 2 FIG.2: UMAP dimensionality reduction for binary datasets. A, B and C. UMAP dimensionality reduction for binary signal with low, moderate and strong dissimilarity (δ intra = 0.9 and δ inter = 0.89, 0.8, 0.5 respectively. D. Distance between hemispheres in the low dissimilarity regime, with simulation (dots) and prediction (line).

∞-∞ 2 Ne

 2 dη e -(η-m) 2 2s -y 2 +η-ρ σ . Let note ρ = √ m -∆ρ. Within UMAP choice of parameters, we have ∆ρ = k log(N D 1/4

  exponential we get the same value by replacing σ2 to σ2 /2 and c by c 2 .

FIG. 3 :

 3 FIG. 3: UMAP dimensionality reduction for 1D continuous signal datasets. A, B, C. Example of a two dimensionality reductions for a noise amplitude equal to 0.05, 0.02, 0.005 respectively. B. Radius of the dimensionality reduction as a function of √ 1 -v, with simulation (dots) and our prediction (line). D. Length of the dimensionality reduction (points: simulation, line: prediction).E. Thickness of the dimensionality reduction (points: simulation, line: prediction).

Acknowledgments

The project leading to this publication has received funding from the "Investissements d'Avenir" French Government program managed by the French National Research Agency (ANR-16-CONV-0001) and from Excellence Initiative of Aix-Marseille University -A*MIDEX.

Appendix A: Calculation of v(u, u ) in the pure noise case

In the case of uniform ρ, the calculation is straightforward. Let first note d 2 ij = 2m + i j. By expanding the exponential terms for large σ with the assumption that δ and ij are small compared to m, we obtain:

This is the equation ( 14) in the main text.

The calculation for non-uniform is similar, except that one has to introduce a δ depending on the considered point δ i .

Appendix B: Code availability

As noted previously, our implementation of UMAP stictly follows the optimization function 1, which is different from the initial implementation [3]. Our implementation has been coded in the Rust programming language and is available at https://gitlab.com/ rouault-team-public/analysis/umaprs.