Toxicokinetic Modeling of Per- and Polyfluoroalkyl Substance Concentrations within Developing Zebrafish (Danio rerio) Populations.

Per- and polyfluoroalkyl substances (PFAS) are pervasive environmental contaminants, and their relative stability and high bioaccumulation potential create a challenging risk assessment problem. Zebrafish (Danio rerio) data, in principle, can be synthesized within a quantitative adverse outcome pathway (qAOP) framework to link molecular activity with individual or population level hazards. However, even as qAOP models are still in their infancy, there is a need to link internal dose and toxicity endpoints in a more rigorous way to further not only qAOP models but adverse outcome pathway frameworks in general. We address this problem by suggesting refinements to the current state of toxicokinetic modeling for the early development zebrafish exposed to PFAS up to 120 h post-fertilization. Our approach describes two key physiological transformation phenomena of the developing zebrafish: dynamic volume of an individual and dynamic hatching of a population. We then explore two different modeling strategies to describe the mass transfer, with one strategy relying on classical kinetic rates and the other incorporating mechanisms of membrane transport and adsorption/binding potential. Moving forward, we discuss the challenges of extending this model in both timeframe and chemical class, in conjunction with providing a conceptual framework for its integration with ongoing qAOP modeling efforts.

Introduction

Per- and polyfluoroalkyl substances (PFAS) have moved beyond the emerging contaminant label and are starting to be noticed outside of the environmental community. The prevalence of PFAS use in various industries and the stability and bioaccumulation of PFAS from both an environmental[1,2] and human health standpoint,[3−5] including influences on important healthcare outcomes such as vaccine efficacy[6,7] and COVID-19 symptom severity,[8] creates an immense problem from bioremediation to risk analysis and assessment. On the risk assessment side, given the thousands of compounds that make up the PFAS chemical class,[9] there is a data gap in the literature linking internal dose and toxicity endpoints (including mixture toxicity[10−12]), as pointed out by the recent perspective of Tal and Vogs.[13] However, even for the more well-studied PFAS of perfluorooctanoic acid (PFOA) and perfluorooctanesulfonic acid (PFOS), there is still doubt as to the statistical rigor of internal dose and adverse effect relationships in environmentally relevant organisms, such as the toxicological model zebrafish (Danio rerio). Addressing current limitations of zebrafish toxicokinetic models could strengthen the link between the internal dose and toxicity endpoints, as well as provide a means to extrapolate the effects of well-studied chemicals to those of data poor emerging contaminants.

The early development zebrafish is a fundamental component used in the PFAS risk assessment endeavor.[14−16] It not only provides a robust animal model to work with in the laboratory, such as its high-throughput compatibility,[17] but it also provides a risk assessment model for both the environment (fish in particular) and for potentially humans. In the recent review of Hoffmann et al.[18] they conclude that the early development zebrafish has a capacity to identify mammalian prenatal developmental toxicants. Such a conclusion reinforces the promising human risk assessment capability of the zebrafish as the field waits for more human data. In the meantime, new predictive tools and modeling paradigms can be introduced to (1) enhance the fidelity of temporal and spatial distributions of PFAS within the zebrafish using toxicokinetic modeling and (2) predict toxicity endpoints by linking a PFAS distribution with toxicity data through quantitative adverse outcome pathway (qAOP) modeling.

Adverse outcome pathways (AOPs) and AOP networks are gaining popularity and provide a qualitative framework to organize data across vast length scales (from molecules to individuals to populations) and levels of biological organization (from cells to tissues to organs). Specifically, AOPs conceptualize correlated biological key events (KEs) that describe how the activity of a molecular initiating event (MIE) can be associated with an organism level adverse outcome (AO).[19−23] As an example that we’re using as a case study, exposure of zebrafish to aqueous PFOA has been linked with iodothyronine deiodinase inhibition, which is the MIE of a thyroid disruptor AOP associated with a reduced young of year survival AO.[24] Understanding the toxicokinetics that link exposure with intracellular concentrations could, therefore, inform future new approach methodologies that leverage the chemically agnostic framework of the AOP[25] to fill information gaps in the risk assessment of specific chemicals and substances found in the environment.[26]

