A High-Resolution Spatial Model to Predict Exposure to Pharmaceuticals in European Surface Waters: ePiE.

Environmental risk assessment of pharmaceuticals requires the determination of their environmental exposure concentrations. Existing exposure modeling approaches are often computationally demanding, require extensive data collection and processing efforts, have a limited spatial resolution, and have undergone limited evaluation against monitoring data. Here, we present ePiE (exposure to Pharmaceuticals in the Environment), a spatially explicit model calculating concentrations of active pharmaceutical ingredients (APIs) in surface waters across Europe at ∼1 km resolution. ePiE strikes a balance between generating data on exposure at high spatial resolution while having limited computational and data requirements. Comparison of model predictions with measured concentrations of a diverse set of 35 APIs in the river Ouse (UK) and Rhine basins (North West Europe), showed around 95% were within an order of magnitude. Improved predictions were obtained for the river Ouse basin (95% within a factor of 6; 55% within a factor of 2), where reliable consumption data were available and the monitoring study design was coherent with the model outputs. Application of ePiE in a prioritisation exercise for the Ouse basin identified metformin, gabapentin, and acetaminophen as priority when based on predicted exposure concentrations. After incorporation of toxic potency, this changed to desvenlafaxine, loratadine, and hydrocodone.

Introduction

Over the past decades, human consumption of pharmaceuticals has steadily increased.[1,2] In combination with continuing improvements in our analytical capabilities,[3,4] this has led to the detection of many active pharmaceutical ingredients (APIs) in surface waters worldwide.[5,6] The environmental presence of 631 different pharmaceuticals has been reported in 71 countries covering all continents,[5] but the actual number of APIs present in surface waters is likely higher due to the self-fulfilling selection bias of many monitoring campaigns.[7]

A crucial step in the environmental risk assessment of chemicals is the determination of their environmental exposure potential. Since there are currently at least 1500 distinct APIs in use,[8,9] monitoring all of them everywhere and continuously is practically impossible. Moreover, APIs under development will not be present in the environment so monitoring will provide no information on exposure of these molecules. There is therefore a need for exposure modeling approaches that can help us prioritize our monitoring efforts, support more robust environmental risk assessment of new APIs, and that can be used to take targeted measures.[10] These should preferably be spatially explicit, acknowledging that geographical variability can lead to substantial differences in the concentrations of APIs across and within regions.[11,12] For example, rankings of APIs established at the continental European level may lead to misguided allocation of resources when adopted at a regional level.[12] Such mismatches between EU-level and regional level prioritization of APIs might, for example, be the result of geographical variation in API consumption, a heterogeneous distribution of emission sources, or spatially varying environmental conditions driving the fate of APIs after emission.

The environmental exposure potential of chemicals is reflected by the measured (MEC) or predicted (PEC) environmental concentrations at which they occur in the environmental compartment of interest. PECs can be derived using multimedia fate models, such as the EUSES model[13] and our previously developed prioritization tool for APIs.[11] These are based on mass-balance equations for interconnected compartments that represent the relevant environmental media (e.g., fresh and salt waters, air, urban and agricultural soils, et cetera), and are therefore especially useful for larger scale (regional, continental) assessments where multiple media might be relevant. However, they are less suitable for answering locally specific questions (e.g., hotspot identification, scenario analyses for optimal mitigation measures), because they assume a homogeneous distribution of chemicals within their compartments and do not account for any spatial variation at that scale.[14,15] This also inherently limits the options for model corroboration with local measurement data.

APIs tend to largely remain in the compartment where they are emitted,[16] implying that the use of single-media models is also an option. Examples of geographically based single-media models for down-the-drain chemicals are GREAT-ER,[17] PhATE,[18] GWAVA,[19] LF2000-WQX,[20] iSTREEM,[21] and the recent unnamed model by Grill et al.[15] Combined, these models have been applied to assess the distribution of APIs in many river basins worldwide. Invariably, they integrate information on API consumption, human metabolism, removal in wastewater treatment plants (WWTPs), and dilution and dissipation in receiving surface waters, to estimate PECs throughout river basins. The characterization of hydrology is broadly done in one of two ways: via gridded approaches incorporating extensive process-based hydrological models,[15,19] or via segmentation of the river network into discrete river segments with calibration against measured hydrology and extrapolation to ungauged sites.[17,18,20,21] Both approaches have their own drawbacks, related to the computational demands of large scale hydrological models, the extensive data collection and processing efforts required for the parametrization of river basins, and the limited spatial resolution determined by the grid-cell size or the length of individual river segments.

