&EPA
United States
Environmental Protection
Agency
EPA/600/R-16/062 | May 2016
www.epa.gov/ord
Matrix Population Model for Estimating Effects
from Time-Varying Aquatic Exposures:
Technical Documentation
1.2 -,
1.0
0.8
0.6
0.4
0.2
0.0
10
15
Year
20
25
30
Glen B. Thursby
Atlantic Ecology Division, NHEERL, ORD
US Environmental Protection Agency
Narragansett, Rl 02882
Office of Research and Development
National Health and Environmental Effects Research Laboratory
-------
Notice
The research described in this report has been funded wholly by the U.S.
Environmental Protection Agency. This report is contribution number ORD-016393 of
the Atlantic Ecology Division, National Health and Environmental Effects Research
Laboratory, Office of Research and Development. This document has been subjected
to USEPA's peer review process and has been approved for publication. The mention
of trade names of commercial products does not constitute endorsement for use.
Abstract
The Office of Pesticide Programs (OPP) models daily aquatic pesticide exposure
values for 30 years in its risk assessments. However, only a fraction of that
information is typically used in these assessments. The population model employed
herein is a deterministic, density-dependent periodic matrix model for integrating
time-varying pesticide exposure effects on the marine invertebrate Americamysis
bahia. The external exposure concentrations are converted to time-varying scaled
internal concentrations by coupling a one-compartment toxicokinetics-
toxicodynamics model with the matrix model. Several exposure scenarios (each with
the same risk as determined by OPP's traditional approach) were created within
which population modeling documented different risk conclusions than assessments
based on the traditional approach. Population modeling incorporates all available
toxicological and exposure data, making a more complete assessment of the
potential risk of time-varying aquatic concentrations.
Keywords: Americamysis bahia, matrix modeling, population level risk assessment,
time-varying exposures
Acknowledgements
Dr. Keith Sappington from the US Environmental Protection Agency's Office of
Pesticides Programs provided valuable insights into the application of these results
for pesticide risk assessments. Drs. Giancarlo Cicchetti, Jason Grear and Diane Nacci
provided technical reviews of the final report. Other comments from Marty Chintala,
Joseph LiVoIsi and Dr. Matthew Etterson also were helpful. Any remaining errors or
misinterpretations of comments are the responsibility of the author.
-------
Preface
Evaluation of different aquatic exposure scenarios is an integral part of EPA's risk
assessment process for the registration or re-registration of pesticides. While the
Office of Pesticide Programs (OPP) has research needs associated with exposure
models (e.g, incorporation of spatial variability in exposure parameters or
development of urban and residential aquatic exposure models), there is also a need
to better link time-varying exposure to population-level effects. The successful
application of population modeling would significantly reduce the uncertainty
associated with OPP's risk assessments. Such modeling can account for cumulative
impacts of multiple "low" exposures to the same chemical (minimizing potential for
under protection), and can account for recovery after a short "high" exposure
(minimizing potential for over protection).
One of the criticisms of population modeling is the limitation placed on this approach
by its data requirements. Current toxicity data requirements for pesticide registration
do not provide enough information for complex models. The specific model herein
was created with this concern in mind. All toxicity parameters are derived from
standard test data for the marine invertebrate Americamysis bahia. The model
successfully integrates acute and chronic toxicity data, and uses all of OPP's time-
varying modeled exposure data.
This report is a product for CSS 18.04.6, Integrated Modeling for Ecological Risk
Assessment, Task 6: Case Study 3 - Pesticides Impacts to Aquatic Endangered Species.
-------
Contents
Notice ii
Abstract ii
Acknowledgements ii
Preface iii
INTRODUCTION 1
The Issue 1
The Potential Solution 2
Risk Analyses 4
MODEL STRUCTURE 6
Matrix Population Model 6
Effects Modeling Constraints 8
Estimating Time-Varying Survival Probability 8
Estimating Time-Varying Reproduction Probability 9
Minimum Data Requirements 10
Model Operation 10
TOXICOLOGICAL FACTORS 13
RESULTS AND DISCUSSION 16
Deterministic 16
Probabilistic 17
Population Modeling 17
Acute "Kinetics" 18
Survival vs Reproduction Effects 20
Spawning Season 22
CONCLUSIONS 24
References 25
Appendix A. Daily Modeled Pesticide Concentrations 27
Appendix B: Explanation of Periodic Model and Density Dependent Factor 28
Survivorship Calculations 31
Reproduction 34
Appendix C: Calculation of Elimination Constant 35
Appendix D: Explanation of Mancini Calculation of Scaled Internal Concentration 38
Appendix E: Mysid Chronic Data for Endosulfan and Model Calibration (Quality Control)
40
IV
-------
INTRODUCTION
For assessing pesticide risks to aquatic organisms, the Office of Pesticide Programs
(OPP) models pesticide spray drift, runoff, and erosion into an agricultural pond with
specified water body and watershed characteristics. Numerous agricultural crop
scenarios represent different crop, regional, climate, watershed, and agronomic
specifications across the country. These crop scenarios are intended to capture
agronomic and regional factors that greatly influence the delivery of pesticides to
surface waters (e.g., precipitation patterns, soil characteristics, pesticide application
timing). With each of these crop scenarios, an exposure model called the Pesticide in
Water Calculator (PWC)1 simulates daily concentrations for 30-year exposure
distributions for surface water, sediment, and interstitial (pore) water. Because the
PWC output depends in part on soil properties, soil and crop management practices
and weather data, different regions of the country and different crop types will have
different aquatic exposure time series for equivalent applications of the same
pesticide. Typically, OPP derives a single acute and chronic estimated environmental
concentration (EEC) from each 30-year time series that represents a conservative
(high end) exposure concentration with an infrequent occurrence (e.g., once in 10
years). The acute and chronic EECs are then compared to available acute and chronic
toxicity data to estimate risk to aquatic organisms. OPP's use of EECs reflects a tiered
process whereby high end estimates of exposure are first used to identify potential
risks, which identifies any need for additional refinements to risk estimates. It is clear
that the EEC distills a large amount of information into a single value and thus greatly
limits the ability to consider temporal aspects of organism exposure and life history
which are known to influence risk. Population modeling, while not currently part of
OPP's standard risk assessment, offers an additional higher tiered refinement—one
that can use all of the exposure data, as well as combine acute and chronic
assessment into a single estimate of risk. The use of population modeling assessment
procedures will reduce the likelihood of over or under protection decisions that may
occur due to an inability to incorporate time-varying exposures in the risk assessment
process.
The Issue
OPP models daily aquatic pesticide exposure values for 30 years (10,350 values), yet
only a fraction of that information is typically used in pesticide risk assessments.
Time-varying exposure information is essentially not used. Depending on the type of
analysis, the daily time series data may be used as is or converted to a 4-d running
average (for comparison against acute toxicity data), or, if chronic toxicity data are
used, converted to either a 21-d (invertebrates) or 60-d (fish) running average. Once
1 http://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/models-pesticide-risk-
assessmenttfaquatic
-------
a time series is selected, the maximum value for each year is noted. The EEC is
determined as the 90th percentile of 30 ranked values of annual maxima and is
subsequently used to calculate a risk quotient (RQ), which is an EEC divided by a
toxicity value of interest. In theory, this EEC reflects the annual maximum with an
expected 1 in 10 year return frequency and roughly corresponds to the fourth highest
annual maximum. In the case described in this report, the toxicity value is a mysid
chronic value. Because such a small fraction of the available exposure data is used,
one can easily imagine a situation whereby several time series have very similar
EECs—thus very similar RQ values—yet the underlying pattern of daily exposures
could be drastically different. In this situation, the currently estimated risk to the
taxonomic group of interest would be similar, but the actual potential risk might be
quite different among the different exposure series. Figure 1 shows a hypothetical set
of exposure scenarios. The three scenarios were constructed so that each has the
same third highest annual maximum. Visual inspection of these scenarios suggests
they should have different effects on a population; however, the current risk
assessment method would indicate the same risk. The challenge is to provide a
procedure that can distinguish among these scenarios, yet be easy to understand and
simple to implement. Ideally the procedure also would require no new data—that is,
take full advantage of all of the currently required toxicity data for a given species.
The Potential Solution
One technique for incorporating effects of time-varying concentrations of pesticides
is toxicokinetics-toxicodynamics modeling—TK-TD (Brock et al. 2010, Ashauer and
Brown 2013). Toxicokinetics relates the time course of the concentration of a toxicant
within an organism to the time course of that toxicant in the external medium. TK
modeling is similar among many researchers. Most publications on the subject make
the simplifying assumption that the organism is a one-compartment model2 (i.e., the
whole body concentration represents the target concentration), and assume uptake
and elimination kinetics are first order (linearly related to external and internal
concentrations, respectively). The biggest objection for using toxicokinetics has been
the need for internal concentration data—not typically available from standardized
tests. This can be overcome by using a scaled internal concentration (Ashauer and
Brown 2013), which does not require estimation of the uptake kinetic parameter, and
may even be the preferred method for determining uptake kinetic parameters over
whole body measurements (Jager et al. 2011). Scaled internal concentration is
described below in more detail in the MODEL STRUCTURE section.
2 This may be because much of the TK-TD modeling to date has been with relatively small
invertebrates.
-------
0.8
,0.6
0.4
0.2
0.0
21-d Running Average: Low
0
1.2 n
1.0
0.8
,0.6
0.4
0.2
0.0
0.8
B>0.6
0.4
0.2
0.0
10
15
Year
20
21-d Running Average: Medium
10 15 20
21-d Running Average: High
10
15
Year
20
25
25
25
30
30
30
Figure 1. Three different time series, each with the same 30 annual maximum values. Time series differ only in the
concentrations represented by the non-maximum values. For the "Low" time series, each of the "Medium" values
were divided by 2, and for the "High" series, each of the "Medium" values were multiplied by 1.5. The plots are
21-d running averages and the original raw daily values are plotted in Appendix A.
-------
While toxicokinetics is similar among many different TK-TD modeling efforts,
toxicodynamics modeling can take several forms. Most of the TD approaches are
different applications of hazard modeling3—as the internal concentration increases,
the probability of an effect increases. The approaches for toxicodynamics differ
primarily in the component to which the toxic effect (usually survival) is related—that
is, whether it relates directly to the internal concentration or to some level of internal
damage (Brock et al. 2010). In its simplest form the hazard rate is proportional to the
internal concentration, and all concentrations, no matter how small, cause some
degree of effect. This simple model also assumes any effect (or any recovery) is
instantaneous and complete. The model described in this report uses this simplified
form for toxicodynamics. Because of this simplification, the model can be
parameterized using the toxicological data typically available to OPP during pesticide
registration or re-registration.
Even though TK-TD models usually are applied to individuals, they can be directly
coupled with population models (Brock et al. 2010), which is what is done herein. For
this population approach, the report builds upon a matrix population model
previously developed for the marine invertebrate Americamysis bahia (Thursby
2009). That model has a one week time step4. This new approach takes that into
consideration and uses a periodic matrix technique (Caswell 2001), which creates 52
sub-matrices5 to represent annual population activity. This method allows
incorporating potential effects due to degree of seasonal exposures overlapping with
seasonal biological events, such as reproduction. The purpose of this report,
however, is not to promote a definitive population model, but rather to use this
model to demonstrate the ability of such models to distinguish among various
exposure time series, as well as among different toxicological features—such as
seemingly minor differences in sensitivity.
Risk Analyses
A comparison of three different types of risk analyses demonstrated how population
modeling provided a more quantitative distinction among different exposure time
series. These analyses were deterministic, probabilistic, and population-level
methods. The deterministic approach is currently the first tier of an assessment
within an OPP regulatory risk determination. As mentioned above, this method uses a
very small portion of a 30-year modeled exposure data series. The probabilistic
approach is occasionally included in a risk assessment and includes all of the
exposure data within a cumulative distribution. The same toxicological endpoints are
3 Survival is considered a stochastic process. This is in contrast to an individual tolerance approach
whereby there is the assumption of a distribution of thresholds for survival. In the latter approach, the
individuals that die are the more sensitive ones; in the stochastic approach, the individuals that die are
just the unlucky ones (Jager et al. 2011).
4 The earlier model also is stochastic and density independent. The periodic version is not stochastic
and is density dependent. Details of the model are in Appendix B.
5 Note, this results in a 364-day year.
-------
used as in the deterministic approach, except this approach calculates the probability
that the exposure data exceeds that endpoint (based on a count of how many
concentration values are greater than the endpoint value). The probabilistic approach
makes no distinction about when a concentration occurs, it only determines how
many times within 30 years that concentration occurs. For example, it does not
distinguish between a high concentration occurring thirty times within a single year
and the same concentration only occurring once a year.
As with the probabilistic approach, the population approach uses all of the exposure
data—integrating potential time-varying effects on survival and reproduction into a
single population response tracked through time. Population modeling also
incorporates recovery when exposure concentrations decline. The approach
eliminates the need to make separate acute and chronic risk determinations. In
addition, it eliminates the need to rely on running averages because daily
concentrations in the external medium are used to estimate daily scaled internal
concentrations, which in turn track daily effects6. As such, and unlike the probabilistic
method, the population modeling method can distinguish among exposure scenarios
that may differ only by when a given concentration occurs. Finally, population
modeling can evaluate other exposure-based issues, including such things as
frequency of stressful exposures and the temporal distribution of stressful exposures
(e.g., clumping vs. uniform spacing of exposure within a time series).
Population models can also distinguish among effects that differ from a toxicological
perspective. For example, the relative ratio of 24 vs. 96 hr LC50 values. With the
current pesticide risk assessment procedures, two species with similar 96 hr LC50
values (a standard acute assessment endpoint) for a given toxicant would be deemed
to have similar acute risks. However, for some species the LC50 may stabilize after 24
hr, and for others, it may continue to change with increased duration of exposure.
These are two different toxicokinetic parameters, and very likely result in two
different population responses to the same exposure time series. Another
toxicological factor is the relative response of survival vs. reproduction. Two species
could have similar final chronic values based on reproduction, but the effect of the
toxicant on survival could be very different. Again, two similar risks by current
procedures, but likely very different population responses.
This report shows toxicokinetics-toxicodynamics modeling and periodic matrix
population modeling are useful tools to account for a variety of exposure and
toxicological factors. This assessment is based, in part, on the mysid's response to
endosulfan, a compound already being removed from all uses within the United
States.
6 Each weekly survival probability is the minimum survival rate of that seven-day period—a mysid can
only die once so the minimum within each week is selected.
5
-------
MODEL STRUCTURE
Matrix Population Model
As stated above, the model builds upon an earlier matrix population model
developed for Americamysis bahia (Thursby 2009). A matrix model allows
independent tracking of effects on different subpopulation groups. This adaptation of
the model retains the basic structure of the original matrix model—subpopulation
groups are age classes and a one week time step. All age classes (ranging from 1 to 13
weeks) are assigned the same sensitivity, differing only in the length of time exposed.
The earlier model is density independent, stochastic, and assumes constant exposure
concentrations. The current model is density dependent7 and deterministic8. For
specifics on the derivation of control survival and reproduction demographics, refer
to the earlier report.
The original mysid matrix model construct is insufficient to deal with time-varying
toxicant concentrations. Because OPP uses modeled 30-year daily time series for
estimating exposure to aquatic organisms, a different matrix configuration is needed.
To accomplish this the original weekly time step was retained, but a separate matrix
was created for each week of the year. This approach was patterned after periodic
matrix models (Caswell 2001)9. In addition to providing a mechanism to track changes
in effects due to variability in exposure through time, this approach also provided a
means for incorporating variability in demographic parameters with time (e.g.,
spawning and non-spawning seasons). Figure 2 demonstrates pictorially how the
model was constructed. Time-series exposure data occur one year at a time, and the
beginning of subsequent years starts where the previous year left off until all 30 years
are processed. Each matrix represents the population's status for its given week;
however, each matrix retains a "memory" of its previous 12 weeks of the exposure
series (see Appendix B for a more complete explanation). The endpoint was the
proportion of weeks within the 30-year time series that the population declined to or
below a given threshold based on weekly counts of total population size.
7 Populations in the natural environment do not grow indefinitely. There must be some sort of
compensatory factor governing overall population growth. The current model uses density
dependence as this factor. The mechanisms for density dependent growth are varied, and generally
not specifically known for most species. Many population modelers default to one of several simplistic
ways to incorporate density dependence. For this report, I have chosen to use the density dependent
mechanism described by Leslie (1948).
o
The demographic parameters (survival and reproduction) are fixed, changing only when exposed to a
toxicant. Therefore, if you run the same exposure time-series more than once, the distribution of
population size over time will be the same each time.
9 Often, periodic matrix models are presented as a product of a sequence of sub-matrices which
themselves represent different periods of time within an annual cycle—for example, spring, summer,
fall and winter. The model presented in this report uses sub-matrices, but does not multiply them
together for a single annual matrix.
-------
x q X
Spawning Season
150 200
Time (da)
Figure 2. Diagram of how the population model interacts with the toxicant time series.
The solid black line represents a hypothetical time series of toxicant concentration in
the water column. The gridded squares across the top represent the positioning of the
periodic sub-matrices in time. Note that there are only 11 matrices pictured for
simplicity. In the model there are 52 such matrices covering each week of the year. The
red line represents the spawning season—user defined.
The original mysid matrix model relies on Population Viability Analysis (PVA) to
provide the effects endpoint (Ginzberg et al. 1982, Morris and Doak 2002).
Population Viability Analysis estimates the probability that a population will fall
below a given threshold within a given period of time (e.g., the probability of a 20%
decline from the current population size within the next 10 years). Initially, the intent
was to use a similar stochastic process to model the annual growth rate based on the
periodic matrix—then apply PVA procedures. However, PVA was not practical after
the analyses were expanded from one matrix to fifty-two10. Even very small changes
in demographic parameters within the control matrices resulted in large variations in
the annual growth rate of the population. For example, a weekly survival as high as
0.99 translated into an annual survival rate of 0.59 (0.99 raised to the 52nd power).
Because of this, the periodic version of the matrix model was simplified to
deterministic rather than stochastic—this also eliminated the need for hundreds, or
even thousands, of runs for each annual series11.
In addition, PVA analysis assumes that the mean population growth rate is density-independent.
Incorporating density dependence in the model complicates estimates of extinction risk (Morris and
Doak 2002).
11 PVA analysis requires not only an estimate of an average growth rate for a time period of interest,
but also an estimate of the variability of that rate due to stochasticity of demographic parameters—
thus the need for hundreds or even thousands of model runs.
-------
The 2009 version of the model does not track population numbers, only changes in
growth rate—so that model uses a simplifying assumption of no density dependence,
that is, it allows exponential growth. The periodic model tracked population size on a
weekly basis for 30 years. Exponential growth made tracking population size
impractical since extremely large population sizes occurred. Density dependence
within the model eliminated the possibility of these excessive population sizes. In
addition, because the periodic model tracked exposures which vary over time, a
mechanism for incorporating recovery was needed when exposure concentrations
declined. Density dependence was one way to accomplish this12. The model applied
the density dependent factor equally across each age class. Appendix B describes
how this factor was calculated.
Effects Modeling Constraints
The Office of Pesticides Programs has specific test guidelines for acute and chronic
tests for the mysid shrimp Americamysis bahia (OPPTS 850.1035 and OPPTS
850.1050). At least for the foreseeable future, ecological risk assessments will likely
rely on the current suite of traditional tests. The approach in this report, therefore,
restricted the parameter estimations for TK-TD modeling to data typically expected
from these standard toxicity tests. This limited some of the types of TK-TD modeling
to couple with the periodic matrix ensemble. Most notably damage and
recovery/re pair models were not considered. Damage and recovery models can
account for delays in effects and different rates of recovery. By definition, damage is
a reduction in the fitness of the organism (Brock et al. 2010), and the effect (usually
survival) is proportional to the amount of damage. The rate of repair or recovery (also
proportional to the amount of damage) has to be taken into account in "damage" TD
models. However, the calculation of recovery rate constants requires measured
internal concentrations (Jager et al. 2011). Recovery, by definition, also has to be
related to sub-lethal endpoints—an individual cannot recover from mortality.
Endpoints such as growth and reproduction often do not have enough time
dependent observations to efficiently estimate organism recovery. For these reasons,
the model herein did not consider a damage and recovery form of TK-TD modeling.
Estimating Time-Varying Survival Probability
To translate the effects of time-varying exposure concentrations into time-varying
survival probabilities, there needs to be a way to estimate how the internal
concentration changes with time. Often this is accomplished assuming a one-
compartment model (Kooijman and Bedaux 1996a)—which treats any internal
concentration as being uniformly distributed within the organism. The concentration
in an organism at any given time depends on an increase based on the concentration
in the water times the uptake rate and a decrease based on the internal
12 The actual mechanism of density dependence is not modeled. For example, the approach does not
distinguish between whether the dependence is a result of intra-specific competition for food or
increase rate of predation as the population increases.
8
-------
concentration times an elimination rate13. The elimination rate can be estimated
from data that is frequently available from traditional toxicity tests. If the LC50 is
calculated for different time intervals (e.g., 24, 48, 72 and 96 hr), the LC50 often
shows an exponential decay in time with a decay constant that can serve as an
estimate of the elimination rate constant14 (Bass et al. 2010). See Appendix C for a
full explanation. The other rate constant (uptake rate) is not so easily estimated,
since concentrations within the animal's body are generally not available.
Kooijman (1983) and Mancini (1983) independently solved the issue by using what is
now referred to as a scaled internal concentration (Jager et al. 2011). Kooijman
(1983) scaled the unknown internal concentration by the bioconcentration factor
(uptake rate divided by elimination rate) and related mortality to this value. Mancini
(1983) scaled the internal concentration by dividing it only by the uptake rate
constant, since it is possible to independently estimate the elimination rate constant
based on LC50 vs. time. Mancini's explanation of the mathematics is easier to follow,
and his explanation is apparently the first to apply the model to exposure scenarios
where the concentrations both decrease and increase over time. A relationship
between the scaled internal concentration and % mortality forms the basis for
tracking survival probability. The slope of this relationship is what Kooijman and his
colleagues more recently refer to as the "killing rate" (Kooijman and Bedaux 1996a,
b; Jager and Zimmer 2012). Mancini's procedure forms the basis for the use of scaled
internal concentration in this report—a full explanation is provided in Appendix D.
The calibration procedure for fine tuning this rate is described in Appendix E, using
chronic toxicity test results for endosulfan.
Estimating Time-Varying Reproduction Probability
Mancini (1983) and others provided ways to evaluate survival probability using data
from standard toxicity tests—much of it relying on acute data. Simple TK-TD models
for sublethal endpoints are not easily calibrated and generally require more data than
provided by traditional test protocols (Ashauer and Brown 2013). Standard toxicity
tests often do not have sufficient time series data for reproductive output in order to
estimate directly the kinetic coefficient for reproduction (i.e., a reproductive "killing
rate"). The most readily available reproduction data will be chronic end-of-test
effects information. For the purpose of the model described herein, the ratio of
chronic survival effect (e.g., 28-d LC50) to chronic reproduction effect (e.g., EC50) was
used. A Survival to Reproduction Ratio (SRR) was calculated and the assumption
13 The elimination rate accounts for many different processes, including actual excretion, internal
binding, internal transformation to another chemical form, etc. Mancini (1983) uses the term
detoxification rather than elimination.
14 The elimination rate constant may not just represent whole body elimination. It is possible that a
portion of the compensating process may represent repair or recovery's rate constant. Some
researchers, therefore, refer to this as the "dominant rate constant" (e.g., Jager et al. 2011). For the
purpose of the current model I will assume elimination constant refers to the combined actions that
eliminate toxicity.
-------
made that this ratio remains constant for any probability of survival. The
reproduction rate factor (RRF) is related to survival probability as: RRF =
EXP[LN(SP)*SRR]. Where SP is the survival probability for a given week. The model
multiplies the control maternity rates by the RRF to estimate the maternity rate for a
given exposure week. The RRF is assumed to be constant for all age classes. The
calibration procedure for fine tuning the SRR, if needed, is described in Appendix E.
Minimum Data Requirements
The period matrix model toxicity data requirements are the same as the minimum
data requirements for the earlier, non-periodic version of the model (Thursby 2009).
Thursby (2009) provides a procedure for estimating default values even if these
minimum data are not available. Briefly, the minimum data requirements are:
1. Acute LC50 values for different time periods—typically these will be 24, 48, 72
and 96 hr. These data are used to calculate the elimination rate constant. This
constant is the only variable needed to convert time-varying external
concentrations to time-varying scaled internal concentrations.
2. Survival probability over time for one or more "constant" concentrations.
These data are often available from daily observations of mortality and are
used to estimate the killing rate—which is assumed to be a constant for any
concentration. The killing rate converts daily scaled internal concentrations to
a daily survival probability. Calibration of the model output using chronic data
can refine the killing rate constant.
3. Chronic LC50/EC50 ratio (SRR). This ratio is used to estimate the weekly
reproduction adjustment factor.
Model Operation
After the killing rate constant, the uptake rate constant, and the SRR are derived, the
model operation is straightforward. An exposure time-series is selected and the
spawning season defined. For this report each week was assigned a spawning factor
of either 1 or 0. A factor of 1 meant that the original fecundity rates were used. A
factor of 0 meant that all fecundity rates were set to zero. The default spawning
season was 39 continuous weeks and began at week 10 (week 1 was the first week of
January). For scenarios evaluating spawning season effect, six different scenarios
were run—a 39-week season starting either week 10, 20 or 30, and a 26-week
season, also starting either week 10, 20 or 30. Figure 3 demonstrates the first three
years of a control population with spawning beginning on week 10 and lasting for 39
weeks. At the beginning of each modeling run, the model runs for several "years"
with no exposure. This allows the population to reach stability with respect to annual
cycling of weekly total population size, as well as stability with respect to the
distribution of individuals among the different age classes. For each series with a
different spawning scenario a separate 30-yr control was created. The modeling
result is a time series showing the weekly change in population size as a percentage
10
-------
of the control response. Figure 4 presents a partial time series as an example of
how data are evaluated. Results are shown from the first ten years. The
response is quantified by counting the number of times the weekly population
size falls below a given threshold—expressed as a fraction of the total number of
weeks in 30 years (1560).
Control
100
Year
Figure 3. Example of a control population time series with
spawning beginning in week 10 and lasting for 39 consecutive
weeks. Only the first three years are shown—all 30 years are
identical.
The risk of a population falling below a threshold obviously is a function of a specific
threshold (the smaller the change from the control, the greater the potential for
observing that change). A problem with this approach has been the selection of a
population threshold (Thursby 2009). This can be overcome by the use of risk curves
in which a range of population thresholds is used and the area under such curves
calculated (Burgman et al. 1993)15. An alternate approach could use already
established thresholds for different degrees of severity in population decline. For
example, the World Conservation Union (IUCN 2012), defines a population as
vulnerable if a 30% decline is observed over a specified amount of time or number of
generations. A population is endangered if there is a 50% decline, and critically
endangered if an 80% decline is observed or estimated16.
15 Risk definition in Burgman et al (1993) uses quasi-extinction, which is not the same as what is
presented herein; however, the concept of total risk being related to area under the curve is similar.
16 The World Conservation Union thresholds are for threatened and endangered species. They are
presented here only as examples of how one might summarize the severity of risk to a population.
11
-------
100
c
o
1/1
c
o
Q.
O
Q.
20 --
Figure 4. Sample model output showing population size as % of control for the
first ten years of a model run. The horizontal dashed lines for the population
thresholds for vulnerable (70%), endangered (50%) and critically endangered
(20%). The number of weeks a population falls below a given threshold
compared to the total number of weeks in 30 years (1560) is a direct estimate
of the susceptibility of the population.
12
-------
TOXICOLOGICAL FACTORS
Besides the above biological scenarios, several toxicological model scenarios also
demonstrated the utility of population modeling. Figure 5 shows two sets of LC50
dynamics with time. These could represent two different species or two different
toxicants with the same species. The point addressed is the effect of the rate at which
mortality occurs as an organism responds to constant exposure. Both scenarios have
the same 96-hr LC50—thus the traditional risk assessment would assign the same
acute risk. However, one scenario clearly responds quicker to exposure (dashed line)
than the other (solid line).
LC50 dynamics
5.0 i
4.5 -
4.0 -
3.5 -
5 3.0 -
D)
^
$ 2-5'
o
2.0 -
1.5 -
1.0 -
0.5 -
0.0 -
Slow Kinetics
Fast Kinetics
Figure 5. Acute toxicity scenarios. One in
which the 24 hr and 96 hr LC50 values are very
similar (dashed line) and the other where the
LC50 changes more slowly with time (solid
line). The elimination rate constant of the
former is 1.25 d"1 and the latter, 0.025 d"1. See
Appendix C for a detailed explanation of the
rate constant calculation.
2 3
Time (da)
A population model can distinguish between these acute scenarios. Figure 6
demonstrates why. The bottom panel displays a hypothetical 30-year time series of
daily exposure concentrations. The top panel displays the daily probability of survival.
The species whose LC50 value is similar at 24 and 96 hr ("fast kinetics") tracks more
closely the peaks of daily exposure. The "slow kinetics" species buffers the exposure,
such that sudden daily increases in exposure pass by before the effect of that
exposure reaches its full potential. See Appendix D for a more detailed explanation.
13
-------
Dailv Probability of •Survival
10
15
20
25
3D
Jciilv -Kpostrc 'values
Year
Figure 6. Daily survival probabilities (top) for two hypothetical species
exposed to a daily exposure time series (bottom). See Figure 5 for description
of fast and slow kinetics.
A second set of toxicological scenarios is shown in Figure 7. These hypothetical data
represent dose response data from chronic tests. The chronic runs both used the
same effect on reproduction, therefore both had the same NOAEC (no observable
adverse effect concentration)—based on reproduction. The only difference was the
relative sensitivity of survival compared to reproduction. In one scenario (Figure 7,
top) survival was assumed to be insensitive relative to reproduction. In the other
scenario, survival had a similar dose response to that for reproduction (Figure 7,
bottom). Because both scenarios have the same NOAEC, both would result in the
14
-------
same evaluation of chronic risk. Yet, logic would clearly dictate these two scenarios
should not have the same level of risk. Again, population modeling will distinguish
between these two chronic scenarios.
2 3
Pesticide Concentration
Figure 7. Chronic toxicity dose response
scenarios. The scenario in the upper plot
shows an NOAEC based on reproduction,
and survival is significantly less sensitive
than reproduction. The lower plots shows
the same NOAEC; however, this time
survival has a similar sensitivity to the
pesticide as reproduction.
2 3
Pesticide Concentration
15
-------
RESULTS AND DISCUSSION
Risk assessments using deterministic, probabilistic and population modeling
techniques show the value of population modeling for exposure time series. The goal
of this work was not to promote a particular model, but to promote the value of
population modeling when specifically applied to risk assessments for pesticide time-
varying exposures. The information presented below demonstrates the value of
population modeling to distinguish among exposure and effects where the traditional
assessment resulted in the same estimate of risk for each scenario. The analyses are
based on the 21-d running averages presented in Figure 1 (daily values for each are in
Appendix A).
Deterministic
All three time-series have the same annual maximum concentrations, and thus the
same 90th percentile value. The summary of these values for the 21-d running
averages (Figure 1) is presented in Table 1. The 90th percentile (0.772 ug/L) is used to
calculate the risk quotient (RQ). For this example, the effect concentration is 0.49
ug/L17, making the RQ 1.58, which is above the LOG of 1 for aquatic invertebrates.
Based on this, we know there is the potential for acute effects on marine
invertebrates—and that potential is the same for all three time series (by design).
Table 1. Annual maximum values from Figure 1.
Year
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Maximum
(ug/L)
0.147
0.185
0.754
0.759
1.135
0.468
0.585
0.724
0.532
0.402
0.531
0.623
0.453
0.663
0.522
Year
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
Maximum
(ug/L)
0.462
0.572
0.691
0.631
0.377
0.478
0.971
0.598
0.890
0.419
0.468
0.561
0.625
0.324
0.575
17 This is the mysid chronic value for endosulfan. A different chronic value would be used for
freshwater invertebrates. A 60-day running average would be used for evaluating chronic effects to
freshwater and saltwater fish.
16
-------
Probabilistic
If we look at a cumulative distribution of the data from Figure 1, then differences
among the data are more easily quantified (Figure 8). Although the deterministic
approach shows the same risk analysis for each of the three exposure scenarios, the
probabilistic approach clearly shows differences among the three. The probability of
exceeding the chronic value for mysids (0.49 ug/L) is 3, 5 and 17%, for the low,
medium and high series, respectively. However, we have no easy way to determine if
any of these exposure scenarios are "bad" enough for concern.
21 Day Running Average
1.0 •
0.8
0.6
0.4
0.2
0.0
0.0 0.2 0.4 0.6 0.8
Concentration (ug/L)
1.0
1.2
Figure 8. Cumulative distributions of the exposure data presented in Figure 6.
The vertical dashed red line is the chronic value for mysids exposed to
endosulfan (0.49 ug/L).
Population Modeling
Population modeling results are presented two ways, as the chance of decline below
one of several thresholds, and as the total area under a risk curve. Figure 9 shows
population results for endosulfan exposure which can be directly compared to the
deterministic (Table 1) and probabilistic results (Figure 8). Whereas the deterministic
and probabilistic methods only used the single chronic value for effect of endosulfan
on mysids (0.49 ug/L), the population approach incorporates all of the dose response
information for both acute and chronic exposures. Using all of the toxicological
information available offers a more complete distinction among the three exposure
scenarios. These distinctions are summarized in Table 2.
17
-------
Endosulfan
0 10 20 30 40 50 60 70 80 90 100
% Decline from Control Response
Figure 9. Summary of the chance of a given % decline in mysid population
size relative to the control response for each of the three exposure scenarios
when the exposures are assumed to be from endosulfan. Spawning season
began on week 10 and lasted for 39 weeks. See Appendices C and E for
summary of the endosulfan toxicity data. Vertical dashed lines represent %
declines corresponding to IUCN categories of vulnerable (30%), endangered
(50%) and critically endangered (80%).
Table 2. Summary statistics for Figure 9. AUC refers to the area under the curve.
Exposure
Scenario
Low
Medium
High
30% Decline*
0.30
0.59
0.80
Threshold
50% Decline
0.07
0.25
0.44
80% Decline
0.00
0.02
0.13
AUC
20.4
35.8
48.9
* Decline data are fraction of time below the threshold.
Acute "Kinetics"
Figure 10 displays the population modeling results from the comparison with
different acute toxicity dynamics—the acute dose responses for all model runs had
the same 96 hr LC50. Model parameters assumed no direct effect on reproduction,
and spawning season began on week 10 and lasted 39 consecutive weeks. The only
difference among the runs were the values for elimination kinetics constant (see
18
-------
Figure 5 for example, and Appendix C for explanation of kinetics constant). These
ranged from 0.01 to 1.25 d"1 ("slow" to "fast" kinetics). Only the data from the
"medium" exposure time series are displayed. The vertical dashed lines correspond
to the declines associated with a population being labeled as vulnerable, endangered
and critically endangered (IUCN 2012). The data for each kinetic scenario are
summarized in Table 3. The area under each risk curve corresponds to the total
expected relative decline in the population size. It is worth reminding that for all of
these kinetic scenarios the risks determination by the deterministic and probabilistic
methods are the same—the differences observed in Figure 10 are cause by
alterations in the kinetics of acute toxicity only. The more similar a 24hr LC50 is to a
96 hr LC50 the quicker a species mortality rate responds to a change in the
environmental concentration. The greater the ratio of 24 to 96 hr LCSOs, the slower
the mortality rate responds to changes in external toxicant concentrations and the
more a species response is buffered against sudden, short-term changes in daily
exposures.
LC50
10
20 30 40 50 60
% Decline from Control Response
Figure 10. Summary of the chance of % decline in population size relative to the
control response for ten different LC50 kinetic senarios. All population model
runs used the medium exposure time series (see Appendix A). Elimination
kinetic constants ranged from 0.01 d"1 (lowest curve) to 1.25 d"1 (upper most
curve). Table 3 lists all of the kinetic contants, along with the summary data for
each curve. Spawning season began on week 10 and lasted for 39 consecutive
weeks. Vertical dash lines represent % declines corresponding to IUCN
categories of vulnerable (30%), endangered (50%) and critically endangered
(80%).
19
-------
Table 3. Summary statistics for Figure 10. Kinetic constant refers to the value for the
elimination constant. Only data for the medium exposure series are shown. Percentage
declines correspond to IUCN categories of vulnerable (30%), endangered (50%) and critically
endangered (80%). AUC is the area under the risk curve.
Kinetic
Constant
0.010
0.025
0.050
0.075
0.100
0.250
0.500
0.750
1.000
1.250
30% Decline*
0.00
0.11
0.26
0.35
0.43
0.59
0.64
0.65
0.66
0.66
Threshold
50% Decline
0.00
0.00
0.06
0.12
0.16
0.25
0.29
0.30
0.31
0.32
80% Decline
0.00
0.00
0.00
0.00
0.00
0.02
0.04
0.05
0.05
0.06
AUC
3.67
12.19
20.66
25.74
28.95
35.83
38.42
39.48
40.11
40.49
*Decline data are fractions of time below the threshold.
Survival vs Reproduction Effects
Figure 11 is similar to Figure 10, except instead of displaying the results from
different acute toxicity scenarios, it displays the results from two different chronic
toxicity scenarios. All six model runs had a spawning season beginning on week 10
and lasting 39 consecutive weeks. In the first scenario (the upper plot in each of the
three graphs), reproduction and survival have similar dose-responses to constant
exposure concentrations (see Figure 7, bottom panel). In the second scenario for
each time series, reproduction is significantly more sensitive to those exposures than
survival (see Figure 7, top panel). The vertical dashed lines correspond to the declines
associated with a population being labeled as vulnerable, endangered and critically
endangered (IUCN 2012). The data are summarized in Table 4. As expected, the
probability of decline for both kinetic scenarios increases as the exposure increases
from low to high. Predictably, the model runs where both survival and reproduction
are effected have a higher probability of decline than those where essentially only
reproduction is influenced. Similar to the model runs for Figure 10, each pair of data
is associated with the same probabilistic exposure summary.
20
-------
(U
(U
o
o
0.4 -
0.2 -
0.0 -
Low
Repro = Survival
Repro > Survival
10
20
30
40
50
60
70
80
90
100
Medium
20
30 40 50 60 70
% Decline from Control Response
80
90
100
Figure 11. Summary of the chance of % decline in population size relative to the control response for
each of the three exposure time series (see Appendix A). For each time series, two different
toxicological scenarios were run. For the first (red line), the dose-response of reproduction and
survival were set to be equal. For the second (blue line), reproduction was kept the same as the first
run; however, survival sensitivity was significantly less than that of reproduction. The elimination
constant was 0.25 d"1 (slow kinetics from Figure 5). Vertical dash lines represent % declines
21
-------
corresponding to IUCN categories of vulnerable (30%), endangered (50%) and critically endangered
(80%).
Table 4. Summary statistics for Figure 11. Low, Medium and High columns are data from three
different exposure scenarios (see Figure 4). Repro = Survival corresponds to the model run where
reproduction and survival dose-response were the same. Repro > Survival is the model run where
reproduction was significantly more sensitive than survival (see Figure 3).
Threshold
Dose-Responses
Exposure Scenario
Low Medium
High
Fraction of Weeks Below Threshold
30% Decline
50% Decline
80% Decline
Total Risk
Repro = Survival
Repro > Survival
Repro = Survival
Repro > Survival
Repro = Survival
Repro > Survival
Repro = Survival
Repro > Survival
0.53
0.19
0.26
0.02
0.03
0.00
Area
34.5
15.2
0.84
0.42
0.59
0.18
0.21
0.01
Under Curve
56.6
28.3
0.91
0.63
0.82
0.32
0.52
0.06
75.5
39.5
Spawning Season
The effect of spawning season start date or length was not as dramatic as the effect
of acute or chronic toxicity scenarios (Figure 12). Only the results for the "medium"
time series are presented. There is little difference among the spawning start weeks
(10, 20 or 30). It is worth noting, however, that the order of decline associated with
the three different start weeks is different within the two different spawning season
lengths. Also, the spawning season lasting 26 weeks generally had lower probability
of decline values than the 39 week season.
Figure 13 displays the results from model runs using two different weekly population
growth rates (lambda). Both runs used the medium exposure time series. The dashed
lines corresponds to a run in which the weekly maximum growth rate was 1.6118—
the growth rate for a control population during spawning season used in all previous
model runs. The solid line represents a model run in which the growth rate of the
control population was reduced to 1.41. This was accomplished by multiplying the
spawning rates for each spawning season matrix by 0.6. The potential of a population
decline increases substantially during an exposure time series when the underlying
growth rate is smaller. This is likely for longer-lived species whose weekly growth
rates should be substantially lower than those calculated for the short-lived mysid.
This is the control's weekly growth rate using the mean demographic parameters from Thursby
(2009).
22
-------
M-jd'u n 39
Figure 12. Summary of the effect of length and timing of spawning season
on the % decline in population size relative to the control. Only the
"medium" exposure scenario (see Figure 5) results are shown. The upper
set of plots are for a spawning season lasting 26 consecutive weeks,
beginning week 10, 20 or 30. The lower set are for a spawning season
lasting 39 consecutive weeks.
Medium
i.o
0.9 -*
0,8
0.7
0,6 -
0.5 -
0.4
0.3
0,2
0.1
0.0
> = !.«!
»-!-!> I.
10 20 .» 40 SO 60 70
% Decline from Control Response
80
100
Figure 13. Summary of the effect of weekly population growth rate on the
% decline in population size relative to the control. Only the "medium"
exposure scenario results are shown. The length of spawning season was
39 weeks, with spawning begin on week 10.
23
-------
CONCLUSIONS
Population models hold great promise for integrating exposure, toxicity and life
history information into meaningful measures of risk. Several circumstances were
presented within which population modeling would result in different risk
conclusions than assessments based on the traditional approach. The traditional first
tier procedure to risk evaluation (deterministic) relies only on the annual maxima.
Three exposure scenarios were presented where the total effect of 30-year time-
varying exposure series would have very different risks to a population, yet the
deterministic assessment indicated the same risk for all series. The probabilistic
approach was a little better. That approach at least uses all of the exposure data, and
showed each exposure time series was distinct. This approach, however, can only
indicate the proportion of time the expected environmental concentration exceeds
the toxicological endpoint in question. The probabilistic method provided no
information on the expected total consequence of those exceedences. Population
modeling, on the other hand, provided an approach that allowed those consequences
to be quantified. Once quantified, it becomes easier to determine whether the
estimated effect is unacceptable or not. Clearly, population modeling provides a
more complete assessment of the potential risk of a time-varying exposure than the
traditional assessment methods.
The population model used in this report is certainly not the only mathematical
model that could have been employed. The specific model used, however, was
selected based on several constraints. These included having sufficient demographic
data on which such a model could be based, using a commonly tested species, and
being able to derive toxicity parameters for the model based entirely on standard test
data. While other models may prove eventually to be more appropriate, the general
conclusions about the relative value of population modeling should still stand.
The acute and chronic scenarios presented clearly indicate the potential for missing
significant effects using only traditional endpoints such as LCSOs, NOAECs and RQs.
The model runs using different spawning seasons, however, did not show as great an
effect on a mysids population's response to a given time series. The preliminary
model runs using a reduced weekly growth rate suggest this lack of an obvious effect
could be because mysid populations in general have a rapid growth rate. A short-
term decline from a large exposure concentration can be rapidly overcome with a
few weeks of lower concentration. It is likely that species with a significantly slower
population growth rate—such as longer-lived invertebrates or fish—may display a
more pronounced spawning season effect.
The selection of IUCN Red List thresholds for decline are consistent with their
assessments for vulnerable, endangered and critically endangered species. These are
presented as examples of how population modeling output might be evaluated based
on observations of population size over 30 years. Red List categories can be avoided
by using the area under a risk curve to establish the effect of a time series on a
population. This, however, does not eliminate the need for judgement. Someone still
24
-------
has to decide where to place the cut off between an acceptable and unacceptable
area. Ultimately, the significance of any difference is a policy decision.
References
Ashauer, R and CD Brown. 2013. Highly time-variable exposure to chemicals—toward an
assessment strategy. Integrated Environmental Assessment and Management.
9(3):e27-e33.
Baas, J, T Jager and B Kooijman. 2010. Understanding toxicity as a process in time.
Sci.TotalEnviron. 408:3735-3739.
Brock, CM, A Alix, CD Brown, E Capri, BFF Gottesburen, F Heimbach, CM Lythgo, R Schulz
and M Streloke. 2010. Linking Aquatic Exposure and Effects: Risk Assessment of
Pesticides. CRC Press. 410 pp.
Burgman, MA, S Ferson and HR Akcakaya. 1993. Risk Assessment in Conservation
Biology. Chapman & Hall. 314 pp.
Caswell, H. 2001. Matrix Population Models: Construction, Analysis, and Interpretation.
2nd Ed. Sinauer Associates. 722 pp.
Ginzberg, LR, B Slobodkin, K Johnson and AG Bindman. 1982. Quasi-extinction
probabilities as a measure of impact on population growth. Risk Analysis 2:171-
181.
IUCN. (2012). IUCN Red List Categories and Criteria: Version 3.1. Second edition. Gland,
Switzerland and Cambridge, UK: IUCN. iv + 32pp.
Jager, T, C Albert, TG Preuss and R Ashauer. 2011. General unified threshold model of
survival—a toxicokinetic-toxicodynamic framework for ecotoxicology.
Environmental Science & Technology 45:2529-2540.
Jager, T and Zimmer, El. 2012. Simplified dynamic energy budget model for analyzing
ecotoxicity data. Ecological Modelling. 225:74-81.
Kooijman, SALM. 1983. Statistical aspects of the determination of mortality rates in
bioassays. Water Research. 17(7):749-759.
Kooijman, SALM and JJM Bedaux. 1996a. Analysis of toxicity tests on Daphnia survival
and reproduction. Water Research. 30:1711-1723.
Kooijman, SALM and JJM Bedaux. 1996b. The Analysis of Aquatic Toxicity Data. VU
University Press. 149 pp.
http://www.bio.vu.nl/thb/research/bib/KooyBeda96.html.
Lee, ET. 1992. Statistical Methods for Survival Data Analysis. John Wiley & Sons. 482 pp.
Leslie, PH. 1948. Some further notes on the use of matrices in population mathematics.
Biometrika. 35:213-245.
25
-------
Mancini, JL. 1983. A method for calculating effects, on aquatic organisms, of time
varying concentrations. Water Research. 17(10):1355-1362.
McKenney, CL, Jr. 1982. Final report for the interlaboratory comparison of chronic
toxicity testing using the estuarine mysid (Mysidopsis bahia). US EPA, Gulf
Breeze, FL. Report to S. Ells, Health and Environmental Review Division. Office of
Toxic Substances. US EPA. 35 pp.
Morris, WF and DF Doak. 2002. Quantitative Conservation Biology: Theory and Practice
of Population Viability Analysis. Sinauer Associates. 480 pp.
Schimmel, SC. 1981. Results: Interlaboratory Comparison—Acute Toxicity Tests Using
Estuarine Animals. US EPA Report # EPA-600/4-81-003.
Sprague, JB. 1969. Measurement of pollutant toxicity to fish—I. Bioassay methods for
acute toxicity. Water Research. 5:793-821.
Thursby, GB. 2009. Americamysis bahia Stochastic Matrix Population Model for
Laboratory Populations: Technical Documentation. EPA Report # EPA/600/R-
09/121. 73 pp.
Widianarko, B and N Van Straalen. 1996. Toxicokinetics-based survival analysis in
bioassays using nonpersistent chemicals. Environmental Toxicology and
Chemistry. 15(3):402-406.
26
-------
Appendix A. Daily Modeled Pesticide Concentrations.
3.0 n
2.5
2.0
§) 1.5
1.0
0.5
0.0
3.0 n
2.5
2.0
1.5
1.0
0.5
0.0
Daily: Low
10 15
Daily: Medium
20
10 15
Daily: High
10
15
Year
20
25
30
25
30
Figure A-l. Three different time series used for the model comparative runs.
Each time series has the same maximum value for each of the 30 separate
years, all other values were either multiplied by 0.5 (low), 1.0 (medium) or
1.5 (high).
27
-------
Appendix B: Explanation of Periodic Model and Density
Dependent Factor
A periodic matrix model is named as such because each cycle of the matrix (often
expressed as a year) is divided into several periods. Periods can be seasons, months,
etc.— and they do not have to be equal in length of time covered. Each period has its
own matrix model with its own demographic parameters (see Caswell 2001, Chapter
13). The matrix product of all sub-matrices within a cycle is the projection matrix for
that cycle, and its dominant eigenvalue is the annual growth rate for the population.
While the model can provide annual growth rates via the eigenvalue method, its use
herein is to following weekly population size.
Density dependence is incorporated into each sub-matrix using a density dependent
factor:
[Mj • nj • DDF = ni+l Eq B-l
Where:
Mj = a 13x13 sub-matrix for week / (/ ranges from 1 to 52),
n, = the population vector for week / (13 age classes),
DDF = a density dependent factor, see Equation B-2, and
n-,+1 = the population vector for week i+1.
The density dependent factor is calculated as (based on Leslie 1948):
DDF= - tr}. , N Eq B-2
n
Where:
A = maximum weekly growth rate in the absence of density dependence
(1.64),
n, = the total number of individuals in the ;'th week,
K = carrying capacity (set at 100).
The use of density dependence requires a carrying capacity. Since the maximum size
of field populations for mysids is not easily known, 100 was chosen as the
maximum— 100%, so the carrying capacity is a relative number. The weekly control
population growth rate (lambda) for a spawning week was 1.64, for a non-spawning
week, 0.49 19. As a population size approaches zero, the weekly growth rate
approaches the maximum rate. As a population size approaches the carrying
capacity, the maximum weekly growth rate approaches 1.0. If a population exceeds
the carrying capacity, the maximum weekly growth rate is less than 1.0.
19 Based on average control demographic parameters from Thursby (2009).
28
-------
Each year of exposure data is evaluated in 91-d running blocks—this is the expected
duration of Americamysis bahia life history (13 weeks). This is visualized in Figure B-l.
Each weekly sub-matrix is a 13x13 matrix, consisting of 13 age classes (Thursby 2009).
Age class 1 is only exposed to the concentrations represented in the current week.
Age class 2, however, has been exposed for 2 weeks—the current week and "reaches
back" into the previous week. That is, the current concentration within individuals
that are 2 weeks old is an integration of the previous 14 days. Each age class reaches
back a different number of days. Age class 13 reaches back the full 91 days. This
results in each age class having a different scaled internal concentration at any given
time since they are integrating over a different number of days. Each sub-matrix
covers a different portion of the annual exposure time series (Figure B-2).
Age Classes
8.0 -
7.0 -
c 6.0
o
1
1 5.0
§.3.0
x
1.0
40
Time (da)
60
80
Figure B-l. Demonstration of "reaching back" for each age class. The light
green vertical bar represents the current exposure week. The horizontal
black lines (solid and dashed) represent the time period to which each age
class has been exposed.
29
-------
-15 -10 -5
10 15 20 25 30 35 40 45 50 55
Exposure Week
Figure B-2. Picture of a single year's worth of exposure evaluation within the
model. Each sub-matrix (represented by a horizontal black bar) integrates
different, but overlapping, portions of the time series.
Figure B-3 demonstrates the cumulative internal concentration (represented by
Q(t)/kin~see Appendix C for a detailed explanation) for each age class for a
hypothetical 91-d exposure block. Survival for a given week is estimated as the lowest
calculated survival probability from among the 7 days represented by the current
week's exposure. The lowest value is used since individuals can only die once. This
survival value is the one used for a given age class within a given weekly sub-matrix.
For each year, 676 survival values are calculated—one for each age class (13) for each
week of the year.
30
-------
•Age 13
•Age 12
Age 11
•Age 10
•Age 9
Age 8
•Age 7
•Age 6
Age 5
•Age 4
Age 3
Age 2
Age 1
20
40 60
Time (da)
80
100
Figure B-3. Calculated Q(t)/kin values for each age class over the entire
length of their exposure. Note each curve does not begin a 0 because
the calculations begin with the end of day 1. As long as there is an in
water concentration greater that 0 on day 1 there will be a value for
Q(t)/kin greater than 1 on day 1.
Survivorship Calculations
The "hazard function" or "hazard rate" approach is used by the model to estimate
survival. Terms like "time-to-event" or "accelerated life testing" are also used in the
literature. In this model there are no sensitivity sub-groups; every individual has the
same probability of dying for any given exposure. Whether or not a particular
individual dies is a random process. For example, if during a particular time interval
the stressor exceeds some threshold, then a portion of the population will die. This
portion is determined by the "killing rate" or proportionality constant. Because the
sensitivities of the individuals that die are the same as the sensitivities of the
individual that do not die, there is no change in the range of sensitivity within the
population of living organisms. Thus, if during the next time interval the stressor level
does not change, then more individuals will die (the same proportion as the previous
time interval).
To understand how we apply the hazard function approach to a series of time-varying
exposures, we begin with raw data such as that shown in Figure B-4. This figure is
hypothetical survival data from a single toxicant concentration monitored for 10
periods (e.g., over 10 days).
31
-------
1.0
0.9-
0.8-
1 0.7-
>
1 0.5
0.3-
0.2-
0.1 -
0.0
3456
Time
Figure B-4. A hypothetical time to death curve for a single toxicant exposure
concentration.
To use the hazard function we need to calculate a new constant that Kooijman and
coworkers (e.g., Kooijman and Bedaux, 1996a,b) call the killing rate, /r^//. We fit a
survivorship curve to the above set of data that has kkm as the only unknown. The
development of this survivorship curve is explained below.
A good explanation for the relationship between the following is in Lee (1992). There
are three main functions that we need to understand:
S(t) = Survival Function—defined as the probability that an individual
survives longer than t. Also defined as l-F(t) where F(t) is the
cumulative density function (cdf) for mortality;
f(t) = Probability Density Function (pdf)—probability of dying in a given time
interval divided by the duration of that interval (also defined as the
derivative of F(t);
h(t) = Hazard Function—probability of dying at time t assuming that the
individual has survived to time t (also referred to as the age specific
mortality rate.
These values are interrelated.
h(t) =
/(O _/(0
\-F(f) S(f) I -cdf
Since the pdf\s the derivative of the cdf:
at
Eq B-3
Eq B-4
32
-------
By substitution into h(t) equation above we get:
S(t) dt
Integration from zero to t and using 5(0) = 1 , we have
Eq B-5
Eq B-6
Solving for S(t)
S(t) = exp
Eq B-7
Widianarko and Van Straalen (1996) use the term hazard rate, which is just the
probability of dying. They define hazard rate as being directly proportional to the
internal concentration of a toxicant. They introduce a proportionality constant 6, and
write the hazard rate as:
20
'(0
Eq B-8
Their 6 is similar to the killing rate /r/y// used by Kooijman and his coworkers. The two
values are related by:
- or, rearranging: 0 = k
m
Eq B-9
Using the above definition of 6, and Eq C-2 (from Appendix C) and B-6 from above we
get:
Eq B-10
in out
Which simplifies to:
Substituting this equation back into Eq B-4 we get:
S(t) = exp
Eq B-ll
Eq B-12
Others (e.g., Kooijman group) essentially use Q^-Q/vtc where Qwfc is the no effect concentration of
the toxicant—in other words the hazard rate is proportional to the degree to which the internal
concentration exceeds some no effect concentration.
33
-------
The solution to this integral is:
S(t) = exp
-kmCw(f-kout^
kout kout
Eq B-13
This equation is fit to the % survival vs. time data (e.g., Figure B-4) to determine /r/y//
which will be the only unknown in the above equation—Cw will be known (and
assumed a constant) for a particular data set. The constant kout is determined
independently (see Appendix C). Once we have an estimate for k^i, we can relate
survival to Q(t/kin.
If we just substitute the definition of 9 into Eq B-8 we have:
~k
h(t) =
Rearranging we get:
_°«!fr
"•I-
Q(f)
Eq B-14
Eq B-15
Applying Eq B-7:
S(t) =
Integrating:
S(t) = exp
''kill
Q(t}
dx)
Eq B-16
* out* kill
Q(t}
Eq B-17
Since the time step is 1 (a single day), we substitute t = 1 in to Eq B-17 and apply this
to the time series for Q(t)/kin (see Appendix D) and derive the time series for S(t).
Reproduction
Data sufficient to evaluate reproduction in the same manner as survival—that is,
enough information to establish an empirical rate constant for reproduction
analogous to kkm (e.g., krepro)—is seldom, if ever available. The model assumes a fixed
relationship (user defined) between survival and reproduction.
34
-------
Appendix C: Calculation of Elimination Constant
The equation used for estimating internal concentration assumes a one-
compartment model where the likelihood of death is associated with the
concentration of the toxicant in some critical compartment. One-compartment
models are assumed by most of the references on this topic. The rate of change of
the quantity of the toxicant in the organisms is defined by:
C-l
dt
Where Q^ = concentration of toxicant inside organism at time t (u,g/g);
Cw = concentration of toxicant in external medium (e.g., water—u,g/L);
k^ = uptake rate into organism (L/g-t) 21;
k0ut = elimination or detoxification rate (1/t);
t = time.
Simply put, this means that the amount that the internal concentration increases by
in a given time period is equal to a constant proportion of the external concentration
minus an elimination (or detoxification) rate that is proportional to the current
internal concentration. The solution to this differential is:
k
Eq C-2
k
out
Since kin and kout are considered constants for a given organism, the only variables
that influence Q^ in the above equation are the concentration of the toxicant in the
water and the duration of exposure to that concentration. If we assume Cw is a
constant, and that exposure is infinite, then the maximum internal concentration for
the toxicant is (k!n/kout)Cw22. Note that k!n/kout is the steady state bioaccumulation
factor.
Let's assume that we are exposing a group of individuals to a concentration that will
kill 50% of them. For this case, let's define x as the median threshold value. Since x is
the median internal threshold, this makes the external concentration that results in x
an LC50 value. If we rearrange Eq. C-2, then we can show:
EqC-3
21 Units of rate constants are somewhat arbitrary—established so that units of final value are correct.
22 Mathematically as t approaches oo, the "e" term in the above equation approaches 0.
35
-------
Which at t = oo, becomes:
\=^Q(x} EqC-4
kin
This, because we have defined Cw as an LC50, simplifies to:
Cw=tfLQw=LCW* EqC-5
11 in
Substituting back into Eq B-4 and rearranging, we can calculate the LC50 for any time
t.
EqC-6
In Figure C-l the water concentration to achieve a given toxic event is plotted vs.
time to that event. Although this can be any measure of toxicity, frequently we will
be plotting median lethality as the event (i.e., tracking the change in LC50 over time).
By fitting Eq B-6 to the data, we can estimate the asymptote (the LC50 at infinity— or
the incipient LC5023) and the value for kout— which is the shape parameter for the
curve, related to the steepness of the curve. Jager et al. (2011) make the point that
determination of the elimination constant using toxicity data is actually preferred to
basing that calculation on measured internal concentrations. Internal concentrations
may include toxicant not located at the primary action site, whereas, effects data are
always responding to the concentration at the action site— even if we do not know
the location of the specific site. Figure C-2 shows the application of Eq C-6 to data for
endosulfan LC50 values for Americamysis bahia. Data for days 1 through 4 are from
Schimmel (1981). Datum for 28 day exposure is from (McKenney 1982).
23 This term appears to have been introduced by Sprague, JB. 1969.
36
-------
Time
Figure C-l. Hypothetical relationship between time and LC50 values.
4.00 -
3.50 -
3.00 -
2.50 -
o>
3. 2.00
1.50 -
1.00 -
0.50
0.00
10 15 20
Time (da)
25
30
Figure C-2. Plot of mysid LC50 values for endosulfan. Curve was fit
using Eq C-6. The asymptotic LC50 is 0.80 ug/L and the kinetic
parameter kout is 0.25 t"1.
37
-------
Appendix D: Explanation of Mancini Calculation of Scaled
Internal Concentration
Toxicity is more closely related to internal toxicant concentration than external
concentration. But, internal concentrations are rarely measured as part of traditional
toxicity testing procedures. However, toxicity can be related to what has been
recently called a scaled internal concentration (Jager et al. 2011). Because of this, all
that is needed in order to estimate the internal kinetics of a given toxicant (using the
one-compartment model) is the external concentration time series and kout. The
latter is calculated using traditional toxicity data (see Appendix C), and the former is
from either direct measurements or, as in our case, modeled.
Even though Mancini (1983) was not the first to introduce the equations for the
single-compartment model, he seems to have been the first to show how we can
calculate internal values based on "environmental" time-varying concentrations.
Ultimately what is needed is an estimate of the % survival for a group of organisms
for each time period of interest. Ideally we would establish a relationship between
internal concentration and % survival (or calculate rates of uptake and elimination: k!n
and kout). However, that would require that we actually measure the internal
concentrations, which may be cumbersome and expensive. Alternately we can
establish a relationship between survival and some estimate of internal kinetics. In
the case of Mancini's paper that kinetic estimator is Q(t)/kjn. Recall from Appendix C
that the LC50ooar\d kout can be estimated using standard toxicity data. The LC50X is
the same as kourQ(x)/kjn. So if we divide the LC50oo by koutwe get Q(x/kjn. This means
survival can be related directly to Q(t)/kin, so only how Q(t)/kin varies with time is
needed and not
Mancini's equation 5 shows us how to calculate the time-varying quantity Q(t)/kjn. To
do this, the previous internal "quantity" loss is calculated for the current time
interval:
decay_of _previous_amount = —(—^--e k""''t Eq D-l
kin
And to this add the amount that this "quantity" would increase by from the current
external concentration—based on Eq C-2 (Appendix C):
k
new_amount = w \\-e °" Eq D-2
out
Where Cw(t) is the concentration in the external environment at time t. At the end of
the current time interval t the new value is:
EqD-3
38
-------
Using Eq D-3, the daily values for exposure concentrations are converted into daily
values for Q(t)/kjn. For a time step of 1 day Eq D-3 simplifies to:
EqD.4
in in
Note: The "t" in Q(t> Q(t-i) and Cw(t) is part of the value name and not a mathematical
value.
Appendix B shows survival can be expressed as a function of Q(f//r/n, therefore, once
the time series for Q(t/kin is known, the time series for survival is known.
39
-------
Appendix E: Mysid Chronic Data for Endosulfan and Model
Calibration (Quality Control)
In order to set the survival kinetic parameter (/CM), survival data overtime are needed
for the same concentration. The acute data available for endosulfan shown below
does not provide a lot of these data (Figure E-l). I could have just used the 3.02 ug/L
data—and that might have been okay. Instead, to increase the amount of data, the
dose response data for each day, along with the 28-d survival dose response data
(Figure E-2) were used to create a new data set of survival vs. time using four
different concentrations (1.5, 2.0, 2.5 and 3.0 ug/L) that more evenly cover the effect
range. These are plotted in Figures E-3 and E-4.
120 i
100
Endosulfan Acute Data
0 1 2Time(da)3 4 5
Figure E-l. Acute data for endosulfan from raw data available from Schimmel (1981).
40
-------
24 hr
48 hr
1.0
2.0 3.0 4.0
Endosulfan(ug/L)
5.0
Endosulphan (ug/L)
72 hr
96 hr
1.0
2.0 3.0 4.0
Endosulfan ug/L)
5.0
1.0
2.0 3.0 4.0
Endosulfan (ug/L)
5.0
28 da
2.0 3.0
Endosulfan (ug/L)
28 da Reproduction
1.0
2.0 3.0 4.0
Endosulfan (ug/L)
5.0
Figure E-2. Dose-response data for endosulphan survival (1, 2, 3 4 and 28 da) and 28-d reproduction dose
response.
41
-------
a.
o
3
V)
1.00
0.90
0.80
0.70
0.60
0.50
0.40
0.30
0.20
0.10
0.00
\ \ N\ ' N
\ \ \j
t|
• l
\ \
\ ', !i
I |
v ! »
i i
1 1
: \i .
^1<
.
\
\
\
|
0
1
Endosulfan (ug/L)
Figure E-3. Calculated survival dose-response data for endosulfan.
3,
Probabili
TO
»_
D
w
1.00
0.90
0.80
0.70
0.60
0.50
0.40
0.30
0.20
0.10
0.00
(
£•
A
• 1.5 ug/L
A2.0ug/L
B » 2.5 ug/L
• 3.0 ug/L
.
A
\J *
D 5 10 15 20 25 30
Time (days)
Figure E-4. Plot of data from Figure E-3 for 1.5, 2.0, 2.5 and 3.0 ug/L.
42
-------
The data from Figure E-4 were combined into a single regression (Figure E-5).
Equation B-17 was used to fit the data. The concentration (Cw) used was the average
for the four data sets (2.25 ug/L). The kinetic parameter for elimination (kout) was
from Equation B-6 fit to data in Figure C-2 (0.25 t ). An estimate of 0.60 t for
gave the best fit (r2 = 0.77).
*kill
10
20
25
30
15
Time (day)
Figure E-5. Plot of survivorship versus time using data from Figure 15. Solid red
line uses average concentration (2.25 ug/L) and all of the data (kk!// = 0.60 f1).
Dashed red line is fit using a kk!n = 1.00 1" .
Setting kinetic parameters for reproduction is not as straightforward as for survival.
This is primarily because there is not as much time-to-reproduction data available so
that Q(t)/kin can be related to reproduction kinetics. As a compromise, reproduction
was considered to be a function of the survival kinetics through the ratio of survival
LC50 to reproduction EC50 from the chronic test. This assumes that the short-term
exposure relationship of survival to reproduction is the same as that relationship in
the chronic data set. Figure E-2 (bottom 2 plots) shows that the chronic dose
responses of survival and reproduction are similar. The survival to reproduction ratio
is 0.89.
For calibration purposes, the model was run using constant concentrations covering
the range of measured values within the available 28-d chronic tests (McKenney
1982). A special density independent version of the model was created for
conducting the calibration24. This version of the model tracks population size;
The assumption is that chronic tests have sufficient food and low enough density of individuals to
make density dependent factors minimal.
43
-------
however, because only the first four weeks of the model are used (corresponding to a
chronic test's 28 days) there was no problem with the model output exceeding
Excel's capacity for number size. As with a chronic test, the model run began with all
individuals (100) assigned to the youngest age class. Calibration run only tracked this
initial cohort through each week as a measure of survival (consistent with a chronic
test). For reproduction, total number of young produced during the four weeks was
totaled separately.
The survival kinetic parameter of 0.60 was used as the starting point in the periodic
model. However, the model output (for the first 4 weeks) under estimated the effect
of endosulfan on mortality. Changing the kinetic parameter from 0.60 to 1.00
resulted in a better fit between the model and the laboratory data (see Figure E-6).
Note that a survivorship curve using 1.00 is not that different from the curve using
0.60 (Figure E-5).
Survival
120
100
2345
Endosulfan (ug/L)
Figure E-6. Comparison of model survival output with data from 28-d laboratory chronic
tests. Solid black points are the data from the 28-d endosulfan chronic tests. Solid red
line is the model output using the 1.00 kinetic parameter. Dashed red lines are plus and
minus one standard deviation from the model runs.
The model calibration runs used 0.89 as the survival-to-reproduction ratio. The
reproduction calibration results (number young per female) are shown in Figure E-7.
The match between test data and model run was reasonably good enough to keep
the 0.89 ratio for use in the final model runs for time-varying exposures.
44
-------
Reproduction
120
2 3
Endosulfan (ug/L)
Figure E-7. Comparison of model reproduction (based on number young per female)
output with data from 28-d laboratory chronic tests. Solid black points are the data from
the 28-d endosulfan chronic tests. Solid red line is the model output using the 0.89 SRR.
Dashed red lines are plus and minus one standard deviation from the 100 stochastic model
runs.
45
------- |