Optimization of Sampling for Monitoring Chemicals in the Food Supply Chain Using a Risk-Based Approach: The Case of Aflatoxins and Dioxins in the Dutch Dairy Chain.

Food safety monitoring faces the challenge of tackling multiple chemicals along the various stages of the food supply chain. Our study developed a methodology for optimizing sampling for monitoring multiple chemicals along the dairy supply chain. We used a mixed integer nonlinear programming approach to maximize the performance of the sampling in terms of reducing the risk of the potential disability adjusted life years (DALYs) in the population. Decision variables are the number of samples collected and analyzed at each stage of the food chain (feed mills, dairy farms, milk trucks, and dairy processing plants) for each chemical, given a predefined budget. The model was applied to the case of monitoring for aflatoxin B1 /M1 (AFB1 /M1 ) and dioxins in a hypothetical Dutch dairy supply chain, and results were calculated for various contamination scenarios defined in terms of contamination fraction and concentrations. Considering various monitoring budgets for both chemicals, monitoring for AFB1 /M1 showed to be more effective than for dioxins in most of the considered scenarios, because AFB1 /M1 could result into more DALYs than dioxins when both chemicals are in same contamination fraction, and costs for analyzing one AFB1 /M1 sample are lower than for one dioxins sample. The results suggest that relatively more resources be spent on monitoring AFB1 /M1 when both chemicals' contamination fractions are low; when both contamination fractions are higher, relatively more budget should be addressed to monitoring dioxins.

INTRODUCTION

Exposure of consumers to safety hazards in food can cause a variety of health effects, and can affect agricultural industries and the national economy (Duret et al., 2019; WHO, 2015b). The WHO foodborne disease burden epidemiology reference group (FERG) reported that except bacteria, chemicals in food can also contribute to foodborne diseases (WHO, 2015a). Preventing these foodborne diseases requires prevention and control as well as monitoring of food safety chemicals along the food supply chain. Food industries have food safety standards, as well as prevention and control systems like Hazard Analysis Critical Control Points (HACCP) in place to manage food safety along the production chain (Trienekens & Zuurbier, 2008). Monitoring food safety hazards is part of HACCP, in particular the validation and verification steps (Surak, 2007). Risk managers in the food supply chain as well as governmental bodies make use of monitoring programs to verify if food safety is under control and to make decisions about food safety measures (Wu & Chen, 2018). Sampling for chemical hazards in the food supply chains are important parts of food safety monitoring programs (European Commission, 2006a, 2006d, 2017). The European commission (EC) has set legal maximum limits for the presence of chemicals in food and animal feed products (European Commission, 2002, 2006b, 2006d, 2011). In addition, prescribed procedures for sampling (sample collection, sample preparation, and sample analysis) and general risk assessment methodology are in place in Europe (European Commission, 2009, 2015, 2017).

Even though food safety sampling schemes have been designed in many countries across the world, the tradeoff between the effectiveness of food safety sampling and the available financial resources is insufficiently studied (Focker, Van der Fels‐Klerx, & Oude Lansink, 2018; Lascano‐Alcoser, Velthuis, Van der Fels‐Klerx, Hoogenboom, & Oude Lansink, 2013; Powell, 2014). Several studies have developed optimization models for designing cost‐effective food safety sampling but these studies focus on one chemical hazard only (Focker et al., 2018; Focker, Van der Fels‐Klerx, & Oude Lansink, 2019a; Lascano‐Alcoser et al., 2014). However, in practice, food safety authorities and food industries often need to monitor multiple chemicals along the food supply chain within the limited monitoring resources (Asselt, Van der Fels‐Klerx, Marvin, Bokhorst van de Veen, & Groot, 2017; Guo, Claassen, Oude Lansink, & Saatkamp, 2014). Therefore, more research should focus on optimal allocation of food safety resources in portfolio hazards monitoring.

The objective of our study was to develop a methodology for optimizing sampling for monitoring a portfolio of food safety chemicals along the supply chain using a risk‐based approach. As a case study, we focused on optimizing the number of samples for monitoring dioxins and aflatoxins along various stages of a hypothetical Dutch dairy supply chain, given their reduction in the risk of foodborne disease burden, and a predefined budget.

METHODOLOGY

In our study, we used a modeling approach to optimize sampling for two chemicals (i.e., aflatoxins and dioxins) along a hypothetical Dutch dairy chain using a risk‐based approach, and in this case risk was defined as estimated reduction in the risk of foodborne disease burden caused by aflatoxins and dioxins. Both aflatoxins and dioxins have similar contamination routes along the dairy supply chain, starting from the animal feeds and spreading to the downstream dairy chain stages. Also, humans are exposed to these two chemicals mainly through food consumption with significant health impacts (Asselt et al., 2017; Gibb et al., 2015). Definitions, legal limits and control measures of aflatoxins and dioxins are summarized below.