Here, we present ePiE (exposure to Pharmaceuticals in the Environment), a new spatially explicit model, developed in the frame of the Innovative Medicines Initiative iPiE project, that can calculate concentrations of APIs in surface waters throughout river basins in Europe. It is designed to strike a balance between generating data on exposure at high spatial resolution while having limited computational and data requirements. It does so by employing FLO1K for the underlying hydrology, a global geographic data set with annual predictions of streamflow metrics (annual mean flow, highest and lowest monthly mean flow) spatially distributed at 30 arc seconds (∼1 km).[22] This is a resolution 10 times higher than the most detailed global hydrological models or land surface models currently available.[23,24] In ePiE, river networks are represented as collections of interconnected nodes describing emission points, river junctions, river mouths and inlets and outlets of lakes and reservoirs. It thus provides a modeling architecture supporting linkage and integration of geographic information in vector format, i.e., the nodes of the river networks, and rasterized information on climatic, hydrological, and geochemical conditions.[25] We developed a custom routing scheme to follow APIs through the river network, along the way accounting for dissipation from the water via the processes of biodegradation, photolysis, hydrolysis, volatilization and sedimentation.

In this article, we present the structure of ePiE and evaluate its performance against measured concentration data from the open literature for a combined total of 35 APIs in two European river basins. Finally, to illustrate the utility of the model, we apply ePiE to rank APIs in the river Ouse basin (UK), based on predicted concentrations in surface waters and predicted risks to fish.

Materials and Methods

Model Structure

Central to ePiE are a set of network nodes derived from the global databases HydroSHEDS[26] and HydroLAKES,[27] and agglomerations and WWTPs from the UWWTD-Waterbase.[28] This latter database contains information on the location and characteristics (i.e., generated load, design capacity and level of treatment) of 30 043 European urban WWTPs and 27 695 agglomerations with generated wastewater loads above 2000 population equivalents (p.e.). After curation of the UWWTD-Waterbase (see Supporting Information (SI) S1), agglomerations and WWTPs were incorporated into the river network based on their proximity to the nearest water body. Direct emissions into the sea were excluded from the model. Finally, gridded information on air temperature, wind speed, slope, and streamflow was extracted to all nodes in the network. To optimize its flexibility and accessibility, ePiE is entirely constructed in the open-source software environment R,[29] and a description of the model construction can be found in SI S2.

The ePiE model has a modular structure based on the georeferenced river basins provided by the global HydroBASINS database[25] which includes basins below of 60°N. Depending on the river basin of interest, a subset of the total network of nodes is geographically selected. As a starting point, ePiE then requires yearly consumption data for the API of interest (kg/year) for all countries the river basin covers. When the API of interest is formed as a metabolite from another API, that is, its prodrug, consumption data for that prodrug are also needed. Yearly emissions into the river network from WWTPs (Ew,wwtp; kg/year) and from agglomerations with incomplete WWTP connectivity (Ew,agg; kg/year) are calculated via eq and eq , respectively. The country-specific yearly consumption data (M) include the prescription of pharmaceuticals in hospitals. This means that hospital emissions are not included as location-specific point sources, but spatially distributed according to the wastewater loads per agglomeration (i.e., a proxy for population density).where M and Mpd are the yearly consumption of the API of interest and its prodrug in the relevant country (kg/year); fpc is the fraction of the administered parent compound excreted/egested unchanged or as reversible conjugates via urine and faeces (−); fmet is the fraction of prodrug metabolized to the API of interest, and subsequently excreted/egested via urine and faeces (−); n is the number of agglomerations j connected to the WWTP (−); fconn,agg, is the level of WWTP connectivity per agglomeration j; fwwtp,agg, is the fraction of agglomeration j connected to the WWTP; frem is the API-specific removal efficiency per WWTP (−); and Vww,agg, and Vww,cnt are the wastewater loads generated per agglomeration j and the total in the relevant country, respectively (p.e.).

