Bayesian-Based Probabilistic Risk Assessment of Fipronil in Food: A Case Study in Taiwan

Fipronil, a broad-spectrum insecticide, is widely used in agriculture and veterinary practices. Fipronil-induced neurotoxicity and potential adverse effects on humans and aquatic organisms have raised health concerns. Monitoring programs have been implemented globally to assess fipronil residues in food, including fruits, vegetables, and animal products. However, previous exposure assessments have often focused on specific food categories or subsets of items, resulting in limited insights into the overall health risks. Additionally, the large number of non-detect fipronil residues in food has introduced uncertainties in exposure assessment. To address these issues, a probabilistic exposure assessment and dose-response analysis were adopted in this study, considering the sample distribution below the detection limit to better characterize uncertainties and population variability in health risk assessments. The estimated fipronil exposure to the general public ranges from 6.38 × 10−6 ± 0.00017 mg/kg/day to 9.83 × 10−6 ± 0.00034 mg/kg/day. Only one out of 200,000 simulated individuals had a fipronil dose exceeding the probabilistic reference dose (0.048 mg/kg/day, pRfD), which aims to protect 99% of the population with effects less than 10% extra risk. By incorporating uncertainties in exposure and dose-response data, a more comprehensive understanding of the health risks associated with fipronil exposure in the Taiwanese population has been achieved.


Introduction
Fipronil, a broad-spectrum phenylpyrazole insecticide, has been widely used in crop protection, pest control, and veterinary practice. Fipronil shows potent neurotoxicity to insects and aquatic organisms such as Daphnia magna [1], shrimps [2], and crustaceans [3] via the inhibition of the gamma-aminobutyric acid (GABA)-gated chloride channel, resulting in convulsions, paralysis, and death [4]. Fipronil exposure can induce oxidative stress, neurotoxicity, hepatotoxicity, and nephrotoxicity in rodents and humans [5][6][7]. Additionally, studies showed that maternal exposure to fipronil could be placentally transferred to the fetus, leading to adverse health effects in newborns [8]. The widespread agricultural uses and the toxic potential of fipronil have raised increasing public concerns about fipronil residues in food.
National regulatory agencies have routinely surveyed fipronil residues in food to ensure the consumer's safety. The fipronil monitoring data in food have been reported in China [9], Europe Union [10], India [11], Japan [12], and Taiwan [13]. Aside from fruits and vegetables, chicken eggs and animal products have been included in the surveillance programs to monitor the misuse of non-approved veterinary drugs, such as fipronil, in poultry farms against red mites [14,15]. Fipronil contamination in eggs has been reported in Europe, South Korea, Hong Kong, and Taiwan since 2017 [16,17]. Even though the fipronil We reconstructed the fipronil concentrations in food based on a modified Markov Chain Monte Carlo (MCMC) method as described by Suzuki et al. [30], using the "rstan" package (version 2.21.2) in R software (version 4.0.4). The R code is detailed in Supplementary Materials File S1. Briefly, the Stan model has 5 components: "data", "transformed data", "parameters", "model", and "generated quantities". The "data component" specifies the number of observed samples (N_obs), the number of censored samples (N_cen), and the value of observed data (Y_obs). The "transformed data" component computes the minimum and maximum values of Y_obs. The "parameters" component defines geometric mean (GM), geometric standard deviation (GSD), reporting limit (RL), and the estimates for non-detected values (Y_cen).
The "model" component specifies the prior distributions for the declared parameters, which plays an important role in MCMC modeling. Among the detectable samples (n = 120), fipronil residue in food follows the lognormal distribution as determined by the Kolmogorov-Smirnov test (α = 0.05). The geometric standard deviation of the detected fipronil concentrations was 4.16. Therefore, we applied the lognormal distribution with GM of 4 and GSD of 2 for the prior distribution of GSD. For the food category with only one detectable sample, we added an extra detectable value at the detection limit (0.002 ppm) to successfully run the MCMC modeling. The normal distributions were used as priors for GM~(Y_obs_max/2, Y_obs_max/4) and RL~(Y_obs_min, Y_obs_min/5).
The "generated quantities" component calculated the posterior distribution of log(arithmetic) density for observed and censored data. The cumulative distribution function (CDF) was used to estimate the likelihood of censored data. Lastly, the arithmetic mean and standard deviation were calculated from the GM and GSD.
For the MCMC runs, we calculated 4 Markov chains with 25,000~50,000 iterations per chain. The ratio of inter-chain variance to intra-chain variance (R) less than 1.1 was used to determine the convergence of 4 MCMC chains. For the MCMC runs reaching convergence, the mean distribution parameters (mean_est, std_est) from the last 50% of iterations were used to define the posterior lognormal distributions of fipronil residues in food.