In general, toxicokinetic models are important tools for understanding the broad mechanisms of toxicant mass transfer in biological specimens. This mass transfer understanding can, for example, better inform individual or population level predictions of acute or chronic effects. Here, we focus on zebrafish embryos exposed to PFOA and PFOS during early development (0–120 h post-fertilization or hpf) and address some of the limitations of their toxicokinetic models in the literature. A physiologically rigorous toxicokinetic model of zebrafish early development should account for two main phenomena: (1) the dynamic growth of structures associated with embryonic morphogenesis (i.e., surface area and volume) and (2) the asynchronous hatching event and subsequent shedding of the chorion by the zebrafish embryo, resulting in a mixed population of embryos, with and without a chorion over the timescale of about 24 h. Current toxicokinetic models do not adequately account for both phenomena,[14,15,27−31] yet they are especially important for modeling organisms in the environment outside of the laboratory setting where physiological manipulation does not take place (e.g., dechorionation). Furthermore, these models, except for Simeon et al.[27] rely on classical, mostly first-order kinetic expressions to describe the time evolution of a PFAS concentration, creating a phenomenological modeling framework that is difficult to apply outside of its calibrated window of both time and exposure concentration.

In this work, we propose several toxicokinetic model advancements to predict the PFAS distribution between the chorion and embryo of the early development zebrafish, providing a suitable level of mass balance rigor. In our description of mass transfer within the organism, we alter standard toxicokinetic rate terms to incorporate PFAS specific mechanistic processes, such as classical transporter-mediated membrane transport and PFAS adsorption/binding to account for bioaccumulation. This approach could inform next-generation predictive toxicology and risk assessment frameworks (e.g., qAOP) and improve upon existing methods, such as that presented by Vogs et al.[14] which we view as a benchmark for predictive toxicokinetics in the early development zebrafish.

Materials and Methods

Auxiliary Models (Dynamic Volume, Surface Area, and Hatching)

To establish a dynamic volume model relevant for embryo morphogenesis, we used the average data found in a study by Simeon et al.,[27] which describes the volume of both the embryo and yolk with apparent high precision and time resolution across the tested exposure duration. Of note, the degree of exposure-induced malformations, including altered length,[32−34] is a function of PFAS concentration, but the literature has not thoroughly explored other important geometrical metrics, such as volume, which may vary as a function of PFAS exposure. To describe the Simeon et al.[27] data set, we developed a system of first-order ordinary differential equations to model the time evolution of population-averaged tissue volumes. In essence, these differential equations were derived from a mass balance of the system with several assumptions (see graphical abstract for compartment spatial relationships): (1) there is a constant mass-transfer rate from the yolk to the embryo, (2) any loss in mass of the yolk is directly transferred to the embryo, (3) the yolk is the sole mass source of the embryo as the zebrafish is not fed during 0–120 hpf,[35] (4) the chorion has a constant diameter before hatching as supported through image analysis,[36−39] and (5) the yolk has a constant density but the embryo does not. Our volume model derivation is given in the Supporting Information.

Although volumetric information is required for a detailed PFAS mass balance, surface areas are required to estimate mass fluxes necessary to estimate tissue-specific accumulation. We developed a surface area model based on the Guo et al.[40] data set, which can be found in the Supporting Information.

Moving on from dynamic volume and surface area, we also modeled the dynamic hatching of a population of developing zebrafish. Hatching is an individual transformative process with the kinetic rates underlying the biochemistry having some distribution, which leads to different hatching times across a population, on the order of about a day.[32] As such, hatching is not instantaneous, and any population measurement (i.e., one that is a function of the hatched state) made after the first individual has hatched will reflect an average of the mixed state. This is especially important for the toxicokinetics discussed in this work as the presence or absence of the chorion (i.e., the development status of the population) greatly affects the mass transfer. In general, a mixed-state solution to a problem also solves the uniform population problem when it arises (i.e., the fraction of one state is set to zero). For our purposes, whenever the chorion is present (i.e., any natural hatching scenario), we have a mixed population, but for dechorionation, we have a uniform population. To account for hatching, we developed a hatching fraction model based on the Hagenaars et al.[32] data set. Unlike our volume model, hatching fraction dependence on the PFAS concentration is recorded, but we did not consider this dependence due to simplicity and limited data. The hatching fraction model can be found in the Supporting Information.