The SimpleTreat 4.0 model[30,31] was incorporated into ePiE to estimate the removal efficiency during wastewater treatment (frem). It requires basic physicochemical properties as input, as well as solids-water partitioning coefficients for primary sewage (Kpps; L/kg) and activated sludge (Kpas; L/kg), and (pseudo)first order biodegradation rate constants (kbio,wwtp; s–1). Removal efficiencies were assigned to individual WWTPs depending on their associated level of treatment, using either the full SimpleTreat 4.0 model for those employing consecutive primary and secondary treatment, or the module for primary treatment only.

After their emission, API residues are followed through the river network using a routing procedure ordered from the most upstream to the most downstream nodes. As such, the contribution of all upstream emissions to local concentrations is considered. Along the way, ePiE accounts for dilution in the water column and five (pseudo)first order loss processes, three being degradation processes, that is, biodegradation, photolysis and hydrolysis, and two being intermedia transport processes, that is, sedimentation and volatilization. Eq calculates concentration C (μg/L) at any node i in the river network; eq calculates concentrations in lakes and reservoirs, following an approach similar to Grill et al.[15] in which they are modeled as single completely stirred tank reactors.where Ew, and Ew, are the emissions into the river network at node i and at node j upstream from node i, respectively (mg/s); n is the total number of nodes upstream from node i (−); d is the distance over the river network between node j and node i (m); k is the average (pseudo-) first order rate constant for loss process m over d (s–1); v is the average river flow velocity over d (m/s); and Q is the total river flow at node i (m3/s), including any discharges.where Ew, is the emission into lake or reservoir i coming from node p (mg/s), which can either be a direct emission source (i.e., a WWTP or an agglomeration), or an inlet point carrying API residues from upstream the river network; n is the total number of nodes emitting into lake or reservoir i (−); HRT is the hydraulic retention time of lake or reservoir i (s); V is the volume in lake or reservoir i (m3); and k is the (pseudo-) first order rate constant for loss process m in lake or reservoir i (s–1).

Individual loss rate constants are extrapolated from test to field conditions by accounting for temperature differences, sorption to suspended solids and dissolved organic carbon,[32] and reduced light intensity.[33] Local sedimentation and volatilization rate constants are implemented via mass transport velocities between media.[34] Detailed information on the extrapolation to field conditions can be found in SI S3.

For characterization of annual mean flow, and highest and lowest monthly mean flow, the recent global FLO1K data set was implemented in ePiE.[22] FLO1K is based on an ensemble of artificial neural networks regressions, with upstream-catchment physiography (area, slope, elevation) and year-specific climatic variables (precipitation, temperature, potential evapotranspiration, aridity index and seasonality indices) as covariates. It provides estimations of flow at a spatial resolution of 30 arc seconds (∼1 km) for the years 1960–2015, which are in good agreement with independent data (global R2 of single-year metrics up to 0.91). An additional comparison with independent data obtained from 1,007 European monitoring stations for the period 2010–2015,[35] showed that year-specific annual mean flow, and highest and lowest mean monthly flow in European rivers are predicted well, with R2 values of 0.97, 0.95, and 0.91, respectively (Figure ).

Figure 1

Validation results for year-specific annual mean flow (A), highest monthly mean flow (B) and lowest monthly mean flow (C). Independent validation data set consisted of yearly measurements (2010–2015) from 1,007 GRDC European stations. The solid line represents perfect model fit (1:1 line) and the dashed lines represent a difference of 1 order of magnitude.

Additional hydrological parameters flow velocity v (m/s) and river depth h (m), were calculated via the Manning’s equation for open channel flow, rewritten under the assumption of a wide rectangular river cross section as proposed by Pistocchi and Pennington.[36] In this approach, river width was related to river flow using their power law equation for European rivers (R2 of 0.87).[36]

Model Evaluation