Dietary Exposure Assessment of Fipronil in Taiwan
The dietary exposure assessment of fipronil was achieved by simulating the dietary profiles of 200,000 individuals in Taiwan. The age-stratified (age 0-100) and sex-stratified demographic data were obtained from the Department of Household Registration, Ministry of Interior (data retrieved from April 2023). The body weight's age-and sex-stratified data (mean and standard deviation) was acquired from the Nutrition and Health Survey in Taiwan (NAHSIT, 2013-2016). The simulated number of 200,000 was set to meet the minimal sample size of 10 for each subgroup (per sex per age). The food intakes of 11 food categories (Table 1) were organized from the National Food Consumption Database [31], a four-level food categorization system that provides food consumption estimates for the general population or consumer-only.
Truncated normal distribution was used to randomly assign the body weight and food intakes (Food intake i ) for the simulated male (n = 98,788) and female (n = 101,212) individuals. Fipronil residues in food (Residue i ) were stochastically generated using the transformed posterior mean and standard deviation (Formulas (1) and (2)), following lognormal distribution.
where Meanlog and SDlog are the mean and standard deviation of the distribution on the log scale. Finally, the aggregate exposure of fipronil was calculated using the following Equation (3): Aggregate fipronil exposure mg kg where a-k denotes the 11 food categories, j = 1,. . ., The lower bound (LB, general population) and upper bound (UB, consumer only) estimates were reported accordingly.