Toxicokinetic Model Development

We adopted a compartment-based abstraction of the zebrafish embryo, in which PFAS transport across tissue or cell types can be modeled by a set of chemical “reactions” that regulate PFAS accumulation within each compartment. This approach reduces the complexity of molecular transport across different membrane barriers and biological fluids (e.g., blood, interstitial fluid, and cytosol) to a small set of kinetic rate constants. This simplicity comes at the cost of broad reductive assumptions. First, we assume the number of PFAS molecules is large enough that concentration fluctuations can be safely neglected. When tissue concentrations are small enough, this rate equation approach breaks down, and corrections are required.[41] Second, we must acknowledge that a rate-limited deterministic chemical kinetics model allows for instantaneous molecular transport between cell and tissue types, which is unrealistic, even in principle. Despite these potential drawbacks, the rate equation approach has been employed before to model zebrafish embryo toxicokinetics. Specifically, we consider the model described by Vogs et al.[14] (referred to hereon as the “Vogs model,” including citation) as a benchmark of predictive PFAS toxicokinetics in the early development zebrafish.

We obtained PFOA and PFOS data sets presented in the Vogs model and adopted their well-mixed two-compartment representation of chorion (subscript C) and embryo (subscript E). Here, “chorion” refers to the perivitelline space, and “embryo” refers to a combined compartment of both embryo and yolk tissues. This chorion/embryo system is exposed to a reservoir of constant aqueous PFAS concentration (subscript W). In the following model equations, we also use subscript X which can either be W or C depending on the hatched state.

Our model differs significantly from the Vogs model. First, we adjusted the mass balance of the system, treating the PFAS mass (or number of molecules) as a conserved quantity. This is appropriate as both a changing compartment volume and mass flux across a compartment lead to a change in the compartment concentration. In contrast, the Vogs model treats the PFAS concentration as a conserved quantity. As such, we converted the concentrations reported by the Vogs model back to original mass measurements by using a volume model described in a study by Brox et al.,[28] which the Vogs model used for their calculations. Also, considering the bioaccumulation process of some PFAS and the lack of metabolism of PFOA and PFOS,[42] we hypothesized that elimination rates are much smaller in value than those of absorption (uptake) and are thus negligible, which is consistent with the rate constants reported by the Vogs model.

One limitation of the Vogs model is that it relies on a different set of parameters for each PFAS exposure concentration. As such, there is no unified description of a physical process responsible for such a concentration dependence, leading to a lack of predictive capability in its present state. We found that saturable kinetics, as modeled with Michaelis–Menten kinetics, can describe this concentration dependence. Finally, we scaled mass fluxes to the growing surface area of the relevant fluid/fluid or tissue/fluid interface. A mass balance for the chorion is thenwherein C is the PFAS concentration (μM), V is the volume (L), A is the surface area (mm2), t is the time (h), k is the maximum kinetic rate constant between two locations (μmol/h/mm2), and K is a concentration for which the flux is half its maximum. As VC is known (governed by our established dynamic volume model), eq is used to numerically solve for CC as a function of time.

Equation describes the time evolution for accumulated PFAS concentrations within the chorion of an unhatched embryo. However, embryos will hatch over some time course as we have discussed with our hatching fraction model in the Auxiliary Models subsection. In the Vogs model, each measurement represents a homogenization of multiple zebrafish, meaning that each zebrafish could theoretically be in a different exposure state depending on when it hatched. We account for a mixed population of hatched and unhatched embryos by modeling the fraction of an embryo population that has hatched at a given time. As a result, the chorion mass contribution for a given measurement is directly modulated by the hatched fraction and can be expressed aswherein Mmeas refers to the mass measurement (μmol) and h is the fraction of hatched embryos.