We performed a model evaluation exercise with measured concentrations for 35 APIs consumed in Europe and covering a wide range of pharmaceutical classes. Excretion, sorption and degradation data were extracted from open literature by cross-referencing a set of reviews on human metabolism, sludge sorption, sediment sorption, biodegradation and photolysis. The data obtained were supplemented with additional API-specific searches. The resulting data set was extensive, containing a total of 430 sorption coefficients and 342 degradation rate constants, but not homogeneously distributed over the 35 APIs. Complete experimental data sets were available for 13 APIs, while 12 were missing data on at least one sorption process and 11 on at least one degradation process. No experimental sorption or degradation data were found for sitagliptin and triamterene. Missing sorption coefficients were substituted by combining default mass fractions of organic carbon for sludge[30] or sediments[37] with QSAR predictions of organic carbon–water partition coefficients.[38,39] Moreover, if only ready biodegradability screening test data were available, APIs were assigned a biodegradation rate constant as proposed by Jager et al.[40] When experimental degradation rate constants were lacking altogether, no degradation was assumed. SI Tables S4.1 and S4.2 show the physicochemical and environmental fate properties of the 35 APIs, respectively.

Predicted environmental concentrations were compared with measured concentrations extracted from a database compiled by the German national environmental protection agency,[5] and a limited number of more recent literature studies. Individual studies were included in the model evaluation if (1) measurements were performed after 2010, (2) measurement locations were provided, (3) at least 10 of our APIs were measured above their limit of detection at least 10% of the time, and (4) multiple consecutive measurements were performed over time. These criteria resulted in the selection of three literature studies, being those by Burns et al.,[41] who measured APIs in the river Ouse basin in the United Kingdom, and by Ruff et al.[42] and Munz et al.,[43] who both measured APIs in the river Rhine basin in Northwestern Europe (Figure ). Burns et al.[41] included a total of 30 of our preselected APIs in a monthly grab-sampling campaign throughout 2016. They reported the coordinates of their 11 sampling locations, of which six were located along the river Ouse and five along its tributary, the river Foss, and we integrated these as such into ePiE. The yearly average of the Burns et al.[41] data set was compared to the PEC obtained under annual mean flow conditions for 2015. Ruff et al.[42] measured a total of 23 of our preselected APIs in a weekly flow-proportional composite sampling campaign during “a remarkably dry period with constant low flow conditions” in the early spring of 2011. To reflect these low flow conditions, we used PECs derived under lowest monthly mean flow for 2011 in the quantitative evaluation of model performance. Out of their 16 sampling locations, ten were sampling stations along the river Rhine, but their coordinates were not reported. We georeferenced these sampling locations based on the proximity of the cities mentioned by the authors to sampling stations in the GRDC Station Catalogue.[35] In addition, they sampled six tributaries of the river Rhine. We assumed these were sampled directly before their confluence with the main river. Finally, Munz et al.[43] included a total of 11 of our preselected APIs in two distinct grab-sampling campaigns in 2013 and 2014. Their 24 sampling locations were split evenly over these two campaigns and were all located directly downstream of WWTPs in Switzerland. Two sampling locations outside the river Rhine basin were excluded from our model evaluation. Similar to Ruff et al.,[42] Munz et al.[43] explicitly chose their sampling times to capture low flow conditions. Therefore, we used PECs derived under lowest monthly mean flow conditions for 2013 (site 1–12) and 2014 (site 13–24).

Figure 2

Overview of studies included in the model evaluation exercise, with numbered sampling locations from Burns et al.[41] (A), Ruff et al.[42] (B), and Munz et al.[43] (C).

For estimations in the river Ouse basin, we used consumption data for 2016 from the Prescription Cost Analysis.[44] For the river Rhine basin, consumption data for The Netherlands were obtained from the Dutch National Health Care Institute.[45] German, French and Swiss consumptions during the years of interest were mostly extrapolated from per capita consumption in other years.[46] Consumption data were not available for five APIs in France, one API in Switzerland, and all APIs in Austria, Belgium, and Luxembourg. In these cases, we averaged the per capita consumption from the basin’s other countries. All consumption data are presented in SI S5.