Dioxins refer to a group of toxic chemicals, including polychlorinated dibenzo‐p‐dioxins (PCDDs) and polychlorinated dibenzofurans (PCDFs), which are by‐products of industrial processes, or produced by natural phenomena. Dioxins can easily enter the dairy chain from the environment and feed, and spread throughout the chain through transportation, primary milk production, dairy factories, wholesalers, and retailers to consumers (Flores‐Miyamoto, Reij, & Velthuis, 2014; Hoogenboom et al., 2010). The most sensitive human health effect caused by exposure to dioxins are prenatal and postnatal hypothyroidy, and prenatal induced reduced sperm production. The disease burden related to dioxin exposure is measured as long time exposure for chronic toxicity risk, instead of daily exposure (WHO, 2015a). In the European Union, the legal limits of dioxins in animal feed and milk are set at 0.75 pg WHO98‐TEQ (toxic equivalents)/g feed1 (which was expressed as pg of TEQ/g feed in the following parts) (maximum level) and 2 pg WHO98‐TEQ/g fat (which was expressed as pg of TEQ/g milk fat in the following parts) (action level), respectively (European Commission, 2002, 2006b, 2011). The EC recommended analytical method for dioxins is a combination method, consisting of a screening method (e.g., CALUX assays [DR CALUX®]) followed by a confirmatory method (e.g., gas chromatography/high‐resolution mass spectrometry (GC/HRMS) (Lascano‐Alcoser et al., 2014; Lascano‐Alcoser et al., 2013).

Aflatoxins (types B1, B2, G1, and G2) are another group of relevant chemicals in dairy supply chain (Asselt et al., 2017). All four types of aflatoxins are equally toxic and carcinogenic (International Agency for Research on Cancer [IARC]), but, from these four aflatoxins, aflatoxin B1 (AFB1) is most prevalent in feed for dairy cows (Ismail, Riaz, Gong, Akhtar, & Sun, 2019; Trevisani et al., 2014). When AFB1 is present in feed, it can be converted by dairy cows resulting in the presence of AFM1 in the milk, which—upon consumption of dairy products—is harmful to human health (Van der Fels‐Klerx & Camenzuli, 2016). Exposure of humans to aflatoxins results from food consumption only (Gibb et al., 2015). Hepatocellular carcinoma (HCC) is the most significant clinical outcome of human exposure to aflatoxins. HCC is present most often in persons with chronic hepatitis B virus (HBV) infection and chronic aflatoxin exposure (Wild & Gong, 2010). The FERG developed a multiplicative model for estimating the effects of aflatoxin exposure and HBV infection on HCC development (Liu & Wu, 2010). In European regulations, the maximum limits for the presence of AFM1 in raw (unprocessed) milk is 0.05 μg/kg, and the maximum limit for AFB1 in all feed materials is 0.02 mg/kg, whereas it is 0.005 mg/kg for compound feed with a moisture content of 12% for dairy cows (European Commission, 2002, 2006c). An incremental analytical method, that is, liquid chromatography combined with mass spectrometry (LC‐MS/MS) with detection limit (0.0125 μg/kg), was recommended for the analysis of the aflatoxins in food products (Focker et al., 2019a).

Overview of the Methodology

In our modeling approach, disease burdens caused by contamination with each of two chemicals, were estimated using the models derived from FERG Disability‐Adjusted Life Years (DALYs) framework. This framework quantifies the burden of foodborne diseases in terms of DALYs, a health gap measure expressing the number of healthy life years lost due to reduction of health (quality of life) and death (Devleesschauwer et al., 2015). The sampling procedure for the two chemicals included the number of production units sampled (npu) (e.g., number of silos in feed mill, dairy farms, or milk trucks), the number of samples that need to be collected from each stage (ns), and analyzed (na) per chemical at each food supply chain stage, such that the contaminated samples (i.e., the concentration of the particular chemical is higher than the decision limits, which were set in Table IV) could be identified. An integrated supply chain was assumed implying that the different dairy chain stages share analytical results. It was assumed that, when one contamination was detected at any dairy chain stage, further tracking and tracing to the contamination sources would be done by corresponding stakeholders to identify and reject all related contaminated batches along the dairy chain; thus, the recall costs and destruction costs would not be included in the model. As a result, the identified contamination then did not result in any DALYs anymore, and these possible reduced DALYs were set as the performance of the monitoring sampling. The optimization model aimed to estimate the optimal number of npu, ns, and na at each supply chain stage such that the performance of the sampling for monitoring portfolio chemicals over the supply chain within predefined budgets was maximized.

Table IV

Model Variables and Input Values for Estimating the Costs and Performance of the Monitoring Sampling

Parameter	Value	Unit	Explanation
Nspuhi (nspu_feed)	7		Number of samples collected per feed mill. EC recommendation 1 , 2 & 3
Nspuhi (Nspu_milk) 10	3		Number of samples collected per milk production unit. EC recommendation 1 , 2 & 3
DL_h,i (DLdio,cf)	0.75	pg of TEQ/g in compound feeds	Decision limit of dioxins in compound feeds. The maximum level of dioxins concentration in compound feeds 4
DL_h,i (DLAFB1,cf)	0.005	mg/kg in compound feed	Decision limit of AFB1 in compound feeds. EU maximum level of AFB1 in compound feeds 5
DL_h,i (DLdio,m)	2	pg of TEQ/g in milk fat	Decision limit of dioxins in milk. The action level of dioxins concentration in milk 6
DL_h,i (DLAFM1,m)	0.05	μg/kg in milk	Decision limit of AFM1 in milk. EU maximum level of AFM1 in milk 7
PR_feed	1		Pooling rate of feed. We set one compound feed sample as one analysis sample 8
CostS	10	euros	Reference value 9
CostA_AFB1/M1	100	euros	LC‐MS/MS 9
Sen_A_h (Sen_ AFB1/M1)	100	%	Sensitivity of aflatoxins analysis. Assumed probability of LC‐MS/MS identifying positive samples
CostAS_dio	100	euros	DR CALUX® 8
CostAC_dio	350	euros	GC/HRMS 8
Sen_A_h (Sen_dio)	100	%	Sensitivity of dioxin analysis by combined method. Assumed probability of combined DR CALUX® with GC/HRMS identifying positive samples

European Commission (2006e).

European Commission (2006a).

European Commission (2009).

European Commission (2006b).

European Commission (2002).

European Commission (2011).

European Commission (2006d).

Lascano‐Alcoser et al. (2014).

Focker et al. (2019a).

Nspu_milk: number of samples collected in each milk production unit at different chain stages are all equal to 3 according to European Commission, 2006a, 2006e).

In the modeling approach (Fig. 1), first, the contamination fraction and concentration of each of AFB1/M1 and dioxins were estimated, starting from the initial contamination at the feed mills (FM) through the dairy farms (DF) up to the milk trucks (MT). It was assumed that the distribution of the chemicals in retail milk and in the milk in the truck tank did not differ (European Food Safety Authority, 2012; Trevisani et al., 2014; Van der Fels‐Klerx & Camenzuli, 2016). The exposure of chemicals in certain a population were calculated based on the contamination of retail milk with the chemical and milk consumption patterns, and then DALYs caused by each chemical were estimated separately (Assunção et al., 2018; Boon, te Biesebeek, de Wit‐Bos, & van Donkersgoed, 2014; Gibb et al., 2015; Liu & Wu, 2010). Samples were collected, pooled, and analyzed based on the EC recommended sampling scheme through the various dairy chain stages (FM, DF, and MT) (European Commission, 2006a, 2006e, 2009, 2017).

Fig 1

Illustration of the modeling procedure to obtain an optimal sampling for monitoring both aflatoxin B1/M1 and dioxins along the dairy supply chain.

Dairy Chain Structure

The modeling approach was applied to a hypothetical dairy supply chain in the Netherlands. The feed stage was composed of 10 FM, with 80 silos each. The FM receive new feed ingredients to produce compound feeds once per two weeks (Van der Fels‐Klerx & Camenzuli, 2016). The FM deliver compound feeds to 11,200 DF in the Netherlands; each dairy farm receives new compound feeds every two weeks, and each feed mill distributed feed to 80 (identical) DF per day during the monitoring period (Lascano Alcoser, Velthuis, Hoogenboom, & Van der Fels‐Klerx, 2011). At the dairy farm stage, the lactation period of the cows was set to 45 weeks (excluding the dry period), and thus, the monitoring period in one year was set to 45 weeks. Milk produced at one dairy farm was collected by trucks three times per week. The milk from four DF was mixed into one milk truck, and sent to one processing plant. The related model parameters are presented in Table I.

Table I

Model Parameters and Input Values for the Hypothetical Dairy Supply Chain in the Netherlands

Parameter	Value	Unit	Explanation
TN_f	10		Assumed input variable 1
S_fm	80		Assumed input variable 1
Time interval of FM get new ingredients	2	weeks	Common situation in Dutch dairy chain 2 , 3
Frc_DCF	85	%	Average value 3
Average number of DF getting compound feeds from one feed mill	80		Assumed input variable 1
TN_DF	11,200		Assumed input variable 1
TI_dfc	2	weeks	Common situation 2
The lactation period per cow	45	weeks	Average value 3
Dry period per cow	4	weeks	Average value 3
Cons_CF_day	4.3	kg	Average value 3
milk_fat	4	%	Nutrient fraction in raw milk 6
Qtruck	20,000	liters	Common situation 6
Qmilk	5,000	liters	Common situation 6
Frq_MTM	3	times/week	Common situation 6

Assumption based on database (FEFAC, 2016; ZuivelNL, 2016).

Lascano Alcoser et al. (2011).

Van der Fels‐Klerx and Camenzuli (2016).

Adamse, Van der Fels‐Klerx, Schoss, de Jong, and Hoogenboom (2015).

Malisch (2017).

Lascano‐Alcoser et al. (2013).