Revisit the Health Reference Dose of Fipronil Using Bayesian Benchmark Dose Modeling
The probabilistic reference dose (pRfD) of fipronil was derived using a web-based system, Bayesian BenchMark Dose modeling (BBMD, Available online: https://benchmarkdose.com/ accessed on 4 July 2023) [28]. The chronic animal data of fipronil were retrieved from the toxicological evaluation report for pesticide residues in food in 2021 [32]. Convulsions in male SD rats were deemed as the most sensitive endpoint and subject to probabilistic doseresponse assessment (dose: 0, 0.019, 0.059, 1.27, and 12.68 mg/kg/day; incidence/total animals: 0/50, 0/50, 3/50, 1/50, and 5/50). The default models for dichotomous data (Logistic, Loglogistic, Probit, Logprobit, Quantal linear, Multistage (2nd order), Weibull, and Dichotomous Hill) and non-informative prior were used to fit the dose-response data, using Markov chain Monte Carlo (MCMC) simulation. Three Markov chains were calculated with 50,000 iterations per chain. The model distributions were estimated using the last 75,000 iterations (25,000 iterations per chain). The random seed was set at 76,316. From the fitted dose-response curves, BMD distributions were acquired using 10% extra risk as the benchmark response (BMR). The model average from 8 dichotomous models using the posterior model weights was used to establish the BMD distribution.
Next, BMD in the rat (BMD rat ) were converted to BMD in human (BMD h ), using the following parameters: human body weight, 70 kg; allometric scaling exponent (mean), 0.7; allometric scaling exponent (standard deviation), 0.0243. We applied probabilistic distribution to describe the uncertainty in inter-species (geometric mean (GM) = 1; geometric standard deviation (GSD) = 1.95) and intra-species (GM = 0.746; GSD = 1.5935) extrapolation. The estimated human dose, where 50% of the population has effects greater than or equal to 10% extra risk (HD 0.5 ), the estimated human dose at which 1% population has effects greater than or equal to 10% extra risk (HD I=0.01 M=0.1 ), and probabilistic reference dose (pRfD, the 5th percentile HD I M ) were reported accordingly.

Risk Characterization of Fipronil Exposure in Taiwan
The fipronil-induced health risk in Taiwan was characterized using a probabilistic hazard quotient (HQ) approach. Briefly, the HQ or probabilistic HQ (pHQ) was calculated by dividing the aggregate fipronil exposure by the ADI or the probabilistic reference dose (pRfD, the 5th percentile of HD I=0.05 M=0.1 ): (p)HQ =

Statistical Analysis and Data Visualization
The MCMC runs and the simulation for the Taiwanese population were performed in R software (version 4.0.4, R Development Core Team, Vienna, Austria), using packages "rstan (version 2.21.2)" and "truncnorm (version 1.0-9)". The bar graphs, line graphs, and scatter plots were created using GraphPad Prism 9 (version 9.5.0).

Fipronil Residues in Food
Among the investigated food categories, the detection rate of fipronil was 0.05% for LF, 0.1% for CF, 0.11% for SF, 0.24% for CB, 0.35% for FV, 0.43% for LV, 0.47% for SP, 0.6% for SV, 2.18% for DB, 2.5% for RV, and 5% for CG (chicken egg, CG, Table 1). The detection limit was 0.002 ppm. The highest concentration was found in mustard green (1.69 ppm). The root and stem vegetables and chicken eggs are likely the major contributors to total fipronil exposure, considering both the intake and detection rates of fipronil in food ( Figure 1).

Figure 1.
Comparative analysis of food intake and detection rate for each food category.

Reconstruction of fipronil residues using MCMC simulation
Using MCMC simulation, we reconstructed the lognormal distribution of fipronil residues in food. Figure 2 illustrates the reconstruction of fipronil residue distribution in RV. The fipronil concentrations of the samples below the detection limit (n = 1716) were replaced by the random values generated from the posterior distribution ( Figure 2A). This approach kept the same proportion of detectable samples in the original dataset (2.5%) yet projected fipronil concentrations in non-detects, following the lognormal distribution ( Figure 2B).

Reconstruction of Fipronil Residues Using MCMC Simulation
Using MCMC simulation, we reconstructed the lognormal distribution of fipronil residues in food. Figure 2 illustrates the reconstruction of fipronil residue distribution in RV. The fipronil concentrations of the samples below the detection limit (n = 1716) were replaced by the random values generated from the posterior distribution ( Figure 2A). This approach kept the same proportion of detectable samples in the original dataset (2.5%) yet projected fipronil concentrations in non-detects, following the lognormal distribution ( Figure 2B).

Reconstruction of fipronil residues using MCMC simulation
Using MCMC simulation, we reconstructed the lognormal distribution of fipronil residues in food. Figure 2 illustrates the reconstruction of fipronil residue distribution in RV. The fipronil concentrations of the samples below the detection limit (n = 1716) were replaced by the random values generated from the posterior distribution (Figure 2A). This approach kept the same proportion of detectable samples in the original dataset (2.5%) yet projected fipronil concentrations in non-detects, following the lognormal distribution ( Figure 2B).
x FOR PEER REVIEW 7 of 13

Probabilistic risk characterization of fipronil in Taiwan
Next, the pRfD was used to characterize the chronic health risk caused by fipronil, adopting a hazard quotient approach ( Figure 5). Regardless of LB or UB exposure scenarios, only one out of the 200,000 simulated individuals reported fipronil exposure exceeding the pRfD value (i.e., HQ > 1). Overall, younger individuals receive greater risks than the elderly.

Probabilistic risk characterization of fipronil in Taiwan
Next, the pRfD was used to characterize the chronic health risk caused by fipronil, adopting a hazard quotient approach ( Figure 5). Regardless of LB or UB exposure scenarios, only one out of the 200,000 simulated individuals reported fipronil exposure exceeding the pRfD value (i.e., HQ > 1). Overall, younger individuals receive greater risks than the elderly.

Discussion
The importance of addressing uncertainty and variability in exposure and dose-response assessments has been highlighted by the World Health Organization International Programme on Chemical Safety (WHO/IPCS) since 2008 [33][34][35]. Current practices in fipronil risk assessment have often considered the uncertainty and variability in exposure assessments [20,36,37], yet there is no study performing probabilistic dose-response analysis for the risk characterization. To the best of our knowledge, this study is the first to characterize the fipronil-induced health risk by adopting both probabilistic exposure assessment and probabilistic dose-response evaluation. This study brings several refinements and essential insights into the dietary exposure assessment, probabilistic dose-response analysis, and risk characterization of fipronil in Taiwan.
We significantly reduced the uncertainty in the dietary exposure assessment of fipronil. The uncertainty in exposure assessment usually results from an insufficient understanding of relevant exposure scenarios, exposure models, and model inputs [33]. Regarding dietary exposure assessment, uncertainties can arise from food consumption estimates or chemical residues in food. The uncertainty associated with food consumption was minimized via conducting lower bound (i.e., the general public) or upper bound (i.e., consumer only) exposure estimations. However, the actual concentrations of non-detects can range from zero to the analytical detection limit (DL), with which the uncertainty is difficult to address. Investigators often tackle the non-detects using a conservative approach: substituting the non-detects with 2 , √2 , or DL, yet simply substituting values for non-detects has been shown to cause biased results [38]. Alternative approaches include robust regression on order statistics (i.e., robust ROS), maximum likelihood estimation (MLE), and Kaplan-Meier (KM) method [39]. Most methods can provide a reliable estimation when the censoring rate is less than 50% but are inapplicable when censoring rates are greater than 80%. Bayesian inference has been increasingly used to reconstruct the sample distribution for datasets with high censoring rates (>80%) in the context of risk assessment [30,40,41]. In this study, the censoring rates of fipronil concentrations (i.e., below DL) in food are generally greater than 95%. Therefore, we adopted the Bayesian MCMC approach to reduce the uncertainty deriving from the non-detects.

Discussion
The importance of addressing uncertainty and variability in exposure and doseresponse assessments has been highlighted by the World Health Organization International Programme on Chemical Safety (WHO/IPCS) since 2008 [33][34][35]. Current practices in fipronil risk assessment have often considered the uncertainty and variability in exposure assessments [20,36,37], yet there is no study performing probabilistic dose-response analysis for the risk characterization. To the best of our knowledge, this study is the first to characterize the fipronil-induced health risk by adopting both probabilistic exposure assessment and probabilistic dose-response evaluation. This study brings several refinements and essential insights into the dietary exposure assessment, probabilistic dose-response analysis, and risk characterization of fipronil in Taiwan.
We significantly reduced the uncertainty in the dietary exposure assessment of fipronil. The uncertainty in exposure assessment usually results from an insufficient understanding of relevant exposure scenarios, exposure models, and model inputs [33]. Regarding dietary exposure assessment, uncertainties can arise from food consumption estimates or chemical residues in food. The uncertainty associated with food consumption was minimized via conducting lower bound (i.e., the general public) or upper bound (i.e., consumer only) exposure estimations. However, the actual concentrations of non-detects can range from zero to the analytical detection limit (DL), with which the uncertainty is difficult to address. Investigators often tackle the non-detects using a conservative approach: substituting the non-detects with DL 2 , DL √ 2 , or DL, yet simply substituting values for nondetects has been shown to cause biased results [38]. Alternative approaches include robust regression on order statistics (i.e., robust ROS), maximum likelihood estimation (MLE), and Kaplan-Meier (KM) method [39]. Most methods can provide a reliable estimation when the censoring rate is less than 50% but are inapplicable when censoring rates are greater than 80%. Bayesian inference has been increasingly used to reconstruct the sample distribution for datasets with high censoring rates (>80%) in the context of risk assessment [30,40,41]. In this study, the censoring rates of fipronil concentrations (i.e., below DL) in food are generally greater than 95%. Therefore, we adopted the Bayesian MCMC approach to reduce the uncertainty deriving from the non-detects.
The selection of prior distribution in MCMC modeling could significantly influence the results of the posterior distribution. Empirically, the chemical residues in food are often assumed to follow a normal or lognormal distribution. In this study, we conducted normality tests (i.e., Kolmogorov-Smirnov and Shapiro-Wilk test) for the fipronil concentrations in fruits and vegetables (n = 118), confirming that the data follow lognormal distribution rather than the normal distribution. Our modeled data fit very well with the original data, as determined by the cumulative and relative frequency plots ( Figure 2B). Compared with the simple substitution method, the sample distribution reconstructed using MCMC does largely decrease the uncertainty associated with non-detects.
Herein, we refined the dose-response analysis to derive the exposure limit of fipronil. Probabilistic dose-response assessment is a relatively new concept incorporated only in a few studies [25,29,42,43]. The unified framework for probabilistic dose-response assessment was developed in 2015, where the lower bound fifth percentile estimate of "target human dose" (HD I M ), instead of the traditional RfD or ADI, was proposed to serve as the exposure limit for risk characterization [23]. Target human dose requires the declaration of the magnitude of critical effect M and target population incidence I. Instead of applying deterministic uncertainty factors in the conventional approach, HD I M incorporates probabilistic factors for interspecies body weight scaling, interspecies toxicokinetic and toxicodynamic differences, and human variability for population incidence I. The lower bound HD I M estimate, or pRfD, is scientifically rigorous and transparent, informing risk management for different decision contexts [25]. In this study, HD I=0.01 M=0.1 was calculated to represent the dose at which 1% population has effects greater than or equal to 10% extra risk. In other words, we aimed to protect 99% of the people from the fipronil-induced neurotoxicity of 10% extra risk. Likewise, we can easily tailor the pRfD with flexible M and I for different decision-making contexts, such as benefit-cost analysis, life-cycle impact analysis, and emergency responses to environmental incidents.
The ADI set by JMPR (0.0002 mg/kg/day) for fipronil is more conservative than the pRfD derived in this study, which is concordant with the results from the previous analysis on 1522 chemicals and endpoints [25]. From the BBMD analysis, approximately 12.6% of the uncertainty is associated with the POD, 57.9% for human variability, 28.4% for interspecies difference, and 1.1% for allometric scaling. The moderate uncertainty associated with POD indicates that the BMDL estimations from the eight models are similar, reporting BMDL estimates from one to 100 mg/kg/day. Overall, the BMDL values determined using BBMD are about 100-fold greater than the NOAEL of fipronil-induced neurotoxicity in male rats, resulting in a larger pRfD than the ADI.
In this study, we conclude that the health risk caused by fipronil exposure in Taiwan is minimal. Less than 0.4% of the population had LB or UB exposure estimates exceeding the traditional ADI. Additionally, only one out of the 200,000 simulated individuals reported a pHQ value greater than one, regardless of the LB or UB exposure scenario. In a closer look, this simulated individual is a 38 years-old male with a body weight of 74.1 kg. He was subject to a highly unusual exposure scenario, having a relatively high RV consumption (258 g/day) and an extremely high concentration of fipronil residue in RV. Nevertheless, the relatively high detection rates of fipronil in CG and RV necessitate continuous post-market investigation in the future.
This study is not without limitations. Firstly, the food consumption data were based on a 24-hour dietary recall questionnaire. Recall bias and misrepresentation of portion sizes may occur while collecting food consumption data. However, we perform the best-case (LB) and the worst-case (UB) exposure assessments to better describe the uncertainty in food consumption data. Secondly, the sample size of CG in fipronil analysis is relatively small (n = 40) compared with other fruits and vegetables. Continuous sampling and determination of fipronil concentrations in CG items are warranted. Lastly, we applied the default probabilistic factors to describe human variability in toxicokinetics and toxicodynamics (GM = 0.746, GSD = 1.59), which may underestimate or overestimate the "true" interindividual variability in fipronil's toxicokinetics and toxicodynamics. If accessible, CASF is highly recommended to replace the default probabilistic factor accounting for human variability. Notwithstanding these limitations, this case study demonstrated the utility of Bayesian inference in reducing the uncertainty in exposure assessment and dose-response analyses. We strongly advocate for adopting a probabilistic framework in the dietary risk assessments of foodborne chemicals.