To assess the predictive accuracy of ePiE, we computed the median symmetric accuracy ξ per study included in the evaluation exercise (eq ).[47] This metric reflects the typical percentage error of the predictions compared to the measurements. For example, a ξ of 100% indicates that predicted concentrations will typically be within a factor of 2 of the measurements. Contrary to metrics based on scale-dependent errors (e.g., root-mean-square error RMSE), ξ assigns equal importance to deviations of the same order rather than the same magnitude. This is especially relevant for our data where concentrations ranged from low ng/L to μg/L levels. In other words, a situation where the PEC is 1 ng/L and the MEC is 10 ng/L (absolute error 9 ng/L) receives an equal penalty to that where the PEC is 100 ng/L and the MEC is 1 μg/L (absolute error 900 ng/L). Moreover, since ξ bases on the median (M) of the accuracy ratios of individual pairs of predictions and measurements, it penalizes under- and overpredictions equally. This is an advantage over the often-applied mean absolute percentage error MAPE, which penalizes overpredictions more heavily.[47]

Additionally, we assessed the prediction bias of ePiE by computing the symmetric signed percentage bias (SSPB) (eq ), which is closely related to the median symmetric accuracy ξ.[47] The SSPB can be interpreted similarly to a mean percentage error, but is not affected by the likely asymmetry in the distribution of percentage error.

Model Application

To illustrate the utility of the model, we applied ePiE to prioritise APIs in the Ouse river basin, the basin with the best model performance and most APIs included. Additional nodes were integrated into the network at evenly spaced one-kilometre distances, enabling a basin-wide prioritisation using geographically homogeneous aggregate statistics. In addition to a ranking based on concentrations, we ranked the APIs based on their potential risks to fish. For this we followed a similar method as Burns et al.,[48] based on the fish plasma model approach.[49,50] We extrapolated concentrations in surface water to concentrations in fish plasma using bioconcentration factors computed according to Fitzsimmons et al.[51] for neutral compounds, and Fu et al.[52] for ionizing compounds. The latter were derived assuming a surface water pH of 7.4.[53] Risk quotients (RQ) for fish were then calculated as the ratio of concentrations in fish plasma over therapeutic concentrations in human plasma, which we obtained from the MaPPFAST database.[54] A risk quotient exceeding 1 thus indicates that the concentration of an API in surface water is expected to cause a pharmacological effect in fish, assuming equivalent pharmacological activity as in humans.[55] Finally, to enable exploration of local concentration and risk patterns, model results were geographically visualized as interactive html-maps, using the leaflet package “leafletR” in the R environment.[56]

Results and Discussion

Out of the 940 predicted values used for model evaluation, 36% were qualified as nondetects in the measurement campaign. We qualified a substance as a nondetect in case it was below the limit of detection (LOD) in at least 40% of the samples taken at that location. Such nondetects are less suitable for a quantitative evaluation of model performance. We did, however, include them in a binary comparison between predicted min-max concentration ranges, resulting from the temporal variation in flow conditions, and measurements in relation to their LOD (Figure ). Assigning comparisons to one of four bins (detected, predicted < LOD; not detected, predicted > LOD; detected, predicted > LOD; not detected, predicted < LOD), there was 94%, 88%, and 90% coherence of predictions and measurements for the Burns et al.,[41] Ruff et al.,[42] and Munz et al.[43] studies, respectively (green bars in Figure ).

Figure 3

Binary comparison of measurements and min-max range of predictions, relative to their limit of detection (LOD). All combinations of location and API from Burns et al.[41] (A), Ruff et al.[42] (B), Munz et al.[43] (C), and all studies combined (D).

For a quantitative assessment of model performance, we included all detects at locations downstream of a WWTP, that is, for which PEC > 0. In case measured values were below the LOD (i.e., always less than 40%), these measurements were replaced by .[48] The resulting comparison of predicted versus measured values (Figure ) revealed a substantial variation between the three studies. Model accuracy was best for predictions in the Ouse river basin, with a typical percentage error of 86% (Figure A; Burns et al.[41]). Predictions in the river Rhine basin had typical percentage errors of 143% (Figure B; Ruff et al.[42]) and 158% (Figure C; Munz et al.[43]). Model performance was similar if data points were included for which PEC > LOD and for which more than 40% of the measurements were below the LOD (SI Figure S6.1).

Figure 4