For the embryo compartment, this hatching phenomenon needs to be reflected at the mass balance level because it receives two different fluxes: one from the reservoir and another from the chorion. Both fluxes change over time and are dependent upon the hatched state of the population. We address this problem by writing a mass balance for the average change in mass of the embryo compartmentwherein Cavg labels a population-averaged PFAS concentration. Here, we have assumed that individual embryos are independent, identical copies distributed across the hatched or unhatched states, with their volume governed by our established dynamic volume model. The embryo mass contribution to a given measurement can, therefore, be calculated asonce eq is used to numerically solve for CEavg as a function of time. Equations and 4 can then be added together and compared against the Vogs model data for parameter fitting.

Modeling Membrane Transport and PFAS Adsorption/Binding Processes

Common practice in toxicokinetic modeling is to abstract the toxicant mass-transfer process across biological compartments into that of reaction-limited chemical kinetics. However, this approach neglects the underlying mechanisms (e.g., diffusive flux processes) responsible for the evolution of a compartment’s toxicant concentration. For many substances, this phenomenological-based abstraction is adequate. However, PFAS exhibit complicated molecular conformations (e.g., micelle formation), and because they are amphiphilic, they exhibit unique interactions with biological structures such as lipid membranes and proteins. We hypothesize that these nontrivial interactions require a more mechanistic modeling approach to better capture and understand the toxicant mass transfer. To test this hypothesis, we adjusted the mathematical framework described above by considering PFAS diffusion across the chorion (a membrane) and outer embryo membrane (i.e., abstracting the embryo boundary as a single membrane and the only resistance to mass transfer, such as the apical membrane of an epithelium) as rate-limiting PFAS transport mechanisms, in addition to PFAS adsorption/binding events within each compartment.

We tested a model of simple diffusion across the boundaries of both compartments (i.e., the mass flux across either the chorion or embryo membrane is only governed by a constant permeability and a concentration gradient), which did not fit the data well at any time. At early times, if simple diffusion dominated the physics, we would expect the mass flux to be proportional to the concentration gradient, and the absence of this lends credence to the idea that PFAS molecules may interact with membrane transport proteins that can become saturated. This results in a nonlinear concentration dependence of the mass flux, and we hypothesize that membrane transporters, such as organic anion transporters, interact with and limit PFAS transport across membranes. This idea is justified on the basis that membrane transporters are generally important in pharmacokinetics[43−45] and PFAS toxicokinetics.[46,47] Perhaps one of the simplest representations of carrier-mediated transport is that of Kolber and LeFevre,[48] wherein a solute–carrier complex forms and diffuses across the membrane, and due to the finite number of carriers, the transport flux rises with solute concentration until all carriers are bound.

In the literature, zebrafish-relevant reports span a number of diverse adsorption/binding events, from nanoparticles[49,50] to PFAS,[51−53] and the data directly support the presence of such events as we observed gross underestimation of system mass at long times for our simple diffusion model. Our consideration of adsorption/binding mechanisms assumes diffusion-limited transport across membranes. In conjunction, we identified the simplest adsorption/binding model consistent with the data (i.e., number of parameters), which led to the use of a linear isotherm to describe the general sum of adsorption (or partitioning) to lipid membrane surfaces and binding to other structures (e.g., proteins) within each compartment. As such, we developed a simplified representation of transporters moving PFAS across the membrane boundary of a compartment followed by the instantaneous adsorption/binding of the PFAS to structures within the compartment, of which the potential adsorption to the membrane does not influence the number of available transporters.

As with our previous model formulation, we estimate the mass contribution of the chorion outside of its mass balance (i.e., h appears in the mass contribution equation rather than the mass balance, see eqs and 2), but we must consider a mass balance for unhatched embryos. We assume that PFAS fluxes across the system membranes are driven by concentration gradients; therefore, fluxes associated with unhatched embryos require an embryo concentration term. Additionally, we describe two PFAS phases for each compartment: the phase in free solution (superscript f) and the phase bound to biological structures (superscript b). The relevant mass balance for the chorion is

