R. Woodrow Setzer National Center for Computational Toxicology Office of Research and Development

1. Dose-Response Modeling for EPA’s Organophosphate CumulativeRisk Assessment: Combining Information from Several Datasets toEstimate Relative Potency Factors R. Woodrow Setzer National Center for Computational Toxicology Office of Research and Development U.S. Environmental Protection Agency

2. Background • Food Quality Protection Act, 1996 Requires EPA to take into account when setting pesticide tolerances “available evidence concerning the cumulative effects on infants and children of such residues and other substances that have a common mechanism of toxicity.”

3. Cumulative Risk (per FQPA): • The risk associated with concurrent exposure by all relevant pathways & routes of exposure to a group of chemicals that share a common mechanism of toxicity.

4. Identifying the Common Mechanism Group: OP Pesticides • U. S. EPA 1999 Policy Paper • Inhibition of cholinesterase • Brain • Peripheral Nervous System (e.g., nerves in diaphragm, muscles • Surrogate/indicator (plasma, RBC)

5. Synergy? • Berenbaum (1989) described lack of interaction in terms of the behavior of “isoboles”: Loci of points in “dose space” that have the same response in multi-chemical exposures. • Non-interaction coincides with linear isoboles.

6. Isoboles Example: 2 chems Dose Chem 2 Dose Chem 1

7. Dose-Response for Non-Interactive Mixture For a two-chemical mixture, (d1, d2), if D1 is the dose of chem 1 that gives response R, D2 is the dose of chem 2 that gives response R, then all the mixtures that give response R satisfy the equation: :line For n chemicals: :hyperplane

8. Special Case • When fi(x) = f(ki x): chemicals in a mixture act as if they were dilutions of each other • Isoboles are linear and parallel • Dose-response function for mixture is f(k1x1+k2x2) • Typically, pick one chemical as index (say 1 here) and express others in terms of that. • Then RPF for 2 is k2/k1

9. Strategy of Assessment • Use dose-response models to compute relative potency factors (RPFs, based on 10% inhibition of brain AChE activity: BMD10) for oral exposures; NOAELs to compute RPFs for inhalation and dermal exposures. • Probabilistic exposure assessment, taking into account dietary, drinking water, and residential exposures on a calendar basis. • Final risk characterization based on distribution of margins of exposure (MOE)

10. Vicki Dellarco Elizabeth Doyle Jeff Evans David Hrdy Anna Lowit David Miller Kathy Monk Steve Nako Stephanie Padilla Randolph Perfetti William O. Smith Nelson Thurman William Wooge Plus Many, Many Others OP CRA Science Team

11. Oral Dose-Response Data • Brain acetylcholinesterase (AChE) (as well as plasma and RBC) • Female and male rats • Subchronic and chronic feeding bioassays • Always multiple studies for compounds • Often multiple assay methods • Ultimately, 33 OPs included • Usually ~ 10 animals per dose group/sex • Control CVs < 10%

12. Database of Acetylcholine Esterase Data • 33 chemicals • 80+single-chemical studies • 3 compartments (brain, rbc, plasma)  2 sexes • multiple durations of exposure, subchronic to chronic • total >1655 dose-response relationships (~ 1300 retained)

13. Study 1 Study 2 Study 3 F M F M F M Data Structure Chemical

14. Study 1 Study 2 Study 3 F M F M F M Experimental Design Chemical

15. 7 7 6 28 5 265 4 1077 # Dose Groups 3 1311 2 1312 0 500 1000 1500 # Data Sets Distribution of # Doses

16. Exposure Duration • Preliminary data analysis showed that subchronic feeding studies reached steady state after about 3 weeks • Multiple time points within a study were treated as independent, nested within study. • Only time points with more than 3 weeks of exposure were included.

17. Issues for Modeling • Use as much of the acceptable data as possible • Different units/analytic methods used • Expect responses to differ among compartments, maybe sexes • Generally small number of dose levels in a single data set (limiting the number of parameters that can be identified)

18. Hierarchical Structure of BMD Estimate • Multiple studies carried out at different times by different laboratories, using different analytic methods, reporting results in different units.

19. Two Modeling Approaches: • Model individual data sets, combining estimates. • Model the combined studies for each chemical  compartment. Combined estimate is the estimate of the mean parameter (current revised risk assessment).

20. Modeling Individual Datasets • Fit a model to each dataset, estimating BMD (and estimated standard error) each time. • Model all three compartments and both sexes • Use the global two-stage method (Davidian and Giltinan, 1995; 138-142) twice, once for each level of variability.

21. 2000 1500 1000 AChE Activity 500 potency 0 0 500 1000 1500 Dose Dose-Response and Potency: Approach 1

22. Sequential Approach to Fitting • Fit full model to all data using generalized nonlinear least squares (gnls) • If no convergence or inadequate fit, • Repeat (until good fit or # remaining doses < 3): • set PB 0 • refit to dataset • drop highest dose

23. Potency Measure • Absolute potency is BMD calculated from fitted model: • Relative Potency: • IFPBI = PBk

24. Estimate dose-response for each dataset:

25. Random Effects Model for BMD • Log(BMD) = μlBMD+ EMRID + ETime in MRID • μlBMD varies between sexes • EMRID ~ N(0,σMRID2) • ETime in MRID ~ N(0,σTiM2) • Error variance proportional to (predicted) mean of AChE activity at that dose; constant of proportionality varied among MRIDs.

26. Combine Potency Estimates: • Combine estimates in two stages: among times within study and among studies • At each stage, suppose q individual estimates lmi with variances si2. Potency estimates () and variance components (2) maximize:

27. Combine Potency (more) • Variances for ln(potency) estimates: • This implements the ‘Global Two-Stage’ method of Davidian and Giltinan, (1995) • This method could apply to anysingle statistic or parameter, or vector statistic with simple modification.

28. Problems • Estimate of m depends on PB. Particularly a problem when we cannot estimate PB. • Would like a formal test of whether PBs differ among chemicals. • Is there a shoulder on the dose-response curve in the low-dose region?

29. Solution • Fit the same model to multiple related datasets, allowing information about DR shape to be shared across datasets • Develop a more elaborate model that takes into account some of the biology to give a better description of the lower dose behavior.

30. Q Body ( C ) b b Urine ( k ) e C a Ingestion (Dose ×BW/24) Q Liver ( C ) l l Metabolism ( V , K ) max m Stage 1: A simple PBPK Model • Two compartments: Liver and everything else. • Oral dosing, assume 100% into the portal circulation • Only consider saturable metabolic clearance and first order renal clearance. • Run to steady state

31. Stage 1 (more) • Solve the system of differential equations implied by the model for steady state. • The concentration of non-metabolized parent OP in the body (idose) as a function of administered oral Dose rate is:

32. Stage 2: Same as Before • But reparametrized:

33. DR with First Pass Metabolism

34. Hierarchical Model: • All datasets for a chemical fitted jointly using nlme in R. • S and D varied only among chemicals • A varied among sex × data set • PB varied between sexes • BMD random (same model as before)

35. 15 10 AChE (U/G) 5 0 0.0 0.4 0.8 1.2 Dose (mg/kg/day) Dose Response

36. F M 0.01 0.10 1.00 BMD Benchmark Dose: Fitting One Dataset at a Time

37. F M 0.01 0.10 1.00 BMD Benchmark Dose: Combining Datasets

38. Overall Quality of Fit: Residuals

39. Relative Potencies

40. Computing a MOE (Margin of Exposure) BMD10(A) = 0.08 mg/kg/day MOE = 0.08 X 1000 μg/kg/day / 0.54 = 148

41. Distribution of Total MOEs

42. 1. Combining Estimates • Keeps dose-response modeling “simple” • Delays problems about heterogeneity (sexes, compartments, studies, etc.) until after the modeling. • Number of dose levels in the “smallest” dataset limits the model used, have to drop data sets with too few doses for the selected model.

43. 2. Combining Datasets • Dose-response modeling is (substantially) complicated • Heterogeneity issues addressed in the modeling • Overall number of dose levels (among other things) limits the model used

44. Is PB a High-Dose Effect? • Maybe, but could also be a consequence of multiple binding sites with different functions, or other aspects of the kinetics of AChE inhibition such as variation in aging among chemicals, which could manifest effects at lower doses as well.

45. Horizontal Asymptotes 0.8 Direct Acting Require Activation 0.6 PB 0.4 0.2 0.0 NALED DIAZINON PHORATE PHOSMET FENTHION TRIBUFOS TERBUFOS ACEPHATE ETHOPROP BENSULIDE MALATHION MEVINPHOS PHOSALONE DISULFOTON FENAMIPHOS DICHLORVOS DIMETHOATE FOSTHIAZATE TRICHLORFON DICROTOPHOS METHIDATHION CHLORPYRIFOS METHAMIDOPHOS AZINPHOSMETHYL PIRIMIPHOSMETHYL METHYLPARATHION OXYDEMETONMETHYL TETRACHLORVINPHOS CHLORPYRIPHOSMETHYL

46. Should We Expect Dose-Additivity? (Not Exactly!) • Low-dose shoulder significantly improves fit in a substantial number of chemicals. At best, expect dose-additivity in terms of target dose. • Horizontal asymptotes differ significantly among chemicals (P << 10-6), so dose-additivity cannot hold exactly.

47. Beginnings of A Theoretical Approach • Through mathematical analysis and in silico experiments, ask: • What features determine the shape of individual chemical dose-response curves, and • what are the features of chemicals (if any) that lead to deviations from dose-additivity in cumulative exposures.

48. Example: A “Toy” OP Model • Three compartments: brain, liver, everything else • Constant infusion into the liver • Metabolic clearance in the liver, Michaelis-Menten kinetics: (Vmax, Km) • AChE inhibition in the brain uses same scheme as Timchalk, et al. (2002): Ki, Kr, Ka. • Sample the 5-dimensional parameter space to make example “chemicals”.

49. AChE Inhibition Scheme ks kI ka E + I EI Bound EI kd kr E = AChE I = OP-like inhibitor

50. Strict Sense Dose Additivity