Aquatic risks from human pharmaceuticals—modelling temporal trends of carbamazepine and ciprofloxacin at the global scale

Despite the worldwide presence of pharmaceuticals in the aquatic environment, a comprehensive picture of their aquatic risk (AR) at the global scale has not yet been produced. Here, we present a procedure to estimate ARs of human pharmaceuticals at a freshwater ecoregion level. First, we predicted country- and year-specific per capita consumption with a regression model. Second, we calculated spatially explicit freshwater concentrations via a combination of mass balance models, addressing the pharmaceutical’s fate in respectively humans, wastewater treatment plants and the environment. Finally, we divided the freshwater concentrations at the level of individual freshwater ecoregions with the regulatory limit value derived from toxicity tests to come to an ecoregion-specific AR. We applied our procedure to model time-trends (1995–2015) of ARs of carbamazepine and ciprofloxacin, two widely detected and regulatory relevant human use pharmaceuticals. Our analysis of carbamazepine and ciprofloxacin showed that ARs, due to exposure to these human pharmaceuticals, typically increased 10–20 fold over the last 20 years. Risks due to carbamazepine exposure were still typically low for the time period assessed (AR < 0.1), although some more densely populated and/or arid ecoregions showed higher ARs (up to 1.1). Risks for ciprofloxacin were found to be much higher with ARs larger than 1 for 223 out of 449 freshwater ecoregions in 2015. Comparison with measured concentrations in ten river basins showed that carbamazepine concentrations were predicted well. Concentrations of ciprofloxacin, measured in four river basins, were, however, generally underestimated by our model with one to two orders of magnitude. We conclude that our procedure provides a good starting point to evaluate ARs of a wide range of human pharmaceuticals at the global scale.