Here, the subscript un refers to the unhatched embryos, M is the mass, Kads is the parameter associated with adsorption/binding affinity (L), Ksc is the constant of the solute–carrier complex (μM), and k now represents the constant of diffusional transport across the membranes between two locations (μmol × μM/h/mm2). The chorion mass contribution to a given measurement can, therefore, be expressed as:

For the embryo compartment, we again incorporate a two-phase representation of PFAS

The embryo mass contribution to a given measurement is then

Model Training and Parameter Identification

We employed a curve-fitting methodology to identify model parameters for both model forms considered. Parameter values minimize a modified least squares objective functional developed with data leveraged from those described by the Vogs model. Details of the fitting procedure are provided in the Supporting Information.

Results and Discussion

Empirical Assessment of Model Performance

Auxiliary models developed to support the toxicokinetic models (i.e., dynamic volume, surface area, and hatching) exhibit good qualitative agreement with data in all cases (R2 > 0.97, see Figure S1). The R2 for all toxicokinetic models (model code found in Supporting Information) using the calibration data is > 0.96, with slightly higher fidelity exhibited by models that incorporate PFAS transporters and adsorption/binding. Other goodness of fit metrics, such as the RMSE, show similar results. Thus, a higher level of mechanistic detail seems to improve the description of PFAS toxicokinetics in conjunction with providing a much deeper mass-transfer understanding, albeit at the cost of two additional parameters. From a solely information metric perspective, the Akaike Information Criterion bias-corrected for small sample sizes (AICc) justifies the two additional parameters for PFOA but not for PFOS, but the Bayesian Information Criterion (BIC) justifies the two additional parameters for both PFOA and PFOS.[54] All model fits are presented in Figure S2 (including all goodness of fit and information metrics), as well as our replication of the Vogs model. Also, the Vogs model equations and parameters (see Table S2) are compiled in the Supporting Information.

In Figure , we present the transporter and adsorption/binding model (referred to hereon as the “transporter model”) fits in a parity plot form for both the Vogs model data and other data found in the literature with no experimental setup restrictions.[16,42,55−59] As expected, the Vogs model data, which we used for calibration, cluster around the unity line more so than the other data found in the literature (evaluation data), as can be seen by a difference in the calculated standard deviation of the data about the unity line (i.e., the scaling factor of the data from unity). The area bounded by a factor of 2 from unity is shaded, and ideally, as reported by the International Programme on Chemical Safety,[60] toxicokinetic models should generally yield estimates within a factor of 2 of the data if they are to be used in risk assessment. Overall, for the Vogs model data comparison, we see that 90.0% and 85.7% of the data fall within a factor of 2 for PFOA and PFOS, respectively. For the comparison against the broader literature, 12.5% and 71.4% of data fall within a factor of 2 for PFOA and PFOS, respectively. It is not entirely clear why the PFOA data are so variable in comparison to PFOS, but we compare (see Table S3) and comment on some differences between the experimental studies in the Supporting Information.

Figure 1

Parity plot presentation of the transporter model predictions against that of both the Vogs model data and other data available in the literature[16,42,55−59] (incorporating both pre- and post-hatching phases) for the chorion/embryo system in the concentration domain. When necessary, data were converted to a mass/zebrafish basis before being converted to a concentration basis using our proposed volume model. We assumed a wet weight of 500 μg[58,59] for the Han et al.[42] conversion. The shaded region indicates an area bounded by a factor of 2 from unity, and the standard deviation lines of the data from unity are included.

Hatching Significantly Affects PFAS Toxicokinetics