Contamination Scenarios and Estimation

Table II presents contamination scenarios for the concentration of AFB1 and dioxins in compound feeds. The initial contamination scenarios at FM were defined by nine combinations of AFB1/M1 contamination fractions and dioxins contamination fractions at FM (, being 1%, 5%, or 10%), and the concentrations of AFB1 and dioxins in compound feeds () were set at 20 ug/kg feed and 0.75 pg WHO98‐TEQ/g feed, respectively. Next, contaminations in the consecutive stages of the chain were estimated based on the assumed dairy chain structure and equations in Appendix A (A1a)–(A1j) with parameters shown in Tables I and II. Dairy cows were assumed to be fully housed indoors and only contaminated with dioxins and aflatoxins through intake of compound feeds, not by other sources. Based on animal feed consumption and transfer rates (carry‐over) for each of the two chemicals, the concentrations of AFM1 and dioxins in raw milk were estimated applying Equations (A1d) and (A1g) (Adekunte, Tiwari, & O'Donnell, 2010; Van der Fels‐Klerx & Camenzuli, 2016). It was assumed that one contaminated milk truck collected milk only from one contaminated farm (mixed with other noncontaminated farms), which is the worst‐case scenario for the contamination fraction. Another assumption was that the concentrations of AFM1 and dioxins in the consecutive chain stages (processing plants and retails) were similar to the concentrations in the bulk milk transported by MT (European Food Safety Authority, 2012; Trevisani et al., 2014; Van der Fels‐Klerx & Camenzuli, 2016). The concentrations of the two chemicals in noncontaminated milk were assumed to be equal to their respective background levels (which equals to the average concentration as determined in the national monitoring program or literature studies, Table II) (Lascano‐Alcoser et al., 2013; Van der Fels‐Klerx & Camenzuli, 2016). The distribution of AFM1 and dioxins in liquid milk was assumed to be homogeneously distributed (European Commission, 2006a, 2006e, 2017). All contamination scenarios were assumed to be stable in one year until the contamination was identified by the monitoring scheme.

Table II

Parameters and Input Values for the Estimation of the Contamination of AFB1/M1 and Dioxins Through the Dairy Supply Chain

Parameter	Value	Unit	Explanation
Fcf_ FMh	1; 5; 10	%	Assumed contamination level for each chemical
Cct_FM_AFB1	20	μg/kg in feed	Assumed input data
Cct_FM_dio	0.75	pg TEQ/g feed	Assumed input data
B_AFM1	0.037	μg/kg in milk	Average concentration of AFM1 in milk without any incidents 1
B_DM	0.5	pg of TEQ/g in milk fat	Concentration of dioxins in noncontaminated milk 2
CO_dio	40	%	Transfer rate of dioxins from animal feed to bovine milk 3 , 4

Van der Fels‐Klerx and Camenzuli (2016).

Lascano‐Alcoser et al. (2013).

Malisch (2017).

Adekunte et al. (2010).

Exposure Estimation

The exposure module calculated the estimated intake of each of the two chemicals via consumption of milk produced by the hypothetical dairy supply chain by 100,000 Dutch people on a yearly basis. The exposure estimation started with the calculation of the fraction of contaminated products (with a concentration at legal limits or higher) of each of the two chemicals in the food products, as compared to the reference situation (concentrations at background level). The exposure estimation2 was calculated in Equations (1) and (2) with parameters explained in Tables II and III and Appendix A (Adekunte et al., 2010; Liu & Wu, 2010). where Extra_EXP is the human exposure to AFM1 by consumption of contaminated milk (ng/kg bw/day); Cct_R_AFM is the concentration of AFM1 in the retail milk (μg/kg) (Appendix A); B_AFM is the background level of AFM1 in milk (μg/kg); Cons_Milk is the average milk consumption of Dutch people per day (g/day); Milk_fat is the fraction of fat in milk (4%); Prob_CnAM is the fraction of contaminated milk with a AFM1 concentration higher than decision limit, which equaled the Fcf_R (contamination fraction in retail milk for chemical h, shown in Appendix A); BW is the average body weight of a Dutch person (kg); Extra_EXP is the human exposure to dioxins by consumption of contaminated milk (pg of TEQ/kg bw/day); Cct_R_Dio is the concentration of dioxins in the contaminated retail milk (pg of TEQ/g in milk fat) (Appendix A); B_DM is the background level of dioxins in milk (pg of TEQ/g in milk fat); Prob_CnDM is the fraction of contaminated milk with dioxins higher than decision limit, which equaled the Fcf_R (contamination fraction in retail milk for chemical h, shown in Appendix A).

Table III

Input Values for Exposure Estimation and DALYs Calculation of Chemicals

Parameter	Value	Unit	Explanation
BW_m; BW_f	Male:85; Female:74	kg	Body weight of Dutch people. The average value calculated from the source 1
Cons_milk	258	g/day	Dutch milk consumption. 7–68 years old 1
Tf	0.5	years	The average value 2
Dwf	0.508		Terminal phase with medication; reference from FERG 3
Pf	92	%	The average value 2
t_nf	5	years	The average value 2
Dwnf	0.294		Diagnosis and primary therapy; reference from FERG
LE_m; LE_f	75.2; 80.5	Years	The life expectancy for Dutch male and female. Average value 4
Age_HCC	60	years	Refence age 5
Sex_dist_f/m	50:50	%	An equal age distribution was assumed 6
Rfa‐	0.01	cases	Cases/100,000/year/ng/kg bw/day aflatoxin exposure for individuals without HBV infection 7
Rfa+	0.3	cases	Corresponding cases for individuals with HBV infection 7
Phbv	0.5	%	Worst case were assumed 8
RR_AFM1_AFB1	10	%	Reference value 9
EXP_DF	1.3	pg of TEQ/kg bw per day	Reference value 10
IRi (IR_Hpre, IR_Hpos, IR_Infer)	0.010968, 1.29498, 0.0109686		Diseases incidences rate caused by exposure estimation of dioxins for: infertility; hypothyroidy due to prenatal exposure; hypothyroidy due to postnatal exposure. Reference value 3
DWi (DW_Hpre, DW_Hpos, DW_Infer)	0.019, 0.019, 03056		Disability weight of diseases caused by dioxins for: infertility; hypothyroidy due to prenatal exposure; hypothyroidy due to postnatal exposure. Reference value 3
DRDi (DRD_Hpre, DRD_Hpos, DRD_Infer)	80, 60, 25	Years	Duration until remission or death for: Infertility; hypothyroidy due to prenatal exposure; hypothyroidy due to postnatal exposure. Reference value 3

Van Rossum et al. (2011).

Assunção et al. (2018).

Devleesschauwer et al. (2015).

Van Kreijl et al. (2006).

Verhoef et al. (2004).

WHO (2015b).

Liu & Wu (2010).

Koopsen et al. (2019).

Wu, Narrod, Tiongco, and Liu (2011).

Boon et al. (2014).

DALYs Calculation

The disease burden, expressed in DALYs, related to human exposure to each of the two chemicals via dairy milk consumption was calculated using Equations (3)–(10) with the estimated contamination and exposure as inputs. To this end, the DALYs calculation module was derived from the FERG's methodological framework together with Assunção et al. (2018); Boon et al. (2014); Gibb et al. (2015); Liu and Wu (2010). Basic DALYs calculations (used in the FERG framework) can be found in the Appendix A (Equations (A2a)–(A2d)).

Aflatoxins

For aflatoxins, HCC is the only considered clinical outcome that is related to aflatoxin exposure in this case (Devleesschauwer et al., 2015). To estimate aflatoxin‐induced HCC cases per 100,000 in the two sets of populations (those with and without chronic HBV infection), the exposure estimates were classified by the corresponding cancer potency factor in Equations (6) and (7). Next, these values were multiplied with the nations’ HBV+/HBV− population proportions to arrive at the total number of aflatoxin‐induced HCC cases in the 100,000 Dutch population. Based on a review of multiple studies examining cancer potency in human populations, JECFA (1998) selected two different cancer potency factors for aflatoxins, being 0.01 cases per 100,000 per year for every ng/kg bw/day aflatoxin exposure for individuals without chronic HBV infection, and 0.30 corresponding cases for individuals with chronic HBV infection (Liu & Wu, 2010). The average age at onset of this disease was assumed to be 49 years. In absence of information on the sex distribution of aflatoxin‐induced HCC, a 50:50 male–female distribution was assumed. The DALYs calculation of AFM1 in milk for 100,000 diseased cases in the Dutch people is presented by Equations (3)–(7), the corresponding parameters and values are shown in Table III.

where DALYs_AM are DALYs caused by the AFM1 in contaminated milk (DALYs/100,000 population); HCC are total cases of HCC per 100,000 exposure population; HCC are total cases of HCC with positive HBV infection per 100,000 exposure population; HCC are total cases of HCC with negative HBV infection per 100,000 exposure population; tf is the duration of disease of fatal cancer (years); dwf is the disability weight (0–1) of fatal cancer, in accordance with Devleesschauwer et al. (2015); pf is the probability of a cancer being fatal (%); tnf is the duration of disease of nonfatal cancer (years); dwnf is the disability weight nonfatal cancer, in accordance with Devleesschauwer et al. (2015); YLL is the life years lost due to premature death to a fatal cancer (years); Age_HCC is the reference onset age of HCC; LE_m is the life expectancy of Dutch male (years); LE_f is life expectancy of Dutch female (years); Sex_dist_f/m is in absence of information on the sex distribution of aflatoxin‐induced HCC, a 50 is50 gender distribution was assumed; the percentage of female/male for this disease (%); Phbv is the proportion of Dutch population which is HBV positive (%); Rf is a risk factor of negative HBV infection in HCC, which means cases per 100,000 per year for every ng/kg bw/day AFB1 exposure for individuals without chronic HBV infection; Rf is a risk factor of positive HBV infection in HCC, which means cases for individuals with chronic HBV infection; RR is a relative risk factor regarding AFB1 to AFM1, which can transfer the toxic potency of AFM1 to the toxic potency of AFB1; Extra_EXP is the human exposure to AFM1 by consumption of contaminated milk (ng/kg bw/day).

Dioxins

Hypothyroidism due to prenatal or postnatal exposure and male infertility were considered the clinical outcomes of human exposure to dioxins through contaminated food, in this case milk. Hypothyroidy was assumed to be lifelong after its onset. The effect on male infertility was assumed to be present in the 20–44 age group (in accordance with Van Rossum, Fransen, Verkaik‐Kloosterman, Buurma‐Rethans, & Ocke, 2011). In absence of information on the sex distribution of dioxin‐induced hypothyroidy, a 50:50 gender distribution was assumed. For male infertility, the entire burden was assigned to males. The disability weight (DW) values for these three clinical outcomes were taken from Devleesschauwer et al. (2015). Due to the lipophilic characteristics of dioxins, daily dietary exposure leads to accumulation of these chemicals in human body fat. Consequently, the dioxin body burden, rather than the daily exposure, is taken as yearly basis exposure for a certain population. Boon et al. (2014) estimated DALYs caused by dioxin exposure via food consumption, considering different contaminated food products (Boon et al., 2014). With that information, the contribution that these DALYs were due to consumption of milk was calculated proportional to the estimated exposure, using Equations (8)–(10). The related parameters are shown in Table III.

where Extra_ EXP is the extra exposure due to the dioxins in contaminated milk (pg of TEQ/kg bw/day); EXP_DF is the exposure estimation for dioxins in all food products (pg of TEQ/kg bw/day); Cct_R_dio is the dioxins concentration in contaminated milk at retail (pg of TEQ/g in milk fat); DALYs_DF are DALYs due to dioxin exposure in all foods in a population (DALYs/100,000 population); DALYs_DM are DALYs caused by the dioxins in contaminated milk (DALYs/100,000 population); is the disease incidence rate caused by exposure estimation with i for types of disease caused by exposure to dioxins (infertility, hypothyroidy due to prenatal exposure, and hypothyroidy due to postnatal exposure); is the disability weight with i for diseases (infertility, hypothyroidy due to prenatal exposure and hypothyroidy due to postnatal exposure); is the duration until remission or death (years), due to disease i (infertility, hypothyroidy due to prenatal exposure and hypothyroidy due to postnatal exposure).

Sampling Procedure

The sampling procedure for monitoring the two chemicals in the dairy supply chain was considered to consist of the three steps of sample collecting, pooling, and analyzing. Incremental samples were collected randomly from FM, DF, and MT of the dairy supply chain. The number of incremental samples (nspu) collected from each production unit (e.g., one feed mill, one dairy farm, or one bulk milk truck) at stage i were set according to the EC requirements (Table IV). Collected incremental samples were transported to the laboratory where a certain number of collected samples from different production units at same dairy chain stage were pooled into one aggregated sample applying a certain pooling rate. The incremental milk samples from different production units at same stage can be mixed together, but the incremental feed samples can only be mixed within the same production unit (one incremental compound feed sample were regarded as one aggregate sample). From the aggregate sample, the aliquot or analytical sample (na) for chemical h, at stage i was obtained and analyzed, and in this case, we assumed that the number of analytical samples equaled the number of aggregate samples (so, one analytical sample derived from one aggregate sample). To identify dioxins, the analytical sample was initially analyzed by the DR CALUX® screening method. In case the screening result of analytical sample was positive, that is, the concentration of dioxins in the analyzed sample was higher than a preset decision limit (Table IV), the suspected positive sample was analyzed using the confirmatory method, being GC/HRMS, and individual samples from positive aggregate sample were all analyzed. In this combined method, the screening method was applied to all laboratory samples, and the use of the (more expensive) confirmatory method depended on the screening results of dioxins analytical samples, and detection limits of these methods were assumed to be smaller than the background level of dioxins in products. To identify AFB1 and AFM1, all aggregated samples were analyzed by LC‐MS/MS and its detection limits were assumed to be smaller than background level of AFB1/M1 in products. The requirements and parameters for the sampling procedure are described in Table IV. It was assumed that the batches at any stage with one contamination would be rejected.

Optimization

The optimization model maximizes the possible reduction of the expected DALYs through a sampling procedure that monitors two chemicals in the dairy supply chain during the monitoring period of one year. The possible reduction of expected DALYs was calculated by multiplying original expected DALYs caused by contaminated milk (estimated in the DALYs calculation) and the probability of monitoring sampling that identify a contamination through the entire dairy supply chain (Equation (11)). The optimization module was developed using mixed integer nonlinear programming in Microsoft Excel 2016 and the Solver command from Frontline System Inc. (2015). Related equations are described below, see Equations (11)–(27), with parameters and values shown in Table IV. The sampling procedure were optimized within a preset budget of 10,000 euros (estimated from the results of Focker, Van der Fels‐Klerx, & Oude Lansink, 2019b and Lascano‐Alcoser et al., 2013) and followed by a cost‐effectiveness analysis. The effectiveness is defined as fraction of DALYs saved by optimal monitoring for each contamination scenario under different monitoring budgets (2,000–12,000 €).

Max is

subject to:

with:

where TSP is the total performance of sampling for monitoring both aflatoxin and dioxins in the dairy supply chain (total reduced DALYs/100,000 population); H is the number of chemicals involved in optimization; I is the stage of the dairy supply chain, that is, FM, DF, and MT; DALY are the DALYs caused by the consumption of milk in which the chemical concentration is higher than the decision limits (e.g., legal limit); is the probability of monitoring sampling that identify a contamination with chemical h along the entire dairy chain; prob ( npu ) is the probability of monitoring sampling that can identify at least one chemical h contamination at stage i, following hypergeometric distribution; the probability of collecting noncontaminated samples with chemical h, at stage i, following hypergeometric distribution; is the sensitivity of analyzing each chemical in food or feed products; Fcf is the contamination fraction at each dairy chain stage for each chemical; is the total production units at each dairy chain stage that produce or collect new products in monitoring period, that is, number of silo in FM (), number of farms in DF (), and number of bulks in MT (; is the total contaminated production units at each dairy chain stage for each chemical in monitoring period; TN_f is total number of FM producing feed for dairy cattle, 10 FM; S_fm is number of silos in each feed mill; TN_DF is total number of DF; TI_dfc is time interval of DF getting new compound feeds (weeks); Qtruck is milk truck capacity; Qmilk is milk collected per each farm per delivery; Frq_MTM is frequency of MT collecting milk (times/week); is the number of production units sampled at stage i for chemical h in the monitoring period (e.g., number of farms in dairy farm stage); is the number of samples collected at stage i for chemical h in the monitoring period; Max_ns is the maximum number of samples collected at stage i for chemical h in the monitoring period, which equals to ; is the number of samples collected from one production unit at stage i for chemical h; is the maximum number of samples collected from one production unit at stage i for chemical h, which was set in Table IV; is the number of samples analyzed at laboratory from stage i for chemical h in the monitoring period; is the maximum number of samples analyzed at laboratory from stage i for chemical h in the monitoring period, which was determined by both and ; are the costs of sampling for chemical h, at stage i, that is, are the costs of sampling AFB1 at stage i and (euros) are the costs of sampling dioxins at stage i; CostS are the costs of collecting one sample (euros); CostA_AFB are the costs of analyzing one AFB1/M1 sample (euros); CostAS_dio are thecosts of analyzing one dioxins sample by screening method (euros); CostAC_dio are the costs of analyzing one dioxins sample by confirmatory method (euros); Budgets are the budgets set by food safety authority to monitoring chemicals (euros); is the pooling rate, the number of collected samples mixed into one analysis sample; is the concentration of chemicals in the contaminated products at stage i; is the background level of chemical h in milk at stage i, for example, B_DM and B_AFM; is the decision limit of chemical h in products at stage i; Int(x) is the function which rounds off the results as integer.

RESULTS

Table V presents, for each of the nine initial contamination scenarios, the spread of the contamination with dioxins and aflatoxins through FM, DF, and MT of the hypothetical Dutch dairy supply chain, the exposure of the Dutch population to these two chemicals, and the related expected DALYs per 100,000 population. With identical initial contamination fraction of the two chemicals, the expected DALYs caused by AFB1/M1 were higher than those for dioxins. Scenario S1 was the scenario with the lowest contamination fractions for each of AFB1/M1 and dioxins. This scenario also resulted into the lowest expected DALYs, being 0.046/100,000 population. In this scenario, expected DALYs caused by AFB1/M1 were about two times higher than those caused by dioxins. Scenario S3 was the worst‐case scenario considered in our study, with the highest initial contamination fractions for both chemicals. S7 was the most common scenario with a low contamination fraction for dioxins and a high contamination fraction for AFB1, which best reflects the current contamination situation. Total expected DALYs of this scenario were 0.328/100,000 population (which was around eight times higher than in S1). The concentration of chemicals in milk became lower in the MT due to mixing of milk from various trucks and related dilution of the contamination.

Table V

Nine Scenarios of Dioxins and aflatoxin B1/M1 (AFB1/M1) with Different Contamination Fraction and Preset Concentrations Along the Dairy Supply Chain and Their Corresponding Exposure Estimation and DALYs for 100,000 Population in the Netherlands

Contamination	FM		DF		MT		Extra Exposure estimation	Expected DALYs/(100,000 population)
Scenarios	CF	AFB1 (μg/kg in feed); Dioxins (pg of TEQ/g in feed)	CF	AFM1 (μg/kg in milk); Dioxins (pg of TEQ/g in milk fat)	CF	AFM1 (μg/kg in milk); Dioxins (pg of TEQ/g in milk fat)	ng/kg bw/day; pg of TEQ/kg bw/day	AFM1‐DALYs; Dioxins‐DALYs; Total DALYs
S1	0.01	20	0.01	0.09	0.04	0.05	1.720	0.0308
	0.01	0.75	0.01	7.5	0.04	2.25	0.009	0.0152
								0.0460
S2	0.05	20	0.05	0.09	0.2	0.05	8.600	0.1567
	0.05	0.75	0.05	7.5	0.2	2.25	0.045	0.0527
								0.2095
S3	0.1	20	0.1	0.09	0.4	0.05	17.200	0.3130
	0.1	0.75	0.1	7.5	0.4	2.25	0.091	0.1050
								0.4180
S4	0.05	20	0.05	0.09	0.2	0.05	8.600	0.1567
	0.01	0.75	0.01	7.5	0.04	2.25	0.009	0.0152
								0.1720
S5	0.1	20	0.1	0.09	0.4	0.05	17.200	0.3130
	0.05	0.75	0.05	7.5	0.2	2.25	0.045	0.0527
								0.3657
S6	0.01	20	0.01	0.09	0.04	0.05	1.720	0.0308
	0.1	0.75	0.1	7.5	0.4	2.25	0.091	0.1050
								0.1358
S7	0.1	20	0.1	0.09	0.4	0.05	17.200	0.3130
	0.01	0.75	0.01	7.5	0.04	2.25	0.009	0.0152
								0.3282
S8	0.01	20	0.01	0.09	0.04	0.05	1.720	0.0308
	0.05	0.75	0.05	7.5	0.2	2.25	0.045	0.0527
								0.0835
S9	0.05	20	0.05	0.09	0.2	0.05	8.600	0.1567
	0.1	0.75	0.1	7.5	0.4	2.25	0.091	0.1050
								0.2617

FM = feed mills.

DF = dairy farms.

MT = milk trucks.

CF = contamination fraction.

Tables VI and VII present optimal outcomes of decision variables given a preset budget of 10,000 euros from different aspects. Table VI shows the optimal monitoring settings, including production units to be sampled (npu), number of samples (ns), number of analysis samples (na) for each contamination scenario, and corresponding expected reduction in DALYs per 100,000 population according to monitoring. Table VII shows the optimal allocation of a budget of 10,000 euros, the probability of identifying the presence of each chemical at each dairy supply chain stage, and the expected remaining DALYs after implementing the optimal sampling procedure. From Table VII, it can be seen that—for each contamination scenario considered—it was optimal to allocate most of the budget to sampling at the MT stage, which are also reflected in Table VI with most collected and analyzed samples at MT. The lower fraction of contamination scenario, the more DALYs would be left after implementing the optimal sampling procedure. In the contamination scenarios S1, S2, and S3, in which the contamination fraction of both AFB1 and dioxins was higher, the expected remaining DALYs were lower (Table VII). Both Table VI and VII show that the expected remaining DALYs for identifying the worst case scenario of S3 was the lowest among all scenarios, and all total DALYs were almost removed by implementing the corresponding optimal sampling procedures. When the contamination fraction of both chemicals was higher (e.g., S1, S2, and S3), the budget allocation moved from monitoring AFB1/M1 to dioxins (Table VII). Correspondingly, Table VI shows that the optimal sampling monitoring collected and analyzed the most samples for AFB1/M1 to reduce DALYs in S1; however, the optimal monitoring for S2 and S3 collected and analyzed more samples for dioxins. In the most common scenario of S7, the optimal results showed that most of the budget (7,415 euros) was spent on monitoring dioxins with the chain probability (62%) of identifying dioxins. When the dioxins fraction was at the low level and AFB1/M1 fraction was higher in series of S1,S4, and S7, the budgets optimally spent on monitoring AFB1/M1 and remaining DALYs caused by AFB1/M1 was lower (Table VII). This is consistent with results of Table VI, in which the optimal number of samples lower for AFB1/M1 and higher for dioxins with the same budgets.

Table VI

Optimal Number of Production Units and Samples Collected at Each Control Point in the Dairy Supply Chain for Aflatoxin B1/M1(AFB1/M1) and Dioxins Sampling Settings, As Well As DALYs/100,000 Population Saved by These Sampling, Given a Preset Budget of 10,000 Euros

Contamination		FM			DF			MT			Total DALYs reduced/
Scenarios		Npuh,i 1	Nsh,i 2	NAh,i 3	Npuh,i 1	Nsh,i 2	Nah,i 3	Npuh,i 1	Nsh,i 2	Nah,i 3	100,000 population
S1	AFB1/M1	1	7	1	4	12	1	40	120	40	0.031
	Dioxins	1	7	1	12	36	3	8	24	8
S2	AFB1/M1	1	7	1	1	3	1	23	69	23	0.21
	Dioxins	1	7	1	12	36	3	10	30	10
S3	AFB1/M1	1	7	1	24	72	6	13	39	13	0.42
	Dioxins	1	7	1	4	12	1	12	36	12
S4	AFB1/M1	1	7	1	7	21	2	23	69	23	0.16
	Dioxins	1	7	1	24	72	6	13	39	13
S5	AFB1/M1	1	7	1	4	12	1	13	39	13	0.36
	Dioxins	1	7	1	8	24	2	14	42	14
S6	AFB1/M1	1	7	1	20	60	5	29	87	29	0.13
	Dioxins	1	7	1	3	9	1	7	21	7
S7	AFB1/M1	1	7	1	3	9	1	13	39	13	0.32
	Dioxins	1	7	1	35	105	9	15	45	15
S8	AFB1/M1	1	7	1	2	6	1	29	87	29	0.07
	Dioxins	1	7	1	12	36	3	8	24	8
S9	AFB1/M1	1	7	1	3	9	1	26	78	26	0.26
	Dioxins	1	7	1	6	18	2	9	27	9

FM = feed mills.

DF = dairy farms.

MT = milk trucks.

Optimal number of production units sampled at each control point for each chemical.

Optimal number of samples collected at each control point for each chemical.

Optimal number of analysis samples at each stage for each chemical.

Table VII

Optimal Expected DALYs/100,000 Population After Implementing the Optimal Sampling Procedure of Monitoring Dioxins and Aflatoxin B1/M1 (AFB1/M1) With the Allocation of Budgets 10,000 Euros at Different Stages of the Dairy Supply Chain

Contamination		Costs of sampling				Probability of identifying chemicals				DALYs left 2
Scenarios		FM	DF	MT	Chain 1	FM (%)	DF (%)	MT (%)	Chain 1 (%)	/100,000 population
S1	AFB1/M1 low	170	220	5780	7190	1	4	80	81	0.0057
	Dioxins low	174	780	1821	2775	1	11	28	37	0.0096
	TOTAL	344	1030	7601	9965					0.0154
S2	AFB1/M1 mid	170	130	3910	4210	5	5	99	99	0.0008
	Dioxins mid	188	1143	4426	5757	5	46	89	95	0.0029
	TOTAL	358	1273	8336	9967					0.0037
S3	AFB1/M1 high	170	1320	2210	3700	10	92	100	100	0
	Dioxins high	205	340	5752	6297	10	34	100	100	0.0001
	TOTAL	375	1660	7962	9997					0.0001
S4	AFB1/M1 mid	170	410	3565	4490	5	30	99	100	0.0006
	Dioxins low	174	1769	5752	5508	1	21	41	54	0.007
	TOTAL	344	2179	9317	9998					0.0076
S5	AFB1/M1 high	170	220	2210	2600	10	34	100	100	0.0002
	Dioxins mid	188	676	6504	7368	5	34	96	97	0.0015
	TOTAL	358	896	8714	9968					0.0017
S6	AFB1/M1 low	170	1100	4930	6200	1	18	69	75	0.0076
	Dioxins high	205	285	3291	3781	10	27	97	98	0.0019
	TOTAL	375	1385	8221	9981					0.0096
S7	AFB1/M1 high	170	190	2210	2570	10	27	100	100	0.0003
	Dioxins low	174	2886	4355	7415	1	30	46	62	0.0057
	TOTAL	344	3076	6565	9985					0.006
S8	AFB1/M1 low	170	160	4930	5260	6	7	71	75	0.0078
	Dioxins mid	194	1154	3381	4729	7	47	84	92	0.0084
	TOTAL	364	1314	8311	9989					0.0163
S9	AFB1/M1 mid	170	190	4420	4780	5	14	100	100	0.0004
	Dioxins high	205	708	4289	5202	10	47	99	100	0.0005
	TOTAL	375	898	8709	9982					0.0009

FM = feed mills.

DF = dairy farms.

MT = milk trucks.

Chain: from the integrated chain level, different stages jointly to identify the contamination.

DALYs left: the expected left DALYs refer to equation (A2d).

Figure 2 presents the effectiveness of sampling in reducing DALYs caused by AFM1 and dioxins in the dairy supply chain, with different total monitoring budgets. Here, the effectiveness was defined as the fraction of reduced DALYs over the original expected DALYs. The results show that—for all contamination scenarios—more DALYs would be reduced with a higher monitoring budget. When the contamination fraction of both AFB1/M1 and dioxins decreased (S3, S2, and S1), the increase in the effectiveness with higher total monitoring budgets became more pronounced. In S3, the total effectiveness reached 100% and remained at this level, after total budgets exceeded 4,000 euros. In S1, the effectiveness increased slowly with higher total budgets, and even at the highest total budgets of 12,000 euros, the total effectiveness only reached at around 75%. In S1, S4, and S7 (low fraction of dioxins with three different fractions of AFB1/M1), the effectiveness of monitoring dioxins contributed only to small extent to the total effectiveness for all three scenarios. Even though the total effectiveness reached a very high level at low budgets in S4 and S7, the effectiveness was still below 100% with the highest considered total budget. For monitoring dioxins and AFB1/M1 at different contamination scenarios, except for S6 and S8, the effectiveness for monitoring AFB1/M1 accounted for the major parts of the total effectiveness. In S6 (with low fraction of AFB1/M1 and high fraction dioxins) and S8 (with low fraction of AFB1/M1 and middle fraction of dioxins), the effectiveness for monitoring dioxins contributed for the major part of the total effectiveness, but the difference between effectiveness of monitoring AFB1/M1 and monitoring dioxins became slightly smaller when the total budgets increased from 6,000 euros to 12,000 euros.

Fig 2

The effect of total monitoring budgets on the fraction of DALYs saved (effectiveness) by optimal sampling for monitoring aflatoxin B1/M1 and dioxins for different contamination scenarios.

Results of the sensitivity analysis are shown in the Appendix (Tables B1 and B2) for the adapted value of the input parameters on the decision limits and sensitivity of detection methods, respectively. To assess the impact of the value for the decision limits on the optimal sampling settings and their performance for identifying different contamination scenarios, the decision limits for AFM1 and dioxins decreased to 0.04 μg/kg in milk and 1.6 pg of TEQ/g in milk fat, respectively. The results for all contamination scenarios show that more DALYs are saved with lower decision limits. Especially for S6 and S8, the saved DALYs increased from 0.126 and 0.071 DALYs/100,000 population (Table VI) to 0.1316 and 0.0761 DALYs/100,000 population (Table B1). The ns at DF and MT became much higher for identifying the same contamination scenario compared with using the input value presented in Table VI. With the same npu or ns at DF and MT, the na in Table B1 became smaller than the value in Table VI. As shown in Table B1, the highest na was still at MT for both chemical chemicals, and the ns and na were lowest for both AFM1 and dioxins at FM. Most of the samples were collected for identifying dioxins at DF, and for identifying AFM1 at MT.

Table B1

The Sensitivity Analysis Results Within 10,000 Euros Monitoring Budgets When the Decision Limits of AFM1 and Dioxins Decreased to 0.04 μg/kg in Milk and 1.6 pg of TEQ/g in Milk Fat

Contamination		FM 1			DF 2			MT 3			Total DALYs reduced/
Scenarios		Npuh,i 4	Nsh,i 5	NAh, i 6	Npuh,i 4	Nsh,i 5	Nah,i 6	Npuh,i 4	Nsh,i 5	Nah,i 6	100,000 population
S1	AFB1/M1	1	7	1	6	18	1	60	180	15	0.04
	Dioxins	1	7	1	30	90	5	8	24	8
S2	AFB1/M1	1	7	1	3	9	1	25	75	7	0.21
	Dioxins	1	7	1	17	51	3	12	36	12
S3	AFB1/M1	4	28	4	3	9	1	13	39	13	0.42
	Dioxins	4	28	4	8	24	2	11	33	11
S4	AFB1/M1	1	7	1	7	21	1	19	57	5	0.17
	Dioxins	1	7	1	54	162	9	11	33	11
S5	AFB1/M1	1	7	1	1	3	1	28	84	7	0.36
	Dioxins	1	7	1	23	69	4	13	39	13
S6	AFB1/M1	1	7	1	9	27	1	59	177	15	0.13
	Dioxins	1	7	1	24	72	4	3	9	3
S7	AFB1/M1	1	7	1	2	6	1	16	48	4	0.32
	Dioxins	1	7	1	60	180	10	13	39	13
S8	AFB1/M1	1	7	1	10	30	1	50	150	13	0.08
	Dioxins	1	7	1	24	72	4	6	18	6
S9	AFB1/M1	1	7	1	15	45	1	32	96	8	0.26
	Dioxins	1	7	1	18	54	3	9	27	9

Optimal number of production units sampled at each control point for each chemical.

Optimal number of samples collected at each control point for each chemical.

Optimal number of analysis samples at each stage for each chemical.

Table B2

The Sensitivity Analysis Results Within 10,000 Euros Monitoring Budgets When the Sensitivity of Agnatical Method Decreased to 98% and 98% for AFB1/M1 and Dioxins Separately

Contamination		FM 1			DF 2			MT 3			Total DALYs reduced/
scenarios		Npuh,i 4	Nsh,i 5	NAh,i 6	Npuh,i 4	Nsh,i 5	Nah,i 6	Npuh,i 4	Nsh,i 5	Nah,i 6	100,000 population
S1	AFB1/M1	1	7	1	11	33	3	38	114	38	0.03
	Dioxins	1	7	1	12	36	3	8	24	8
S2	AFB1/M1	1	7	1	39	117	10	15	45	15	0.20
	Dioxins	1	7	1	16	48	4	8	24	8
S3	AFB1/M1	4	28	4	36	108	9	8	24	8	0.42
	Dioxins	1	7	1	20	60	5	7	21	7
S4	AFB1/M1	1	7	1	56	168	14	14	42	14	0.16
	Dioxins	1	7	1	20	60	5	11	33	11
S5	AFB1/M1	1	7	1	32	96	8	8	24	8	0.36
	Dioxins	1	7	1	16	48	4	11	33	11
S6	AFB1/M1	1	7	1	12	36	3	32	96	32	0.12
	Dioxins	1	7	1	8	24	2	6	18	6
S7	AFB1/M1	3	21	1	32	96	8	9	27	9	0.32
	Dioxins	1	7	1	19	57	5	16	48	16
S8	AFB1/M1	1	7	1	1	3	1	30	90	30	0.07
	Dioxins	1	7	1	12	36	3	8	24	6
S9	AFB1/M1	1	7	1	48	144	12	14	42	14	0.26
	Dioxins	1	7	1	12	36	3	7	21	7

Expert opinion (personal communication).

Lascano‐Alcoser et al. (2013)

Optimal number of production units sampled at each control point for each chemical.

Optimal number of samples collected at each control point for each chemical.

Optimal number of analysis samples at each stage for each chemical.

Table B1 shows the effects of lower sensitivity of detection methods (98%) on optimal sampling settings and their performance. The DALYs reduced by implementing the optimal sampling procedures generally decreased, and more samples needed to be collected (ns) and analyzed (na) as compared to the original detection methods settings. Table B2 shows that more samples were collected and analyzed for identifying AFM1 than dioxins at both DF in most scenarios except scenario S8. No major changes were seen in the optimal number of collected and analytical samples for monitoring dioxins at each dairy chain stage.

DISCUSSION

In our study, we developed an optimization model to determine optimal sampling settings for monitoring different chemicals along the food supply chain based on a risk‐based perspective. We considered the disease burden in the population associated with the chemicals as the risk related to presence of the chemical in the dairy chain. Several previous studies developed methods for optimal allocation of resources in food safety monitoring programs (Focker et al., 2018). However, this is the first study that quantifies the monitoring for multiple food safety chemicals along the food supply chain and the first study that used the DALYs reduction as approach for risk‐based monitoring. In our study, the focus was on monitoring different chemicals using an integrated chain approach, that is, multiple chain stages share their monitoring responsibility. We considered such a chain approach since it has been shown to be more cost‐effective than monitoring at individual food chain stages (Lascano‐Alcoser et al., 2014). FERG developed a methodological framework to estimate global burdens of foodborne diseases, encouraging countries to use DALYs estimates for cost‐effective analysis of control measures (Devleesschauwer et al., 2015). Based on the framework of FERG, our model estimated the DALYs caused by chemical contaminations in the dairy chain, and these DALYs were used as inputs in the optimization procedure to optimize sampling for the considered chemicals. To our best knowledge, it is the first study that optimized the portfolio chemical hazards monitoring along the food supply chain from both epidemiological and economic aspects.

We applied the model to the sampling for monitoring AFB1/M1 and dioxins in a hypothetical dairy supply chain in the Netherlands. The results of our case study showed that the expected DALYs caused by the two chemicals increased with the contamination fraction. When the level of both chemicals in milk were similar, the expected DALYs caused by AFM1 were higher than those for dioxins (Table V). For all contamination scenarios, the higher the DALYs caused by the chemical in the baseline situation, the more DALYs would be reduced by implementing their corresponding optimal sampling procedure. Even though most of the optimal number of collected and analyzed samples were distributed differently in DF and MT stages for all scenarios (Table VI), the model outcomes suggest that largest part of the budget are allocated to MT stages in the optimal sampling settings (given the 10,000 euro budget) (Table VII). Since compound feeds were regarded as one mixed sample and the number of samples collected from one production unit (silo) in FM was higher than that from other considered chain stages, the costs of collecting and analyzing one production unit at later stages would be lower than that at FM. This is also the reason why the smallest number of samples and budgets were used at FM (Tables VI and VII). This conclusion is consistent with Trevisani et al. (2014) who found that collecting milk samples in transport tankers was more economical than in upstream stages. However, Lascano‐Alcoser et al. (2014) and Focker et al. (2019b) showed that it is cost‐effective to monitor dioxins at early stages (e.g., FM) of the food chain, which may be explained by the fact they included tracing costs in the monitoring programs. In terms of using different budgets for different scenarios, the effectiveness of optimal sampling for monitoring AFB1/M1 accounted for the major parts of the total effectiveness, except for monitoring S6 and S8 (Figure 2). The values of decision limits and the sensitivity of analytical methods will also influence the performance of optimal sampling and their optimal settings. For example, a lower sensitivity of both analytical methods (from 100% to 98%) would reduce the DALYs saved from 0.0306 DALYs/100,000 population to 0.0293 DALYs/100,000 population in S1, and a decrease in the decision limits (from 0.05 μg/kg in milk and 2 pg of TEQ/g in milk fat to 0.04 μg/kg in milk and 1.6 pg of TEQ/g in milk fat separately) will improve the performance of optimal sampling to 0.0355 DALYs/100,000 population in S1. Focker et al. (2019b) showed that uncertainty of analytical method increases in analyzing maize samples and the sampling plan would also have a worse performance, which is consistent with our results.

Based on DALYs estimations, human exposure to AFB1/M1 via milk consumption can cause more DALYs as compared to dioxins in most of the scenarios, except for S6 and S8. This could explain why—in the cost‐effectiveness analysis—the optimal sampling reduced more DALYs caused by AFB1/M1 than for dioxins for most of the scenarios (Figure 2). S7 was set as the most common situation with a high contamination fraction for AFB1/M1 and a low fraction for dioxins in the dairy supply chain (Asselt et al., 2017). Previous studies concluded that the presence of dioxins in the global food supply chain is not linked to human deaths, and aflatoxins were associated with most of the DALYs, which is in line with our findings (Boon et al., 2014; Devleesschauwer et al., 2015; Van Kreijl, Knaap, & Van Raaij, 2006). The major part of the budgets was allocated to the MT stage in most optimal monitoring schemes. This is because the concentration of the chemical diluted in the milk from DF to MT, and the number of incremental samples mixed into one na was smaller at MT as compared DF, which resulted into more na at MT. Another reason is that the contamination fraction and the probability of identifying chemicals at MT was highest among all three stages for all scenarios, and the monitoring at MT would reduce more DALYs than other stages.

In the practical, the monitoring situation would be different from our theoretical model in terms of different decision limits and food safety hazards. The decision limits for both chemicals were set according the European legal limits in both feeds and milk, but some food companies may use lower limits in practice and this might influence the results for the optimal sampling settings as described in the sensitivity analysis. Risk managers in food business can change the parameter values of the model to compute the optimal sampling settings for their own business. In the exposure estimation, we assumed a 100,000 population who consumed milk from this dairy supply chain within one year, which can also be transferred to the entire Dutch population. Input values on consumption of milk and body weight of the population were taken from average values from historical data. Our study focused on AFB1/M1 and dioxins because they are the two most relevant chemical food safety chemical groups in the European dairy supply chain, based on their occurrence, toxicity and health impacts, and both of them enter the dairy supply chain through the animal feed intake (Asselt et al., 2017; WHO, 2015b). The model can be easily applied to other food supply chains and can be easily be to other chemicals or other pathogens. For example, except for aflatoxins and dioxins, the FERG's framework has estimated global foodborne disease burdens of 29 hazards, and all of them are available for studying and monitoring food safety at the national level.

The model was calculated in a deterministic way because of data limitations. The exposure estimation and DALYs calculation were built based on available studies and the results of these procedures were used to optimize the sampling settings in the optimization module. Also we assumed that milk truck only gets contamination form one contaminated farm and this would cause the largest contamination fraction in MT stage; however it may be different from the realistic contamination and would cause relatively less collected samples to identify these contamination. It is possible to consider variation and/or uncertainties for other model parameters and to compute results for the full ranges of possible input values (all situations) by using a stochastic approach, for example, Monte Carlo simulation (Adekunte et al., 2010; Van der Fels‐Klerx & Camenzuli, 2016). In our study, sensitivity analyses were done to investigate model results with changing two important model parameters, namely the sensitivity of the analytical method and decision limits. In addition, scenario analyses were performed to obtain insights into the range of outcomes from different initial contamination fractions.

Future research can use the current approach and extend it to optimizing the allocation of resources for monitoring various safety chemicals in food chains. Also, we recommend calculating the shadow price of limited budgets, since that can help risk managers to know the DALYs change with 1 euro additional budget. Food safety authorities and/or food industries can use their own historical data on contamination rates and available resources for monitoring, as well as their own constraints for budgets spent at different food chain stages, such to adapt the model to their case. Our approach was used to estimate risks for human health and food safety, quantify the appropriate monitoring sampling settings to control the risk, and provide stakeholders with the risks and optimal resources allocation for tackling with different chemicals at same time.

CONCLUSION

The risk‐based approach, defined by possible reduction in DALYs caused by chemicals, showed to be appropriate to define a portfolio monitoring sampling for multiple chemicals in the supply chain. In most contamination scenarios, the more focus should be given to monitoring AFB1/M1 in the dairy supply chain, because it reduced more DALYs caused by AFB1/M1 with sampling small amounts of samples. However, when contamination fractions of both chemicals increase, monitoring should give more focus on monitoring dioxins. Overall, more DALYs could be reduced, if the optimal sampling was implemented for higher contamination fraction scenarios within same monitoring budgets.

The results of this research could help food industry and food safety authority to allocate resources in food safety monitoring. The methodology of this research could be used in a wide range of applications to optimize hazard monitoring schemes considering different impacts not only from the perspective of food safety but also from the perspectives of animal health, social welfare, epidemiology, and economics.