Introduction
Active pharmaceutical ingredients (APIs) are consumed in large quantities (e.g. Van Boeckel et al 2014) that are not completely removed during their passage through the human body and wastewater treatment plants (WWTPs). Consequently, their residues can be transported into the environment, and many APIs have been detected in surface waters globally (Hughes et al 2013, Guang et al 2014, aus der Beek et al 2016. Adverse effects on aquatic species have been observed at environmentally relevant concentrations for a wide range of (pharmaceutical classes of) APIs, an extensive overview of which is provided by Brausch et al (2012). In addition, recent studies have shown that chronic exposure to carbamazepine, an anti-epileptic drug, leads to altered feeding behaviour and reduced egg viability in zebrafish (da Silva Santos et al 2018), and lowers the reproductive output in crustaceans (Oropesa et al 2016). Antibiotics are a specific group of APIs of potential ecological concern since all major nutrient cycles in the environment depend on the activity of bacterial communities. Additionally, antibiotics entering WWTPs can affect bacterial communities used for biological degradation and consequently decrease their efficiency to remove other pollutants from the water (e.g. Schmidt et al 2012).
The aquatic risk (AR) of pharmaceuticals has been acknowledged as an emerging environmental problem requiring regulatory and scientific attention (Küster and Adler 2014). Because of technological and practical limitations, monitoring cannot provide a complete overview of the global situation. Therefore, fate models have been applied to estimate concentrations of APIs in several river basins in Europe . These approaches are, however, generally limited in their spatial coverage, and despite the worldwide presence of APIs in the aquatic environment a comprehensive picture of ARs at the global scale has yet to be produced.
Here, we present modelled time-trends  of ARs of two APIs, i.e. carbamazepine and ciprofloxacin, in 449 aquatic ecoregions worldwide (Abell et al 2008). We chose to model carbamazepine and ciprofloxacin, a second-generation broad-spectrum antibiotic, because they are widely detected in surface waters (Hughes et al 2013, aus der Beek et al 2016), and because they were both identified as candidates for the Watch List of the European Water Framework Directive (Carvalho et al 2015). Our model predicts API emissions and aquatic concentrations at a 0.5°by 0.5°grid-cell resolution, and ARs are derived using publicly available predicted no effect concentrations (PNEC).

Consumption
Global and temporal coverage of API consumption was achieved through development of linear mixed effects models (LMMs) that incorporate socio-economic and demographic predictors expected to influence pharmaceutical consumption (Kookana et al 2014). Our LMMs extend current approaches that extrapolate between countries and years based on population dynamics alone They were built for each API using a database containing 227 and 567 data points for carbamazepine and ciprofloxacin, respectively, each representing per capita total consumption for a unique combination of country and year, spanning from 1995 to 2015. Data were collected from open literature and online databases. Additionally, data on ciprofloxacin consumption were received from The European Surveillance System-TESSy, released by (European Centre for Disease Control (ECDC). The final database contained data on carbamazepine consumption in 69 countries covering all continents (32 in Europe,19 in Asia,4 in Central-/North-America, 4 in Africa, and 2 in Oceania), and data on ciprofloxacin consumption in 49 countries covering all continents except for Africa (32 in Europe, 5 in Asia, 7 in South-America, 4 in Central-/North-America, and 1 in Oceania). The supporting information provides a description of the data collection procedure including a list of countries and the years with consumption data available (section SI1).
Based on their supposed correlation with API consumption, predictors selected were the Human Development Index (United Nations Development Programme 2016), the ageing of the population (The World Bank 2016), and the time since first marketing of API of interest. The models were constructed using the 'lme' function of the package 'nlme' (Pinheiro et al 2016) in the R environment (R Development Core Team 2017). They were based on the log-transformed per capita consumption of carbamazepine and ciprofloxacin, incorporating fixed effects from all predictors including their quadratic terms. Random effects were included per country on the regression's intercept and on the slope over the years. Inclusion of these random effects was justified based on the decrease in Akaike Information Criterion score (AIC) compared to the full fixed-effects model. We employed the function 'dredge' of the package 'MuMIn' in the R environment (Barton 2015) to select the most parsimonious models based on their AIC score. Independence of predictors was checked and Variance Inflation Factors were all below 2, well below the threshold of 5 for unacceptable multi-collinearity as proposed by Zuur et al (2009). The model structure of the final LMMs is given by equation (1) In which M is the per capita consumption of the API of interest (mg/cap/yr), HDI is the Human Development Index (-), AGE is the ageing of the population (% of population65 years), YR_M is the time since first marketing of API of interest (years), and (1+YR_M|Country) reflects the random effect structure on the regression's intercept and slope over YR_M per country. Coefficients of predictors are represented by β, with subscripts 0 for the intercept, 1 for linear terms, and 2 for quadratic terms. Error term ε represents the residual variance not explained by the model.
Model performance was assessed via their marginal and conditional pseudo-R 2 values, calculated according to Nakagawa and Schielzeth (2013) and Johnson (2014). Moreover, the predictive capacity Q2 of the fixed-effects part of the models was computed via a leave-one-cluster-out cross validation (Fang 2011). While the Q 2 was smaller than the marginal pseudo-R 2 values for both models (0.34 and 0.29 versus 0.27 and 0.15 for carbamazepine and ciprofloxacin, respectively), we considered this decrease sufficiently small to rule out unacceptable overfitting for both models. Table 1 contains the characteristics of the models, including their regression coefficients. Moreover, the supporting information contains partial dependence plots per combination of predictor and API, to aid in model interpretation (figure SI1.1 is available online at stacks.iop.org/ERL/14/034003/ mmedia).
In our model calculations for API consumption, each unique combination of country and year was preferably assigned a consumption value extracted from our database. If this was not possible, we applied the conditional models to estimate the consumption, provided that our database contained data for the country of interest in at least one other year. For countries not included in the database, we predicted consumption using the fixed-effects part of the models only.

Aquatic concentrations
The hydrological model PCRGLOBWB (PCRaster GLOBal Water Balance model) calculates the global hydrology and water resources at a grid-cell resolution of 0.5°by 0.5°(Van Beek and Bierkens 2008). Here we coupled the discharge from PCRGLOBWB with consecutive spatially explicit emission modules for consumption, human metabolism, WWTP-treatment, and environmental fate of APIs to calculate predicted environmental concentrations (PEC) for the years 1995-2015. Steady-state concentrations in surface water were estimated per spatial grid cell i in PCRGLOBWB via equation (2): In which E i is the emission into cell i (kg d −1 ; see equation (3)), Q i is the mean annual discharge in cell i (m 3 d −1 ), and V i is the water volume in cell i (m 3 ), to derive a cell-specific downstream PEC (PEC i ; g l −1 ). Five loss processes j are accounted for via cell-specific first order rate constants (k j, i ; d −1 ), being the three degradation processes biodegradation, photolysis and hydrolysis, and the two intermedia processes sedimentation and volatilization. Individual degradation rate constants were derived from test conditions, and translated into field conditions by accounting for temperature differences, sorption to suspended solids and dissolved organic carbon (e.g. Honti et al 2016) and reduced light intensity (Schwarzenbach et al 1993). Local sedimentation and volatilization rate constants were implemented via mass transport velocities between media, similar to Fantke et al (2017) and Margni et al (2004). Detailed information on the derivation of all rate constants can be found in the supporting information (section SI2).
Emissions into surface water in cell i (E i ; kg d −1 ) were calculated via equation (3): In which the per capita API consumption (M i ; kg/ cap/d) and the local urban population (P urb,i ; cap), are combined to derive the total API consumption in the urban part of cell i. PCRGLOBWB considers the nonurban rural population to have on-site treatment (e.g. local septic tanks), from which no emissions to surface water take place. After consumption, APIs are partly absorbed from the intestinal tract and subsequently metabolized. The remaining fraction of the administered dose is eliminated unchanged as parent compound from the body, either egested via faeces ( f pc,f ; -) or excreted via urine ( f pc,u ; -). Unchanged API might Table 1. Characteristics of linear mixed effects regression models for per capita consumption of carbamazepine (CBZ) and ciprofloxacin (CIP).
be transported with wastewater towards a WWTP, depending on the fraction of the urban population connected to wastewater treatment ( f conn,i ; -). In cells with a WWTP-connectivity below 100%, the fraction excreted via urine for the non-connected part was assumed to be discharged directly into surface water. The fate of APIs during wastewater treatment was estimated using the multimedia model SimpleTreat 4.0 (Struijs 2014). SimpleTreat 4.0 distinguishes between WWTPs that apply primary treatment only and those that combine this with secondary treatment with activated sludge. In equation (3), this is incorporated with Boolean parameter a i indicating the presence or absence of merely primary treatment or primary treatment combined with secondary treatment in cell i, and with the corresponding removal fractions f rem,i . Finally, the chemical inflow into cell i from the upstream cell (E up,i ; kg d −1 ) is calculated as:

Comparison with measured concentrations
Concentration estimations for carbamazepine and ciprofloxacin were compared with measured concentrations included in the global database from the Umwelt Bundesamt, the German national environmental protection agency (aus der Beek et al 2016). From this database, all records were extracted reporting on the presence of carbamazepine or ciprofloxacin in fresh surface waters, provided that the sampling year(s) and sampling location were reported. To enable a meaningful comparison with model predictions, we allocated individual records to hydrologically distinct river basins. Records were excluded if they could not be allocated to a specific river basin based on their reported sampling locations. The original literature sources were consulted for those records representing aggregate average concentrations over multiple sampling points or over multiple years. Where possible, concentrations for unique combinations of sampling location and year were extracted from the original sources. Otherwise, the records were excluded from the analysis. Finally, river basins were only included in the validation exercise if measured concentrations were available for at least five locations in the river basin, with at least 20% of them above their limit of detection. For these river basins, box-whisker plots were constructed with the left-censored data substituted via non-parametric robust regression on order statistics (ROS) using the 'cenboxplot' function of the package 'NADA' in the R environment (Lee 2017). These were compared with box-whisker plots on modelled concentrations for grid cells in the respective river basins. Cells were excluded for which zero concentrations were predicted, because they generally represent upstream areas where measurements are not performed.

Aquatic risks (ARs)
We calculated yearly ARs for carbamazepine and ciprofloxacin over a period of 21 years , for 449 freshwater ecoregions (Abell et al 2008). Ecoregions are distinct geospatial units that represent unique patterns of environmental and ecological variables. As such, ecoregions are considered particularly useful for risk assessment and conservation planning at a larger scale (Schäfer et al 2013, Szöcs et al 2017. The AR for each freshwater ecoregion b (AR b ) was calculated as the average risk-quotient over its total water volume (equation (5)).
In which w iεb is the weight of grid cell i in ecoregion b, based on the water volume of the cell. The larger the water volume, the larger the weight of that cell. The PNEC is the PNEC which is 500 ng l −1 for carbamazepine As an additional measure of risk, we also calculated the percentage of the water volume per ecoregion exceeding the PNEC (equation (6)).

Consumption
Applying the regression models from 1995-2015, we found that global carbamazepine consumption increased from an estimated 742 tonnes in 1995, to an estimated 1214 tonnes in 2015 (figure 1). Our estimation of global ciprofloxacin consumption shows a much larger increase: by 2015, 2318 tonnes of ciprofloxacin were consumed worldwide, more than seven times the estimated consumption of 298 tonnes in 1995. It should be noted, however, that these numbers represent total tonnages, and are a combined result of increases in per capita consumption as well as in population size. Indeed, the strong increase in total ciprofloxacin consumption in Asia ( figure 1(b)), is less pronounced when consumption is expressed on a per capita basis.

Concentrations
Despite its larger consumption volume and its lower degradation in the human body (table SI2.1), concentrations of ciprofloxacin, aggregated at the ecoregion level, were similar to those of carbamazepine (figure 2). Moreover, aquatic concentrations of carbamazepine and ciprofloxacin in the typical (median) ecoregion showed the same temporal trend, increasing at an average rate of 13% and 16% per year, respectively. Interestingly, ecoregions have converged over the years with respect to their carbamazepine concentrations, as indicated by the interquartile range in figure 2(a). This pattern was not found for ciprofloxacin, for which the lower bound of the interquartile range stayed below 10 −8 ng l −1 over the whole period 1995-2015 ( figure 2(b)).
The Sankey diagrams in figure 3 show the global mass flows of both APIs in 2015, with the thickness of respective ribbons reflecting the absolute amounts (tonnes yr −1 ). Dark shaded ribbons represent transport flows and light shaded ribbons represent actual removal and retention. Similar figures for the years 1995-2014 are included in the supporting information (section SI.3). Figure 2(B) shows that ciprofloxacin has a relatively large in-stream removal. This is the result of its higher removal in surface waters, mainly due to photolysis (table SI2.1), while sedimentation contributes very little to dissipation of ciprofloxacin from surface waters (<1%). In fact, no more than 4% of all ciprofloxacin emitted into fresh surface waters reaches the sea, compared to 79% of all carbamazepine. As a consequence, carbamazepine concentrations vary less throughout the river systems than do concentrations of ciprofloxacin.
Our modelled carbamazepine concentrations show good agreement with measurements in nine river basins throughout Europe and the Mississippi basin in North America, as is shown by the overlap between box-whisker plots of measured versus modelled concentrations (figure 4). In contrast with the ten river basins for carbamazepine, only four river basins were identified where ciprofloxacin was detected in more than 20% of the measurements ( figure 4). This comes as no surprise considering that typical detection limits for ciprofloxacin are in the 1-100 ng l −1 range (aus der Beek et al 2016). Two of these basins were located in Asia (Yangtze and Chao Phraya basins), one in Europe (Ebro basin), and one in North-America  (Mississippi basin), and in all four of them predicted ciprofloxacin concentrations were generally one to two orders of magnitude lower than the measurements.

Aquatic risks
In line with the increasing consumption (figure 1) and concentrations (figure 2), ARs at the ecoregion level increased over time for both APIs ( figure 5). The ARs of carbamazepine stayed, however, below 1 for all but one of the 449 ecoregions ( figure 5(b)). The ecoregions with relatively large ARs, i.e. 15 ecoregions with a value above 0.1, mainly represented densely populated areas such as western and central Europe, or water-scarce arid regions such as the Arabian and Californian peninsulas. Moreover, throughout the 21 year period the % of water volume exceeding the PNEC in any ecoregion was never larger than 50%.  . Box-whisker plots comparing measured concentrations (yellow) and modelled concentrations (green: carbamazepine CBZ; purple: ciprofloxacin CIP) in several river basins. Boxes represent interquartile range, with black horizontal line the median concentration; whiskers represent minimum and maximum. Left-censored data (<LOD) were substituted via non-parametric robust regression on order statistics (ROS) using the 'NADA' package for R (Lee, 2017). Red lines represent the largest LODs reported for that river basin. River basins were visualized using the HydroSHEDS drainage network (Lehner et al 2008).
In contrast with carbamazepine, the AR of ciprofloxacin exceeded 0.1 in 196 ecoregions and 1 in 112 ecoregions in 1995 ( figure 5(c)). These numbers further increased over time, leading to 282 ecoregions with an AR higher than 0.1 and 214 ecoregions with an AR higher than 1 in 2015 ( figure 5(d)). In 6 ecoregions, 90% or more of the freshwater volume exceeded the PNEC for ciprofloxacin. These were small and often relatively densely populated ecoregions, such as the coastal region of Israel, the Canary Islands off the African coast, and the area around Rio de Janeiro in Brazil.

Discussion
We showed that it is indeed possible to predict temporal changes in spatially explicit ARs at the global scale for human pharmaceuticals. Our modelling is, however, not without uncertainty. General practice in chemical risk assessment is to follow the consecutive steps of use, emission, environmental fate and effect. Previous research has shown that the first and last of these steps contribute most to the uncertainty in the final risk estimation (Oldenkamp et al 2016). Below, we reflect on the uncertainties encountered during the emission, effect and fate estimations, respectively.
One major hurdle to overcome when performing a sound global environmental risk assessment of APIs, is to reliably estimate their country-and year-specific consumption. Not only are consumption data scarcely and spatially not homogeneously available, consumption also varies substantially between countries and years (Van Boeckel et al 2014, aus der Beek et al 2016, Klein et al 2018). To make optimal use of the available data, we extrapolated per capita consumption to all countries globally using mixed effects regression models (table 1). The models describing carbamazepine and ciprofloxacin consumption were based on a set of readily available predictor variables (equation (1)).
The regression coefficients in table 1 and the partial dependence plots in the supporting information (figure SI1.1), indicate that consumption of both APIs increases with HDI, although this trend levels off above HDI values of ∼0.8. In addition, their consumption follows traditional trends of initial increase and subsequent decline in the years after introduction. In accordance with Zhang and Geißen (2010), carbamazepine consumption increases with ageing of the population.
Fixed effect predictors explained 29%-34% of the variance in per capita consumption, indicating that extrapolation to countries without any available consumption data is particularly uncertain. Furthermore, the uneven geographical coverage of consumption data might have resulted in additional uncertainties in our estimations for low-and middle-income countries outside of Europe. Nevertheless, our estimation of total global consumption of carbamazepine in 2007, based on the country specific estimates (968 tonnes; figure 1) is in accordance with actual sales data from IMS Health for 96% of the global market that year (942 tonnes; Zhang et al 2008). Moreover, Van Boeckel et al (2014) and Klein et al (2018) report that total antibiotic consumption was stable or had moderately decreased in high-income countries between 2000 and 2010. Klein et al (2018) further estimated an increase in the consumption rate of fluoroquinolones, the subgroup of antibiotics to which ciprofloxacin belongs, in low-and middle-income countries, while it decreased in high-income countries. This supports our estimation of global ciprofloxacin consumption, which levels off towards 2015 (figure 1). It should be pointed out, however, that our estimations of ciprofloxacin consumption in the years after 2010 were largely based on data originating from European countries as provided by ECDC. While this might have introduced a bias in our estimation of global ciprofloxacin consumption, similar temporal trends of  (Rivas and Alonso, 2011).
While PNECs should reflect the concentration at which adverse ecological effects do likely not occur (EMA, 2006), they are generally based on a limited number of ecotoxicity data derived under laboratory conditions. To account for uncertainties in the necessary interspecies and lab-to-field extrapolations, PNECs are derived via application of assessment factors to the lowest available ecotoxicity endpoint, with the magnitude of these assessment factors depending on data quality and quantity. This approach has led to a very strict PNEC for ciprofloxacin at 0.15 ng l −1 , based on a chronic No Observed Effect Concentration for the sensitive bioluminescent bacterium Vibrio fischeri, combined with an assessment factor of 10 (Załęska-Radziwiłł et al 2014). For comparison, the PNEC for ciprofloxacin when based only on the sensitivity of the standard taxonomic groups of fish, crustaceans, and algae, would be 100 000 times less strict (at 15.6 μg l −1 ; Załęska-Radziwiłł et al 2011). In fact, while our calculations indicate almost half of all ecoregions worldwide to have an unacceptable AR due to chronic exposure to ciprofloxacin, the PNEC based on the three standard taxonomic groups is not exceeded anywhere.
Finally, the comparison of predicted with measured concentrations provides confidence that our model can properly estimate concentrations of carbamazepine, but to a lesser extent ciprofloxacin. The underestimation of measured ciprofloxacin concentrations should at least partially be attributed to analytical limitations and non-representative sampling locations. An additional explanation may be that we did not consider four potential emission sources of ciprofloxacin in our modelling. Non-prescription use of ciprofloxacin might account for 19%-100% of the total consumption of antibiotics outside northern Europe and North America (Morgan et al 2011). Furthermore, veterinary use of ciprofloxacin, while not approved in either the USA (FDA, 2016) or Europe (EMA, 2016), can be widespread in other countries. Indeed, ciprofloxacin has been detected in manure from livestock (i.e. chicken, cow, and pig) originating from several Chinese provinces (Hu et al 2008, Zhao et al 2010). Moreover, ciprofloxacin concentrations can be substantial in biosolids applied to agricultural soils (Verlicchi and Zambello 2015). Biosolids can thus form a significant emission pathway when a considerable part of the excess sludge is disposed of this way (Oldenkamp et al 2013a). For example, 75% of all Spanish excess sludge was applied on agricultural soils in 2012 (Eurostat 2017). However, conservative calculations we performed with the multimedia box model SimpleBox 4.0 (Hollander et al 2016), indicated that transfer of ciprofloxacin from agricultural soils to surface waters was negligible, as ciprofloxacin adsorbs strongly to soils. Finally, we assumed that APIs emitted into rural wastewater do not reach surface waters, resulting in 100% rural retention (figure 3). In reality, however, leakage of pharmaceuticals into ground water from septic tanks and pit latrines may not be negligible (Del Rosario et al 2014, Phillips et al 2015, Schaider et al 2016, Hanamoto et al 2018, showing that this could form an additional route into surface waters that is not yet understood well enough to be accounted for at a global scale. Photodegradation is the most important environmental degradation pathway for various APIs (Wang and Lin 2014), and the extensive in-stream removal as estimated for ciprofloxacin can be largely contributed to this process (figure 3). The photodegradation rate constant for ciprofloxacin (1.41×10 2 d −1 ; table SI2.1) is multiple orders of magnitude higher than that of carbamazepine, explaining the difference in environmental removal between the two APIs. We obtained this rate constant as the average of 12 experimentally derived values. Although data were relatively abundant, values ranged from 5.76 d −1 (Belden et al 2007) to 7.85×10 2 d −1 (Babic et al 2013). As a consequence, we might have overestimated the actual environmental photodegradation of ciprofloxacin, providing another possible explanation for the underestimation of measured concentrations.
To further strengthen confidence in the environmental fate calculations, a more extensive validation covering a wider geographical range should be performed, ideally focusing on densely populated and arid ecoregions with relatively high environmental concentrations, as identified in our study.

Conclusions
Time trends of yearly ARs from human pharmaceuticals can be derived with reasonable confidence for all river basins worldwide, by combining spatially explicit chemical fate and effect modelling with predictions of pharmaceutical consumption. To our knowledge, this is the first time such an exercise has been performed for specific APIs at a global scale. Our analysis has shown that ARs due to carbamazepine exposure are typically low, but that risks due to ciprofloxacin exposure resulting from human use are expected to be extensive and widespread. Applying our method to a wider range of (classes of) APIs will provide further insight in the general applicability of our approach.

Disclaimer ECDC
'The views and opinions of the authors expressed herein do not necessarily state or reflect those of ECDC. The accuracy of the authors' statistical analysis and the findings they report are not the responsibility of ECDC. ECDC is not responsible for conclusions or opinions drawn from the data provided. ECDC is not responsible for the correctness of the data and for data management, data merging and data collation after provision of the data. ECDC shall not be held liable for improper or incorrect use of the data'.