Predicted concentrations (i.e., > 0) versus detects (i.e., < 40% of the measurements below LOD), separately for data from Burns et al.[41] (purple; A), Ruff et al.[42] (golden; B), Munz et al.[43] (green; C), and for all studies combined (black; D). Concentrations predicted under annual mean flow conditions (A) or lowest monthly mean flow conditions (B and C). Solid line represents 1:1 relationship; dashed lines represent 1:10 and 10:1 relationships. ξ: median symmetric accuracy; SSPB: symmetric signed percentage bias.

The worse performance of ePiE in the river Rhine basin might relate to the quality of the consumption data used in the calculations. First, Swiss and German consumption data were often reported as “greater-than” values instead of exact amounts.[46] Second, we extrapolated the consumption in 2009 to that in the actual years of sampling (2011–2014), based on changing demographics and the assumption of a constant per capita consumption over the years (SI Table S5.1). However, actual per capita consumption has increased significantly for at least some pharmaceuticals, e.g., antidiabetics like sitagliptin[57] or antidepressants like venlafaxine.[58] These were therefore underestimated by ePiE due to the temporal extrapolation. In addition, errors might have been introduced when sampling sites from Ruff et al.[42] were allocated to the river network, because limited geographical detail was available on their specific locations. Inaccuracies may also be due to the fact that HydroSHEDS does not provide the real geometry of a river network in a basin, but most likely flow paths between individual cells according to flow accumulation. Similarly, errors might have been introduced during the allocation to the river network of the WWTPs sampled by Munz et al.[43] These were all located at smaller streams in the upper Swiss catchment of the Rhine river basin, without other upstream emission sources. In such smaller upstream catchments, proximity-based allocation is more prone to errors because the main stream within the floodplain is less easily identified. Nevertheless, the ξ values and the scatterplots in Figure indicate that concentrations were typically predicted within a factor of 2–3, with approximately 95% of predictions within a factor of 10.

Concentrations measured by Burns et al.[41] were typically underestimated by ePiE, with a symmetric signed percentage bias (SSPB) of −44% (Figure A). From the scatterplot in Figure A, underestimations seem to be more prominent at lower concentrations. This can at least partly be explained by the fact that measured concentrations have a lower bound in the form of their LOD, while model predictions do not. As a consequence, underestimations are more likely than overestimations in the vicinity of that LOD, since nondetects are excluded from the comparison. Indeed, model performance slightly improved if data points were included for which PEC > LOD, and which had more than 40% of the measurements below the LOD which were replaced by (SI Figure S6.1). Additionally, the reliability of measured concentrations decreases closer to the LOD. This complicates the evaluation of model performance, because any difference between predicted and measured concentrations might then be attributed to errors in either of them. Finally, inputs from tourism, specific point sources (e.g., hospitals), operation of combined sewer overflows at selected times of the year and use of over the counter medicines may also explain the slight mismatch between measurements and predictions in the river Ouse basin.

In contrast to the river Ouse basin, concentrations measured in the river Rhine basin were typically slightly overestimated, with SSPB values of 30% and 5% (Figures B and 4C). When we ran ePiE under annual mean flow settings, these values dropped considerably to −70% and −313%, respectively. This indicates that actual streamflow during sampling was probably somewhere between lowest monthly mean flow and annual mean flow conditions.

Ratios of predicted over measured concentrations (PEC/MEC ratios) provide further insights into the performance of ePiE (Figure ). PEC/MEC ratios are grouped according to study and sampling location, numbered as in Figure . Similar graphs grouped according to API are included in the Supporting Information (SI Figure S6.2). Figure A shows that the spread around predictions in the river Ouse (locations 1–6) is smaller than around those in its tributary river Foss (locations 7–11). This indicates that ePiE predicts concentrations in larger rivers better than in smaller ones. While concentrations in larger rivers reflect an accumulation of APIs over a larger upstream catchment area, concentrations in smaller rivers and streams are more directly influenced by specific local conditions, that is, water extraction and retention or small scale discharges. Indeed, comparison of predicted and measured mean annual flow at two gauging stations, i.e. one in the river Ouse and one in the river Foss (SI Table S6.1), shows that our flow prediction is less accurate for the smaller river Foss. The impact of local conditions can furthermore be observed at the most upstream location on the river Foss (location 7), where multiple APIs were detected but ePiE predicted zero concentrations for all of them. This deviation was likely due to the presence of a small upstream WWTP not included in the UWWTD-Waterbase because its size was below the reporting threshold of 2,000 p.e. National consumption data and default WWTP characteristics might thus not always suffice to estimate concentrations in locally influenced rivers. The same likely holds for the tributaries of the river Rhine sampled by Ruff et al.[42] (locations 11–16) and by Munz et al.[43] However, the pattern is less obvious here, probably due to errors introduced by the aforementioned incoherent flow conditions, consumption data, and geographical detail on sampling locations and emission sources. One option to improve predictions in upstream tributaries is to extend the UWWTD-Waterbase with WWTPs smaller than 2000 p.e.