The Vogs model suggests that it is important to model the chorion due to its role as a mass-transfer barrier. Although the chorion is, indeed, a physical barrier to PFAS flux into the embryo, it is not clear whether it is also a rate-limiting step for PFAS toxicokinetics or accumulation within embryo tissues. To justify inclusion of the chorion, the Vogs model points to a dynamic, biphasic PFAS uptake pattern, as demonstrated by an increased PFAS accumulation in the period just after hatching, which suggests a net-positive uptake rate for embryos without their chorion. However, our reanalysis of these data does not support a definitive conclusion regarding a biphasic pattern in the data (see Figure S3). In particular, it is difficult to distinguish positive or negative trends from fluctuations in data, which is not unusual given a limited number of data points, especially considering that net uptake is a derived measurement (numerical derivative) from the internal mass measurement. Given limitations in the data, a modeling argument should be made to support claims of biphasic uptake.

In Figure , we show our model results for embryo flux (net uptake) compared to those of the Vogs model. In contrast to our model, wherein we show a distinct increase in PFAS accumulation around the hatching transition for all exposure conditions, the Vogs model does not, instead predicting an abrupt decrease in embryo accumulation once chorions are shed for some conditions. This stark, qualitative difference in predictions can, perhaps, be traced back to how each approach incorporates hatching. Hatching is, approximately, an abrupt physiological change for an individual zebrafish embryo, but the onset of individual hatching can vary over many hours. At the population level, this leads to a situation where some embryos have hatched while others have not, and since PFAS accumulation is measured based on the homogenate of multiple individuals, chorionated individuals will contribute less PFAS accumulation than dechorionated ones. In our transporter model, this imbalance can most readily be seen for PFOA, where a distinct increase in accumulation for the population begins at the onset of hatching and continues as the hatched fraction increases. In contrast, the hatching transition of the Vogs model treats all individuals identically, and a population average therefore adopts the abrupt transitional symmetry of the individual hatching event. Thus, our transporter model predicts a longer timescale to resolve the long-term behavior of the system after hatching. For example, at the end of the hatching period for PFOA, we begin to see a decrease in the embryo flux in contrast with PFOS. This is due to the greater bioaccumulation potential of PFOS (as evident from parameter differences such as the adsorption/binding constant in Table S1), while PFOA begins to approach chemical equilibrium for lower exposure concentrations, as discussed later.

Figure 2

Transporter model results for embryo PFAS net uptake (flux) compared against that of the Vogs model. Fluxes from the Vogs model were estimated by multiplying their change in the concentration and volume predictions.

The “Vogs model” is not a single predictive model, but rather a set of models, each with kinetic parameters specific to a tested exposure condition. Unfortunately, this constrains the ability of the Vogs model to interpolate between measured data and to extrapolate beyond fitted conditions of exposure and duration. As shown in Figure , we extrapolate model predictions (including volume) past 5 days post-fertilization (dpf) to that of 10 dpf. In the Vogs model, PFOA and PFOS exhibit a roughly constant accumulation beyond the hatching event, where a constant uptake flux dominates over the elimination flux, resulting in a relationship for PFAS mass that is approximately linear in the exposure duration. However, our transporter model predicts that lower concentration exposures of PFOA will approach chemical equilibrium after approximately 5 dpf. This steady state is something that the Vogs model cannot predict with constant PFAS accumulation at a given exposure concentration; thus, its fidelity for extrapolation is greatly reduced. The nonlinear accumulation characteristics predicted by our transporter model begs further investigation of PFAS bioaccumulation as we expect PFAS with higher bioconcentration factor (BCF) values (e.g., PFOS) to rapidly accumulate as fluid phase concentrations increase (as captured in our model) but limited by the number of adsorption/binding sites that emerge as the fish matures. The ability to capture such long-term behavior as chemical equilibrium and potential accelerated accumulation puts our model in the unique position for full life cycle fish studies.

Figure 3

Transporter model results compared against that of the Vogs model for chorion/embryo system mass after extrapolating all models past 5 dpf to 10 dpf for three different exposure concentrations.