Figure 5

Ratios of predicted over measured concentrations (PEC/MEC), reported by Burns et al.,[41] (A), Ruff et al.,[42] (B) and Munz et al.[43] (C). Colored dots are individual combinations of API and location, measured above the LOD; black bars represent 95th percentile and median over all measurements per location (numbered as in Figure ). Concentrations predicted under annual mean flow conditions (A) or lowest monthly mean flow conditions (B and C). * = The PEC/MEC ratio of location 7 in panel A equals zero.

Figure A shows that predicted concentrations in the river Ouse basin were highest for metformin, gabapentin and acetaminophen, mainly resulting from their large consumption volumes, high excretion fractions and/or relatively poor degradation (SI S4.2 and S5). The prioritisation of APIs shifts when based on potential risks to fish instead of concentrations (Figure B). Metformin, gabapentin and acetaminophen drop down the list and are replaced by other more pharmacologically active APIs. Desvenlafaxine, loratadine, and hydrocodone (highlighted in Figure A) then become APIs of particular interest. Their risk quotients for fish were larger than 0.1 in one or more locations in the river basin, with risk quotients for desvenlafaxine and loratadine even exceeding 1 in ∼26% and ∼10% of the river length, respectively. Interestingly, desvenlafaxine is formed as a metabolite of its prodrug venlafaxine but is not administered as a separate medication in the United Kingdom. This provides a strong argument for more focus on active metabolites in the environmental risk assessment of pharmaceuticals. Finally, Figure shows that higher risks are mainly found in more densely populated areas, for example, around the city of Leeds. The geographical distribution of surface water concentrations and risk quotients for all APIs is visualized in interactive html-maps in SI S7.

Figure 6

Ranking of all APIs modeled with ePiE in the Ouse river basin, based on concentrations (A) and risk quotients for fish (B) predicted throughout the river basin, excluding zero concentrations. Boxes indicate interquartile range including median; whiskers indicate 1st–99th percentile range for the total river length. Red boxes: RQ exceeds 1 at least somewhere in the river network; amber boxes: RQ exceeds 0.1 at least somewhere in the river network; green boxes: RQ below 0.1 throughout the river network.

Figure 7

Spatial distribution of risk quotients for the three top ranked APIs in the river Ouse basin (UK): desvenlafaxine (A), loratadine (B), and hydrocodone (C). Panel D depicts the spatial variation in population density in the river Ouse basin (individuals/100 m2).

Our model evaluation showed that ePiE generally predicts concentrations in surface waters within 1 order of magnitude of measured concentrations for a wide range of pharmaceutical classes. While other models have been shown to predict PECs of APIs to within a factor 2–15 of measured concentrations,[59] none of these models have been evaluated using such an extensive data set on a diverse range of APIs. To further strengthen confidence in the model, future model development and evaluation should extend toward additional, more hydrologically and climatically diverse river basins. As part of the IMI funded project iPiE, we are currently monitoring additional river basins in Denmark, Germany, Spain, and the UK to develop a broader data set against which to evaluate the model. Because of its flexible setup and the use of global high-resolution gridded streamflow,[22] ePiE can be extended to new basins worldwide in a relatively straightforward way. Our model results also showed that a proper assessment of model performance requires measured concentrations derived under the same conditions as those modeled. This means that further model development should ideally be supported by long-term annual sampling efforts. In addition, incorporation of local consumption patterns, point sources (e.g., hospitals and pharmaceutical production plants), WWTP characteristics, and environmental conditions, would be especially relevant for adequate estimation of concentrations in smaller river stretches.