Overall, our hatching representation provides a rigorous way to deal with the shedding of the chorion from a population of fish either in the laboratory or environmental setting. With this basis, we can address any setting where the chorion is or is not present. For instance, as shown in Figure S4, we compare transporter model predictions for the chorion and dechorionated cases. We see that the chorion acts as a barrier to mass transfer (as shown before in our biphasic uptake discussion) since a greater uptake is seen for the dechorionated cases, resulting in greater system mass. Also, these model predictions indicate that these PFAS have a greater affinity for the embryo over the chorion. If these PFAS had a greater affinity for the chorion over the embryo, we would see the chorion curves show a greater mass than those of the dechorionated curves. In conjunction, if the affinity for the chorion was dominant, we would also see a decline in system mass in the Vogs model data set after the hatching event. We do not see this, which is different than the results reported by Brox et al.[61] for other compounds outside of PFAS (e.g., clofibric acid and mitribuzin), where greater chorion affinity and a decrease in mass after hatching is observed, which points to the chorion acting as more of a mass sink for these compounds. As such, if the mass-transfer role of the chorion can be characterized, such as with our modeling framework, the need for laboratory dechorionation is eliminated. The internal concentration responsible for a given MIE as a function of the exposure concentration and time can always be deduced, regardless of the presence of the chorion.

In Figure S5, we show a compartment-by-compartment analysis of both PFOA and PFOS mass in each phase represented by our transporter model. Given the bioaccumulative nature of some PFAS, we see how the majority of PFAS mass resides in the embryo bound phase. Such a result aligns with the conclusion that the chorion is more of a mass-transfer barrier rather than sink for the PFAS investigated. Also, for PFOS, we see a greater percentage of system mass reside in the bound embryo phase compared to that of PFOA. This result highlights the BCF difference between the two PFAS as will be discussed later but also speaks to the much greater toxicity of PFOS when considering the same toxic effect,[14] perhaps due to a much greater adsorption/binding potential.

Regardless of chorion versus embryo affinity, the chorion acts as a mass-transfer barrier and mass sink, with one phenomenon having the ability to be dominant depending on the chemical of interest. Overall, the parameters of our model (i.e., adsorption/binding constants) can be tuned to account for either barrier or sink dominant situations, allowing our model to be extended beyond PFAS applications. Furthermore, our proposed model provides a more physiologically relevant framework for not only developing zebrafish in the laboratory that could be dechorionated but also fish that hatch naturally in the environment.

PFAS Bioconcentration in a Population of Growing Embryos

Figure illustrates the exposure concentration–response curves for both PFOA and PFOS, as calculated from our transporter model and compared against exposure data reported by the Vogs model. In the transporter model, sigmoid concentration dependence of the relevant transfer fluxes model a PFAS fluid and bound phase phenomenology that manifests with a concentration–response curve that is insensitive to higher PFAS exposures. Here, a competition between PFAS influx/efflux predicts thresholds of insensitivity on the order of 1000 μM and 100 μM for PFOA and PFOS, respectively. Future experiments can test these predictions, and discrepancies could point to transporter and/or membrane integrity disruption not accounted for in our model, in conjunction with a breakdown of the infinite adsorption/binding site assumption used here. Furthermore, more data could motivate more complicated adsorption/binding kinetic and isotherm models that capture the nuance of long-range electrodynamic interactions between polar PFAS molecules and the tissue environments.

Figure 4

Chorion/embryo system PFAS mass expressed over a range of exposure concentrations and times predicted using our transporter model (curves). Vogs model data are shown with symbols.

Our model can also be used to investigate both the temporal and concentration dependence of a PFAS accumulation factor (AF). Here, we define the AF to be a time-dependent quantity that is mathematically synonymous with the BCF, but we avoid referring to the AF as the BCF, as the BCF is normally defined around chemical equilibrium (i.e., a time-independent quantity). The AF encodes important information that is relevant for the MIEs of AOPs that begin well before chemical equilibrium is reached, and the literature also provides some single AF (referred in the references as BCF) values for the early development zebrafish not at chemical equilibrium.[14,42]Figure illustrates that the time evolution of the AF tracks that of system mass (which is requisite considering the definition of the AF), highlighting the need to exercise caution when interpreting BCFs reported in the literature for just a single time point not at chemical equilibrium. Shown here are AF predictions at several different exposure concentrations, contrasted against values reported or derived from the literature.[16,42,55−59] As Figure shows, and as Chen et al.[62] discussed for long-chain PFAS, BCFs can decrease as the exposure concentration increases, which is consistent with our model that incorporates PFAS adsorption/binding and saturable membrane transport.

Figure 5

Our transporter model predicts some nonmonotonic AF profiles that evolve over time (curves), illustrated here for several different constant exposure concentrations. Here, we calculate AF as the ratio of the chorion/embryo system concentration to that of the constant exposure concentration. AF data (open symbols) reported or derived from the literature[16,42,55−59] are also shown with a similar calculation strategy as with Figure . Closed symbols (as depicted in the legend) represent our model predictions for those same data points, whereby giving a sense for the magnitude of the literature exposure concentrations and our model’s fidelity with the literature in the AF domain.

Figure shows that AF predictions for PFOS can be up to about an order of magnitude larger than AF predictions for PFOA, consistent with other studies that present an AF (referred in the references as BCF) at 120 hpf.[14,42] As the Vogs model suggests, internal concentrations of several PFAS are a better predictor of AOs than exposure concentrations. With respect to the AOP framework, internal concentrations are relevant at the level of the MIE, which is correlated to an AO via biological KEs. Thus, if two chemicals act through the same AOP, then knowledge of the whole AF profile could improve or refine the fidelity of the toxicity prediction. Even if two chemicals do not act through the same AOP, knowledge of the whole AF profile can still provide some toxicokinetic learnings separate from being used as an input into two different AOPs. As an example, long-time behavior of the PFOA AF is nonmonotonic for lower exposure concentrations, where it decreases from shortly after the hatching event to the end of the exposure experiment. We hypothesize that for smaller exposure concentrations, an increase in volume of the developing zebrafish dominates the PFOA flux into the zebrafish, acting to dilute tissue concentrations and decrease the AF. Since environmental exposures are much lower than that used in most studies, the physiology of the embryo, such as tissue growth, could play a much larger role in the toxicokinetics when compared to the flux of a chemical into the organism.

Overall, our proposed modeling framework provides a solid foundation for the continuing evolution of toxicokinetic modeling of PFAS in zebrafish by addressing the rigor of system volume changes, hatching, and mass transfer. Given the generality of these modeled phenomena, we could potentially apply this model to embryos of other fish species and potentially to a larger number of chemicals or chemical classes given data availability. As AOP development and qAOP modeling become more prevalent, it becomes more of a necessity to capture as much physiological detail and as many specific chemical transport mechanisms and interactions as possible to provide a reasonable internal concentration estimate at the level of the MIE. As such, future toxicokinetics steps would include adding more compartments to the model to allow more MIEs to be addressed that require more spatial resolution, as well as extending the timeline that the model covers, which would need to include new uptake (and perhaps elimination) mechanisms as the physiology progresses.[52,63,64] In the literature, it has been shown that separating yolk from embryo greatly increases the accuracy of embryo concentration predictions,[65] and such a simple separation can also shed light on AOs that would require enhanced modeling fidelity to distinguish yolk from embryo tissues.[66] Ideally, radiolabeling[67] and tomography could be used to track the PFAS temporal and spatial distribution across the lifespan of a zebrafish and generate multi-compartment data to support improved model development.

For our future qAOP model, we propose an approach that includes modular “response–response” KE linkages, potentially as described by Foran et al.[68] and Song et al.,[69] yet also models the mean and variance of a response across a population experiencing a MIE, KE, or AO (in a similar vein to the multistate model of Simeon et al.[70]). More precisely, we propose that Bayesian networks[71] can parameterize statistical distributions that capture relationships between the MIE, KEs, and AO that are conditioned over a wide range of PFOA exposure concentrations. Thus, we anticipate that suitable application of Bayes’ theorem to every linkage of the AOP will produce conditional distributions that are chemically agnostic, making the final model relevant to a wide range of PFAS and other compounds. We believe that the free solution PFOA concentration predicted by our transporter model is an appropriate concentration to use for predicting the activation of the MIE (the qAOP model input) for our AOP case study.